#  Curves with prescribed rational points

Katerina Santicola  Mathematics Institute
University of Warwick
CV4 7AL
United Kingdom<EMAIL_ADDRESS>
###### Abstract.

Given a smooth curve $C/\mathbb{Q}$ with genus $\geq 2$, we know by Faltings’
Theorem that $C(\mathbb{Q})$ is finite. Here we ask the reverse question:
given a finite set of rational points $S\subseteq\mathbb{P}^{n}(\mathbb{Q})$,
does there exist a smooth curve $C/\mathbb{Q}$ contained in $\mathbb{P}^{n}$
such that $C(\mathbb{Q})=S$? We answer this question in the affirmative by
providing an effective algorithm for constructing such a curve.

###### Key words and phrases:

Diophantine equations, rational points

###### 2010 Mathematics Subject Classification:

Primary 11D41, Secondary 11D61

The author’s research is funded by a doctoral studentship from the Heilbronn
Institute for Mathematical Research

##  1\. Introduction

The genus $g$ of a smooth projective curve $C/\mathbb{Q}$ impacts the
possibilities for its set of rational points. If $g=0$ then
$C(\mathbb{Q})=\emptyset$ or $C(\mathbb{Q})$ is infinite. If $g=1$ then
$C(\mathbb{Q})=\emptyset$, or by the Mordell–Weil theorem, the set
$C(\mathbb{Q})$ can be given the structure of a finitely generated abelian
group. Mazur’s torsion theorem [mazur] restricts the possible torsion
subgroups of elliptic curves. As a consequence, for example, there are no
genus $1$ curves $C/\mathbb{Q}$ with exactly $11$ rational points. For $g\geq
2$, Faltings’ theorem [faltings] asserts that $C(\mathbb{Q})$ is finite. For a
general curve $C/\mathbb{Q}$ of fixed genus $g\geq 2$, little is known about
the possible size of $C(\mathbb{Q})$ beyond finiteness. However, the famous
Uniformity Conjecture of Caporaso, Harris and Mazur [uniformity] predicts that
$\\#C(\mathbb{Q})$ is uniformly bounded as $C$ ranges through curves
$C/\mathbb{Q}$ of fixed genus $g\geq 2$. It is therefore natural to ask about
the possible restrictions on $\\#C(\mathbb{Q})$ as $C$ ranges through all
curves $C/\mathbb{Q}$ of all possible genera $g\geq 2$. In this paper, we show
that there are no restrictions on $\\#C(\mathbb{Q})$, or indeed on the set
$C(\mathbb{Q})$ itself, beyond finiteness.

######  Theorem 1.

Let $n\geq 2$. Given any finite set of rational points
$S\subset\mathbb{P}^{n}(\mathbb{Q})$, there exists a smooth curve
$C/\mathbb{Q}$ contained in $\mathbb{P}^{n}$, of genus $g\geq 2$, such that
$C(\mathbb{Q})=S$.

Our proof gives an effective algorithm for computing the curve $C$ from the
given set of rational points $S$. It builds on earlier constructions due to
Gajović [gajovic] and the current author [me], who respectively prove the
following theorems.

######  Theorem 2 (Gajović).

Let

$\mathscr{P}_{\mathbb{Z}}\;=\;\\{a^{n}\;:\;a\in\mathbb{Z},\quad n\geq 2\\}$

be the set of perfect powers in $\mathbb{Z}$. Let $S$ be a finite subset of
$\mathscr{P}_{\mathbb{Z}}$. Then there is a polynomial $f_{S}\in\mathbb{Z}[X]$
such that $f_{S}(\mathbb{Z})\cap\mathscr{P}_{\mathbb{Z}}=S$.

######  Theorem 3 (Santicola).

Let

$\mathscr{P}_{\mathbb{Q}}\;=\;\\{a^{n}\;:\;a\in\mathbb{Q},\quad n\geq 2\\}$

be the set of perfect powers in $\mathbb{Q}$. Let $S$ be a finite subset of
$\mathscr{P}_{\mathbb{Q}}$. Then there is a polynomial $f_{S}\in\mathbb{Z}[X]$
such that $f_{S}(\mathbb{Q})\cap\mathscr{P}_{\mathbb{Q}}=S$.

We will show how the construction of the polynomial $f$ in those papers can be
adapted to give a construction of a suitable curve $C$ in our context.

##  2\. Superelliptic Curves

The proof of Theorem 1 involves constructing suitable "superelliptic" curves,
and then gluing them together. By a superelliptic curve $C$ over $\mathbb{Q}$
[arithmetic] we mean a smooth projective curve associated to the affine plane
curve given by

$C:\;y^{m}=f(x)$

where $f$ is separable of degree $d\geq 3$ and $m\geq 2$ is an integer. If
$m=2$ and $d=3$ then $C$ is an elliptic curve, whereas if $m=2$ and $d\geq 5$
then $C$ is hyperelliptic. When $m=d$, the genus $g$ of the curve is
$g=(d-1)(d-2)/2$; thus if $d\geq 4$ then $g\geq 2$ and so $C(\mathbb{Q})$ is
finite by Faltings’ theorem [faltings].

######  Definition 1.

Let $S=\\{P_{1},\dots,P_{n}\\}\subset\mathbb{A}^{2}(\mathbb{Q})$ be a finite
set of distinct rational points in the affine plane. We say that $S$ is
acceptable if the following conditions are satisfied:

  1. (a)

if $x(P_{i})=x(P_{j})$ then $P_{i}=P_{j}$,

  2. (b)

$x(P_{j})\neq 0$ and $y(P_{j})\neq 0$ for all $j$,

  3. (c)

$x(P_{i})$, $y(P_{i})\in\mathbb{Z}$ for all $i$.

Here $x(P)$ and $y(P)$ denote the usual coordinates.

For any such $S$, we can find a suitable superelliptic curve $C$ whose set of
affine rational points is exactly $S$, and such that $C$ has no rational
points at infinity.

######  Theorem 4.

Given an acceptable set $S\subset\mathbb{A}^{2}(\mathbb{Q})$, and a finite
subset $\mathcal{B}\subset\mathbb{C}$, there exists a separable polynomial
$f_{S,\mathcal{B}}(X)\in\mathbb{Q}[X]$ of some degree $d\geq 4$, not vanishing
at any point in $\mathcal{B}$, such that $C(\mathbb{Q})=S$, where $C$ is the
smooth projective plane curve associated to the affine model

$C\;:\;y^{d}=f_{S,\mathcal{B}}(X).$

Now say we have a finite set of rational points in
$\mathbb{P}^{n}(\mathbb{Q})$. By finiteness of this set, we can always find a
hyperplane on which none of the points lie. After a suitable change of
coordinates, all the points lie in an affine chart. In the following Lemma, we
assert that we can find a suitable change of coordinates so that none of the
coordinates are zero, and all of the coordinates are distinct from each other.

######  Lemma 5.

Let $S=\\{P_{1},\dots,P_{t}\\}\subseteq\mathbb{A}^{n}(\mathbb{Q})$ be a finite
set of points, which are distinct and nonzero. Then there exists a change of
coordinates matrix $A\in\operatorname{GL}_{n}(\mathbb{Q})$ such that

  1. (a)

for all $i$, $j$, we have $x_{i}(AP_{j})\neq 0$,

  2. (b)

for all $i$, $j$, $k$, $\ell$ with $i\neq k$ or $j\neq\ell$, we have
$x_{i}(AP_{j})\neq x_{k}(AP_{\ell})$.

###### Proof.

We write $M_{n\times n}$ for the affine space of dimension $n^{2}$ considered
as the space of $n\times n$ matrices. Consider the following subvarieties of
$M_{n\times n}$:

  * •

$U=\\{A\in M_{n\times n}:\det(A)=0\\}$,

  * •

$V_{i,j}=\\{A\in M_{n\times n}:x_{i}(AP_{j})=0\\}$ with $1\leq i\leq n$ and
$1\leq j\leq t$,

  * •

$W_{i,j,k,\ell}=\\{A\in M_{n\times n}:x_{i}(AP_{j})=x_{k}(AP_{\ell})\\}$ with
$1\leq i,k\leq n$ and $1\leq j,\ell\leq t$ and either $i\neq k$ or
$j\neq\ell$.

It is easy to see that these are all proper subsets of $M_{n\times n}$ which
are closed in the Zariski topology. As the rational points on $M_{n\times n}$
are Zariski dense, it follows that there is some $A\in M_{n\times
n}(\mathbb{Q})$ that does not lie on any of the $U$, $V_{i,j}$,
$W_{i,j,k,\ell}$.

∎

Using Lemma 5, we are able to deduce Theorem 1 from Theorem 4.

###### Proof of Theorem 1.

We may assume from Lemma 5 that the set $S=\\{P_{1},\dotsc,P_{t}\\}$ satisfies
the conclusion of Lemma 5, scaling the points if necessary to ensure
integrality. For $2\leq j\leq n$, let

$S_{j}=\\{(x_{1}(P_{1}),x_{j}(P_{1}))\,,\,(x_{1}(P_{2}),x_{j}(P_{2}))\,,\,\dotsc\,,\,(x_{1}(P_{t}),x_{j}(P_{t}))\\}\;\subset\;\mathbb{A}^{2}(\mathbb{Q}).$

Observe that these are acceptable sets. We construct sequences of finite
subsets of $\mathbb{C}$:

$\mathscr{B}_{2}\subset\mathscr{B}_{3}\subset\cdots\subset\mathscr{B}_{n}$

and a sequence of separable polynomials
$f_{2},f_{3},\dotsc,f_{n}\in\mathbb{Q}[x_{1}]$ as follows. We start by letting
$\mathscr{B}_{2}=\emptyset$. We let $f_{2}=f_{S_{2},\mathscr{B}_{2}}$
constructed as in Theorem 4. For the inductive step of the construction, let
$\mathscr{B}_{j}$ be the union of $\mathscr{B}_{j-1}$ and the set of complex
roots of $f_{j-1}$, and we let $f_{j}=f_{S_{j},\mathscr{B}_{j}}$. Write
$d_{j}=\deg(f_{j})$. Let $C\subset\mathbb{P}^{n}$ be the projective closure of
the affine variety

$C\;=\;V\left(\,x_{2}^{d_{2}}-f_{2}(x_{1})\,,\,x_{3}^{d_{3}}-f_{3}(x_{1})\,,\,\cdots\,,\,x_{n}^{d_{n}}-f_{n}(x_{1})\,\right)\;\subset\mathbb{A}^{n}.$

Note that $C$ has function field

$\mathbb{Q}(C)=\mathbb{Q}\left(\sqrt[d_{2}]{f_{2}(X)},\sqrt[d_{3}]{f_{3}(X)},\dotsc,\sqrt[d_{n}]{f_{n}(X)}\right)$

which is a finite extension of $\mathbb{Q}(X)$. Thus $C$ is an absolutely
irreducible curve. It follows from Theorem 4 that
$C(\mathbb{Q})=\\{P_{1},P_{2},\dots,P_{t}\\}$. Moreover, we note that the
construction ensures that $f_{1},\dotsc,f_{n}$ have pairwise distinct complex
roots. It follows from this, using the Jacobian criterion, that the curve $C$
is smooth.

∎

The rest of this paper is devoted to the proof of Theorem 4.

##  3\. Rank zero elliptic curves

A consequence of Faltings’ theorem is the celebrated theorem of Darmon and
Granville [darmongranville]. They show that the generalised Fermat equation

$Ax^{p}+By^{q}=Cz^{r}$

where $A$, $B$, $C$, $p$, $q$, $r\in\mathbb{Z}$ and $1/p+1/q+1/r<1$ has
finitely many proper integer solutions. There is much progress toward showing
certain cases have no integer solutions, for more see [samirsurvey] and
[ivorrakraus]. To prove that the curve we construct whilst proving Theorem 4
has exactly the required set of rational points, we will need to first prove
that the following four equations:

(2)
$ax^{6}+by^{6}=cz^{3},\qquad(a,b,c)=(\ell,1,1),(\ell,1,2),(\ell,1,\ell),(\ell,1,2\ell)$

have no proper integer solutions when $\ell$ is a prime satisfying certain
conditions.

If $\ell$ is a prime such that 2 is not a cube residue modulo $\ell$, then the
equations $(\ell,1,2)$ and $(\ell,1,2\ell)$ have no nontrivial solutions. The
remaining two equations define genus 4 curves which admit maps to elliptic
curves. In particular, we look at

$\phi_{1}:C\longrightarrow E_{1}:y^{2}=x^{3}-a^{3}c^{2}b,$

and

$\phi_{2}:C\longrightarrow E_{2}:y^{2}=x^{3}-ab^{3}c^{2}.$

Then the equations $(\ell,1,1)$ and $(\ell,1,\ell)$ both define curves
admitting maps to

$E_{-\ell^{3}}:y^{2}=x^{3}-\ell^{3}.$

This is a quadratic twist of the Mordell curve $y^{2}=x^{3}-1$ by the prime
$\ell$. Corresponding to the 2-torsion point $(\ell,0)$ we have a 2-isogeny
$\phi:E_{\ell}\rightarrow E^{\prime}_{\ell}$, where

$E_{\ell}^{\prime}:y^{2}=x^{3}-6\ell x^{2}-3\ell^{2}x.$

Using the classic technique of 2-descent we have the resulting exact sequence

${0}$${E^{\prime}(K)/\phi(E(K))}$${\operatorname{Sel}^{(\phi)}(E/K)}$${{\Sha(E/K)[\phi]}}$${0}$

where $\operatorname{Sel}^{(\phi)}(E/K)$ is the $\phi-$Selmer group and
$\Sha(E/K)[\phi]$ is the $\phi-$torsion of the Tate-Shafarevich group. We have
a similar exact sequence for the dual isogeny $\hat{\phi}$. These exact
sequences allow us to bound the rank of $E(\mathbb{Q})$ in terms of the size
of the Selmer groups:

######  Lemma 6.

[phillipe] Let $E/\mathbb{Q}$ be an elliptic curve with a rational point of
order 2. Let $\phi:E\rightarrow E^{\prime}$ be an isogeny of degree $2$, with
$\hat{\phi}:E^{\prime}\rightarrow E$ the dual of $\phi$. Then

$\operatorname{rank}(E(\mathbb{Q}))\leq\dim_{\mathbb{F}_{2}}\operatorname{Sel}^{(\phi)}(E,\mathbb{Q})+\dim_{\mathbb{F}_{2}}\operatorname{Sel}^{(\hat{\phi})}(E^{\prime},\mathbb{Q})-2.$

Juyal, Moody and Roy [phillipe] use this Lemma to show that for the elliptic
curve $E_{p}:y^{2}=x^{3}-p^{3}$, the rank is 0 if $p\equiv 5\mod 12$, and
conditionally on the Parity Conjecture that the rank is 1 if $p\equiv 11\mod
12$. They also calculate the Selmer groups of $E_{p},E_{p}^{\prime}$ when
$p\equiv 1\mod 12$ if $p$ satisfies some additional conditions. We will choose
a prime $\ell$ that satisfies these additional conditions, so that we can
ensure the rank is 0.

Note if $\ell$ is a prime such that $\ell\equiv 1\mod 12$, then the Legendre
symbol

$\left(\frac{\sqrt{3}}{\ell}\right)$

is well defined, since it does not depend on the choice of sign for
$\sqrt{3}$.

######  Proposition 7.

If $\ell$ is a prime such that $\ell\equiv 1\mod 12$ and

$\left(\frac{2\sqrt{3}}{\ell}\right)=-1,$

then the elliptic curve

$E:y^{2}=x^{3}-\ell^{3}$

is such that

$E(\mathbb{Q})=\\{(\ell,0),\infty\\}.$

###### Proof.

Let $\phi:E\rightarrow E^{\prime}$ be the 2-isogeny corresponding to the
2-torsion point $(\ell,0)$, and $\hat{\phi}:E^{\prime}\rightarrow E$ its dual.
From Proposition 3 in [phillipe], we have

$\operatorname{Sel}^{(\phi)}(E,\mathbb{Q})=\langle-3\rangle\text{ if
}\ell\equiv 1(12)\text{ and }\left(\frac{3+2\sqrt{3}}{\ell}\right)=-1.$

For the dual isogeny, by Proposition 2 in [phillipe] we have

$\operatorname{Sel}^{(\hat{\phi})}(E^{\prime},\mathbb{Q})=\langle
3\rangle\text{ if }\ell\equiv 1(12)\text{ and
}\left(\frac{-6-2\sqrt{-3}}{\ell}\right)=-1.$

Both of these conditions are equivalent to

$\ell\equiv 1(12)\text{ and }\left(\frac{2\sqrt{3}}{\ell}\right)=-1.$

This follows from the
identities111https://math.stackexchange.com/questions/4736110

$2\sqrt{3}\left(3+2\sqrt{3}\right)=\left(\sqrt{3}\right)^{2}\left(1+\sqrt{3}\right)^{2},$

and

$2\sqrt{3}\left(-6+2\sqrt{-3}\right)=\left(\sqrt{-1}\right)^{2}\left(\sqrt{3}\right)^{2}\left(1+\sqrt{-1}\right)^{2}\left(1+\sqrt{-3}\right)^{2}.$

By Lemma 6 we have that

$\operatorname{rank}(E(\mathbb{Q}))\leq 1+1-2=0$

as required. The torsion structure of $E$ is well known [silverman], and is
$\mathbb{Z}/2\mathbb{Z}=\\{(\ell,0),\infty\\}$. ∎

As a consequence, the equations (2) have no proper integer solutions when
$\ell$ is a prime such that $\ell\equiv 1\mod 12$ and $X^{4}-12$ factors into
two irreducible quadratic polynomials modulo $\ell$.

##  4\. Proof of Theorem 4

We start by proving the following Lemma.

######  Lemma 8.

Let $S=\\{(a_{1},b_{1}),\dots,(a_{r},b_{r})\\}$ be an acceptable set of points
of size $|S|=r$. Let $d$, $n$, $m$ be positive integers and $n\geq r$. Then
there exists a separable polynomial

$h(X)=\sum_{i=1}^{n}\alpha_{i}X^{i}\in\mathbb{Q}[X]$

satisfying the following conditions

  1. (i)

$h(a_{i})=b_{i}^{d}$ for all $(a_{i},b_{i})\in S$;

  2. (ii)

$h(X)$ has degree $n$

  3. (iii)

$m\cdot h(X)\in\mathbb{Z}[X]$ where

$m=\prod_{j<k}^{r}(a_{j}-a_{k}).$

###### Proof of Lemma 8.

Let $L(X)$ be the Lagrange interpolation polynomial passing through the points
$S^{\prime}$, where

$S^{\prime}=\\{(a_{i},b_{i}^{d})\mid(a_{i},b_{i})\in S\\},$

so we have

$L(X)=\sum_{j=1}^{r}b_{j}^{d}\prod_{\begin{subarray}{c}0\leq k\leq r\\\ k\neq
j\end{subarray}}\frac{X-a_{k}}{a_{j}-a_{k}}.$

Note that $m\cdot L(X)\in\mathbb{Z}[X]$. Now let

$b(X)=\prod_{\begin{subarray}{c}i=1\end{subarray}}^{r}(X-a_{i}).$

Since $S$ is acceptable, $b_{j}\neq 0$ for all $j$, so $L(X)$ and $b(X)$ are
coprime. Choose $c(X)\in\mathbb{Z}[X]$ to be coprime to $L(X)$, of degree
$n-r$. Define

$R(X,Y)=L(X)+b(X)c(X)Y\in\mathbb{Q}(Y)[X].$

This is irreducible in $\mathbb{Q}(Y)[X]$. By Hilbert’s irreducibility theorem
[serre], there is a thin set $A\subset\mathbb{P}^{1}(\mathbb{Q})$ such that if
$t\notin A$, then $R(X,t)$ is irreducible over $\mathbb{Q}$. Let $N$ be a
positive integer satisfying $Nb(X)c(X)\in\mathbb{Z}[X]$. As $A$ is thin, it
follows that there exist infinitely many integers $y\notin A^{\prime}$
satisfying $N\mid y$. Choose any such $y$ and let $h(X)=R(X,y)$. As $N\mid y$
we see that $yb(X)c(X)\in\mathbb{Z}[X]$, and so $mh(X)\in\mathbb{Z}[X]$.

∎

For an acceptable set
$S=\\{(a_{1},b_{1}),\dots,(a_{r},b_{r})\\}\subseteq\mathbb{A}^{2}(\mathbb{Q})$,
Lemma 8 constructs a separable polynomial such that $h(a_{i})=b_{i}^{d}$ for
all $(a_{i},b_{i})\in S$. We now have to modify the polynomial to ensure that
the resulting curve has no other rational points.

Our general argument will not apply to the case $|S|=1$. In this case, we can
take any curve $C/\mathbb{Q}$ of genus $g\geq 2$ with one rational point. A
suitable change of coordinates will send the point in $S$ to the point on $C$.

From now on we assume that $S$ is acceptable and $|S|=r>1$. Let

(3) $m=\prod_{j<k}^{r}(a_{j}-a_{k}).$

Let $p_{1},...,p_{t}$ be the distinct prime divisors of $m$, and let

$k=\text{LCM}(6,p_{1}-1,....,p_{t}-1).$

Now set

$n=2r+e,\qquad d=2kr+n$

where $e=0,2$ or $1$ depending on whether $2r(k+1)$ is $0,1$ or $2\mod 3$
respectively. This guarantees that $d\equiv 0\mod 3$. Now apply Lemma 8, with
$d,n,m$ as above. Observe that $mh\in\mathbb{Z}[X]$ from the construction in
the proof of Lemma 8. Let

$H(X):=\prod_{i=1}^{r}(X-a_{i})^{k}(h(X)-1).$

In order to prove that the curve we will construct is smooth and that its set
of rational points is exactly $S$, we will choose a prime $\ell$ that
satisfies a list of conditions. These conditions are luckily all compatible,
insofar that there is a positive density of primes that satisfies them.

######  Lemma 9.

Let $N\in\mathbb{Z}_{>0}$ such that $12\mid N$. Then there exists a positive
density of primes $p$ such that:

  1. (1)

$p\equiv 1\mod N$,

  2. (2)

$X^{3}-2$ is irreducible modulo $p$,

  3. (3)

$X^{4}-12$ factors into two irreducible quadratic polynomials modulo $p$.

###### Proof.

Let $K_{1}=\mathbb{Q}(\zeta_{N})$, let $K_{2}$ be the splitting field of
$X^{3}-2$ over $\mathbb{Q}$, and let $K_{3}$ be the splitting field of
$X^{4}-12$ over $\mathbb{Q}$.

We first compute that $\operatorname{Gal}(K_{2}K_{3}/\mathbb{Q})\cong
C_{3}\rtimes D_{4}$ in Sage [sage] by comparing character tables. Given that
$\mathbb{Q}(\zeta_{12})\subset K_{1}$ and $\mathbb{Q}(\zeta_{12})\subset
K_{2}$, we have that $\mathbb{Q}(\zeta_{12})\subset K_{1}\cap K_{2}K_{3}$.
Then since $(C_{3}\rtimes D_{4})^{\text{ab}}\cong C_{2}^{2}$, and $K_{1}$ is
an abelian extension, we have that $K_{1}\cap
K_{2}K_{3}=\mathbb{Q}(\zeta_{12})$. So, we have that

$\operatorname{Gal}(K_{1}K_{2}K_{3}/\mathbb{Q})\cong\\{(\sigma,\tau)\in\operatorname{Gal}(K_{1}/\mathbb{Q})\times\operatorname{Gal}(K_{2}K_{3}/\mathbb{Q})\mid\sigma(\zeta_{12})=\tau(\zeta_{12})\\}$

Now consider $(\sigma,\tau)$ where $\sigma=id$ and
$\tau=\operatorname{Frob}_{13}$. Note that $p=13$ satisfies (2) and (3). Since
$\sigma,\tau$ both act trivially on $\zeta_{12}$, we have
$(\sigma,\tau)\in\operatorname{Gal}(K_{1}K_{2}K_{3}/\mathbb{Q})$. By
Chebotarev’s Theorem [neukirch] there is then a positive density of primes
such that $\operatorname{Frob}_{p}=\tau$.

∎

By Lemma 9 we may choose an odd prime $\ell$ such that $h$ is separable modulo
$\ell$, and

  * •

$\ell\equiv 1\pmod{m}$, where $m$ is (3),

  * •

$2$ is not a cube reside mod $\ell$,

  * •

$X^{4}-12$ factors into two irreducible quadratic polynomials modulo $\ell$.

  * •

$\ell$ does not divide the numerator or denominator of any coefficient of
$h(X)$, nor of $H(X)$.

We let

(4) $g(X)=\ell q^{k}\prod_{i=1}^{r}(X-a_{i})^{k}+1,\qquad
f_{S}(X)=g(X)((h(X)-1)g(X)+1)$

where $q\neq 2$, $\ell$ is a prime, which we will choose shortly. Let

(5) $C:y^{d}=f_{S}(X).$

We note that

$f_{S}(a_{i})=h(a_{i})=b_{i}^{d}$

as required. Moreover,

$\deg f_{S}(X)=2\deg g(X)+\deg h(X)=2kr+n=d.$

Note that the leading coefficient of $f_{S}(X)$ is

$a_{d}=\ell^{2}q^{2k}\alpha_{n}$

where $\alpha_{n}$ is the leading coefficient of $h(X)$, and by the choice of
$\ell$ in Lemma 8, $v_{\ell}(a_{d})=2$. There are therefore no rational points
at infinity, as $\sqrt[d]{a_{d}}\notin\mathbb{Q}$.

We want the zeros of $f_{S}(X)$ to lie outside the given finite set
$\mathcal{B}\subset\mathbb{C}$. If $x$ is a root of $f_{S}(X)$, then $x\neq
a_{i}$ (since $f_{S}(a_{i})=b_{i}^{d}\neq 0$), and any $x\neq a_{i}$ is a root
of $f_{S}(X)$ for finitely many values of $q$. Thus we may choose a prime $q$
so that the zeros of $f_{S}(X)$ lie outside $\mathcal{B}$.

To check the curve is smooth, we need to show $f_{S}(X)$ is separable. We are
able to show this through the theory of Newton polygons.

######  Lemma 10.

$f_{S}(X)\in\mathbb{Q}[X]$ is separable.

###### Proof.

Note that $g(X)$ and $(h(X)-1)g(X)+1$ are coprime, thus it is enough to show
that each is separable. Applying Eisenstein’s criterion with prime $\ell$ to
the reverse polynomial of $g(X)$, shows that $g(X)$ is irreducible.

We now show that $(h(X)-1)g(X)+1$ is separable. We may write this as

(6) $(h(X)-1)g(X)+1\;=\;h(X)\;+\;\ell
q^{k}H(X)\;=\;\sum_{i=0}^{kr+n}c_{i}X^{i}.$

Recall that $\ell$ was chosen so that it does not divide the prime appearing
in the numerator or denominator of $h$, nor of $H$. Therefore

$v_{\ell}(c_{i})=\begin{cases}0&0\leq i\leq n,\\\ 1&n+1\leq i\leq
kr+n.\end{cases}$

Thus the Newton polygon of the polynomial (6) consists of two segments of
lengths $n$ and $kr$ and slopes $0$ and $1/kr$ respectively. By the theory of
Newton polygons, [Gouvea]*7.2.3 (6) will factor as the product of two
irreducible polynomials in $\mathbb{Z}_{\ell}[X]$:

$(h(X)-1)g(X)+1\;=\;F_{1}(X)F_{2}(X)$

corresponding to the two line segments of the Newton polygon. Here $F_{1}(X)$
has degree $n$ and its roots have valuation $0$, and $F_{2}(X)$ has degree
$kr$ and its roots have valuation $1/kr$. In particular, $F_{1}$, $F_{2}$ are
coprime. Moreover, $F_{1}\equiv h(X)\pmod{\ell}$. However $h(X)$ is
irreducible modulo $\ell$ by our choice of $\ell$. Thus $F_{1}$ is irreducible
in $\mathbb{Z}_{\ell}[X]$. The second factor $F_{2}$ is irreducible because it
is pure of slope $1/{rk}$. This means their product $F_{1}(X)F_{2}(X)$ is
separable. ∎

Now that we have shown the curve is smooth, we want to show its set of
rational points is exactly $S$. So, we want to show that if $f_{S}(x)=y^{d}$
for some rational $y\in\mathbb{Q}$, then $(x,y)\in S$. Write

(7) $x=\frac{u}{v},\qquad u,~{}v\in\mathbb{Z},\quad v\geq 1,\quad\gcd(u,v)=1.$

Let

(8) $A\;=\;q\prod_{i=1}^{r}(u-a_{i}v),\qquad
B\;=\;\frac{A}{\gcd(A,v^{r})},\qquad w\;=\;\frac{v^{r}}{\gcd(A,v^{r})}.$

Note that

(9) $\gcd(B,w)=1,\qquad w\geq 1.$

Assuming $r=|S|>1$, we note that

  * •

either $q\mid v$, in which case $q\,||\,A$ and hence $q\nmid B$, $q\mid w$,

  * •

or $q\nmid v$, in which case $q\mid B$ and $q\nmid w$.

Recall that we would like to show that $x=u/v=a_{i}$ for some $1\leq i\leq r$.
This is equivalent to $A=0$, which is equivalent to $B=0$. Thus our proof will
be complete on showing that $B=0$.

We note that

(10) $g(x)\;=\;\frac{\ell A^{k}+v^{kr}}{v^{kr}}\;=\;\frac{\ell
B^{k}+w^{k}}{w^{k}}.$

Now $h(X)$ is a polynomial of degree $n$ as in Lemma 8. So we have

$h(x)=\frac{M}{mv^{n}}$

where $M$ is some integer and $m$ is as above. So

(11) $(h(x)-1)g(x)+1\;=\;\frac{(M-mv^{n})(\ell
B^{k}+w^{k})+mv^{n}w^{k}}{mv^{n}w^{k}}.$

Recall that $f_{S}(x)=g(x)\cdot((h(x)-1)g(x)+1)$. We wish to compare the prime
divisors of these two factors.

######  Lemma 11.

$\gcd(\ell B^{k}+w^{k},mv)\mid 2\ell$.

###### Proof.

Let $p$ be a prime, not equal to 2 or $\ell$, dividing both $\ell B^{k}+w^{k}$
and $v$. From (8) we see that $\ell B^{k}+w^{k}$ divides $\ell A^{k}+v^{kr}$.
Thus $p\mid(\ell A^{k}+v^{kr})$ and so $p\mid A$. From the product definition
of $A$ in (8), and since $\gcd(u,v)=1$ from (7), we must have that $p=q$.
Since $r>1$ and $q\mid v$, we must have that $q\mid w$, and so $q\mid B$, a
contradiction.

Now let $p$ divide both $\ell B^{k}+w^{k}$ and $m$. Since $\gcd(B,w)=1$ (9),
we have that $p\nmid w$ and $p\nmid B$, so $B,w\in\mathbb{Z}_{p}^{*}$. Since
we have chosen $k$ so that $p-1\mid k$, and $\ell$ so that $\ell\equiv 1\mod
m$, Fermat’s little theorem gives $\ell B^{k}+w^{k}\equiv\ell+1\equiv 2\mod
p$. Since $p$ is odd, this is a contradiction.

Next, note that since $\gcd(B,w)=1$ and $k$ is even, if $\ell
B^{k}+w^{k}\equiv 0\mod 2$, then $\ell B^{k}+w^{k}\equiv 6\mod 8$, so
$4\nmid\ell B^{k}+w^{k}$. Similarly if $\ell B^{k}+w^{k}\equiv 0\mod\ell$,
then $\ell\mid w$, and so if $\ell B^{k}+w^{k}\equiv 0\mod\ell^{2}$, we would
have that $\ell\mid B$, which contradicts coprimality.

∎

######  Lemma 12.

If $f_{S}(x)=0$ then $B=0$.

###### Proof of 12.

Note that $g(x)=0$ is not possible since $g(x)>0$. So if

$f_{S}(x)=g(x)((h(x)-1)g(x)+1)=0$

we must have that $(h(x)-1)g(x)+1=0$, in which case

$(\ell B^{k}+w^{k})\mid-mv^{n}w^{k}$

From (11), if $p\mid mv^{n}w^{k}$ and $p\mid\ell B^{k}+w^{k}$ then $p=2$ or
$p=\ell$, so $\ell B^{k}+w^{k}\mid 2\ell$. We must have $w=1$ and $B=0$. ∎

###### Proof of 4.

Recall that we have supposed that $f_{S}(x)=y^{d}$ for some rational $y$, and
to complete the proof of Theorem 4, it is enough to show that $B=0$. Lemma 12
establishes this if $f_{S}(x)=0$. Thus we may suppose $f_{S}(x)\neq 0$. Let
$p$ be a prime, $p\neq 2,\ell$, dividing $\ell B^{k}+w^{k}$. Then, as before
$p\nmid w$ since $\gcd(B,w)=1$, and $p\nmid v$ from Lemma 11. From the
expression for $(h(x)-1)g(x)+1$ in (11) we see that $p$ divides neither the
numerator nor the denominator, and hence
$\operatorname{ord}_{p}((h(x)-1)g(x)+1)=0$. However
$g(x)((h(x)-1)g(x)+1)=f_{S}(x)=y^{d}$, so

$\operatorname{ord}_{p}(\ell
B^{k}+w^{k})\;=\;\operatorname{ord}_{p}\left(\frac{\ell
B^{k}+w^{k}}{w^{k}}\right)\;=\;\operatorname{ord}_{p}(g(x))\;=\;\operatorname{ord}_{p}(y^{d})\;\equiv\;0\pmod{d}.$

As this is true for every odd prime dividing $\ell B^{k}+w^{k}$ we have $\ell
B^{k}+w^{k}=2^{i}\ell^{j}z^{d}$ where $i,j\in\\{0,1\\}$.

Since $6\mid k$, $3\mid d$, and $\ell$ satisfies conditions (2) and (3) in
Lemma 9, then by Section 3, these equations have no nontrivial solutions. This
means $B=0$, as required.

∎

## Acknowledgements

The author would like to thank Samir Siksek for his continuous support, as
well as Ross Paterson and Michael Stoll for helpful conversations.

## References


# Modified Poisson-Boltzmann equations and macroscopic forces in inhomogeneous
ionic fluids

Yury A. Budkov<EMAIL_ADDRESS>School of Applied Mathematics, HSE University,
Tallinskaya st. 34, 123458 Moscow, Russia  G.A. Krestov Institute of Solution
Chemistry of the Russian Academy of Sciences, 153045, Akademicheskaya st. 1,
Ivanovo, Russia  Andrei L. Kolesnikov<EMAIL_ADDRESS>Institut
für Nichtklassische Chemie e.V., Permoserstr. 15, 04318 Leipzig, Germany

###### Abstract

We propose a field-theoretical approach based on the thermodynamic
perturbation theory and within it derive a grand thermodynamic potential of
the inhomogeneous ionic fluid as a functional of electrostatic potential for
an arbitrary reference fluid system. We obtain a modified Poisson-Boltzmann
equation as the Euler-Lagrange equation for the obtained functional. Applying
Noether’s theorem to this functional, we derive a general mean-field
expression for the stress tensor consistent with the respective modified
Poisson-Boltzmann equation. We derive a general expression for the macroscopic
force acting on the dielectric or conductive body immersed in an ionic fluid.
In particular, we derive a general mean-field expression for the disjoining
pressure of an ionic fluid in a slit pore. We apply the developed formalism to
describe three ionic fluid models of practical importance: nonpolarizable
models (including the well-known Poisson-Boltzmann and Poisson-Fermi
equations), polarizable models (ions carry nonzero permanent dipole or static
polarizability), and models of ion-dipole mixtures (including the well-known
Poisson-Boltzmann-Langevin equation). For these models, we obtain modified
Poisson-Boltzmann equations and respective stress tensors, which could be
valuable for different applications, where it is necessary to estimate the
macroscopic forces acting on the dielectric or conductive bodies (electrodes,
colloids, membranes, etc.) together with the local electrostatic potential
(field) and ionic concentrations.

##  I Introduction

Ionic fluids (plasma, electrolyte solutions, molten salts, and room
temperature ionic liquids) have recently attracted a lot of researchers’ and
chemical engineers’ attention. This was mostly caused by the use of ionic
fluids in various applications, including lipid and ion-exchange membranes,
biomacromolecules, colloids, batteries, fuel cells, supercapacitors, etc. In
all these examples, ionic fluids, which either interact with the charged
surface of the membrane, macromolecule, colloid, or electrode surface or are
confined in charged nanopores, are strongly inhomogeneous. Due to the
inhomogeneity of the ionic fluid, the local electrical neutrality is violated,
so that its description requires numerical solution of one or another self-
consistent field equation for the electrostatic potential with appropriate
boundary conditions: the classical Poisson-Boltzmann (PB) equation or its
generalizations. The latter are known in the literature as modified PB
equations.

The modified PB equations obtained up to now take into account the individual
properties of ionic species and the chemical specifics of the surface
interacting with them. In this case it is necessary to especially emphasize
the consideration of the non-zero size of ions (steric interactions) Kralj-
Iglič and Iglič (1996); Borukhov, Andelman, and Orland (1997); Maggs and
Podgornik (2016); Kornyshev (2007), their specific and structural interactions
Goodwin, Feng, and Kornyshev (2017); Blossey, Maggs, and Podgornik (2017),
polarizability/quadrupolarizability and permanent dipole moment of ionic
species and solvent (co-solvent) molecules Frydel (2011); Budkov  _et al._
(2020); Budkov, Zavarzin, and Kolesnikov (2021); Budkov, Kolesnikov, and
Kiselev (2015, 2016); Budkov  _et al._ (2018); Coalson and Duncan (1996);
Abrashkin, Andelman, and Orland (2007); Iglič, Gongadze, and Bohinc (2010);
Buyukdagli and Ala-Nissila (2013); Slavchov (2014), surface charge regulaton
Podgornik (2018), specific ion-electrode interactions Budkov  _et al._ (2018);
Uematsu, Netz, and Bonthuis (2018) and electrostatic correlations Buyukdagli
(2020); de Souza and Bazant (2020); Bazant, Storey, and Kornyshev (2011).
However, despite the variety of the modified PB equations, to the best of our
knowledge, there is still no first-principles approaches allowing us to derive
them from the unified standpoint and, simultaneously, take into account the
aforementioned effects. Moreover, for applications, where it is required to
calculate the macroscopic forces acting on conductive or dielectric bodies
immersed in the ionic fluid (such as electrodes, colloid particles, membranes,
micelles, dust particles, etc.), it is necessary to have a stress tensor
consistent with the modified PB equation.

In this paper, we will present such a field-theoretical approach based on the
thermodynamic perturbation theory allowing us to derive the grand
thermodynamic potential of an inhomogeneous ionic fluid as a functional of the
electrostatic potential for an arbitrary reference system, modified PB
equation (as the Euler-Lagrange equation for this functional), and a stress
tensor (with use of Noether’s theorem). We will apply this approach to
derivation of the modified PB equations and stress tensors for three types of
models that can be employed in the applications: (1) unpolarizable ionic fluid
models, (2) polarizable ionic fluid models, and (3) polarizable models of the
ion-dipole mixture (solvent in salt models).

##  II Mean-field theory of ionic fluids

###  II.1 Model I. Nonpolarizable ionic fluid mean-field models

Let us consider an ionic fluid – a set of two types of ions with electric
charges $q_{\pm}$ confined in a volume $V$ at a temperature $T$ contacting
with an infinite reservoir with ions (bulk). We assume that the ions are
immersed in a dielectric background of constant permittivity,
$\varepsilon_{b}$. The ions interact with each other, in addition to the
Coulomb interactions, via the short-range repulsive potentials
$U_{\alpha\gamma}(\mathbb{r})$, where indexes $\alpha,\gamma=\pm$ enumerate
the ions types. Thus, the grand partition function can be written as follows

$\Xi=\sum\limits_{n_{+}=0}^{\infty}\sum\limits_{n_{-}=0}^{\infty}\frac{z_{+}^{n_{+}}z_{-}^{n_{-}}}{n_{+}!n_{-}!}Q_{n_{+}n_{-}},$
(1)

with the fugacities $z_{\pm}=\Lambda_{\pm}^{-3}e^{\beta\mu_{\pm}}\theta_{\pm}$
(expressed via the chemical potentials, $\mu_{\pm}$, the de Broglie thermal
wavelengths, $\Lambda_{\pm}$, and internal partition functions,
$\theta_{\pm}$, of the ions) of the ionic species and configuration integrals

$Q_{n_{+}n_{-}}=\int d\Gamma_{n_{+}n_{-}}\exp\left[-\beta
U_{n_{+}n_{-}}\right],$ (2)

where $\beta=(k_{B}T)^{-1}$ ($k_{B}$ is the Boltzmann constant, $T$ is the
temperature) and

$\int
d\Gamma_{n_{+}n_{-}}(\cdot)=\int..\int\prod\limits_{j=1}^{n_{+}}{d\mathbb{r}_{j}^{(+)}}\prod\limits_{k=1}^{n_{-}}{d\mathbb{r}_{k}^{(-)}}(\cdot).$
(3)

The total potential energy, $U_{n_{+}n_{-}}$, of the interionic interactions
can be written as a sum of several basic terms:

$U_{n_{+}n_{-}}=U_{n_{+}n_{-}}^{(el)}+U_{n_{+}n_{-}}^{(ev)}+U_{n_{+}n_{-}}^{(ext)}-U_{self}.$
(4)

The first term in eq. (4) describes the total potential energy of
electrostatic interactions

$U_{n_{+}n_{-}}^{(el)}=\frac{1}{2}\int d\mathbb{r}\int
d\mathbb{r}^{\prime}\hat{\rho}(\mathbb{r})G(\mathbb{r}-\mathbb{r}^{\prime})\hat{\rho}(\mathbb{r}^{\prime})=\frac{1}{2}\left(\hat{\rho}G\hat{\rho}\right),$
(5)

where
$G(\mathbb{r}-\mathbb{r}^{\prime})=1/4\pi\varepsilon_{b}|\mathbb{r}-\mathbb{r}^{\prime}|$
is the Green function of the Poisson equation;

$\hat{\rho}(\mathbb{r})=\hat{\rho}_{i}(\mathbb{r})+{\rho}_{ext}(\mathbb{r}),$
(6)

is the local charge density containing microscopic charge density of the ions

$\hat{\rho}_{i}(\mathbb{r})=\sum\limits_{\alpha=\pm}q_{\alpha}\hat{c}_{\alpha}(\mathbb{r})$
(7)

with the microscopic ionic concentrations

$\hat{c}_{\alpha}(\mathbb{r})=\sum_{j=1}^{n_{\alpha}}\delta\left(\mathbb{r}-\mathbb{r}_{j}^{(\alpha)}\right)$
(8)

and the local charge density ${\rho}_{ext}(\mathbb{r})$ of the external
charges. The second term in eq. (4) describes the total potential energy of
the excluded volume interactions, i.e.

$U_{n_{+}n_{-}}^{(ev)}=\frac{1}{2}\sum\limits_{\alpha\gamma}\left(\hat{c}_{\alpha}U_{\alpha\gamma}\hat{c}_{\gamma}\right)=\frac{1}{2}\int
d\mathbb{r}\int
d\mathbb{r}^{\prime}\sum\limits_{\alpha\gamma}\hat{c}_{\alpha}(\mathbb{r})U_{\alpha\gamma}(\mathbb{r}-\mathbb{r}^{\prime})\hat{c}_{\gamma}(\mathbb{r}^{\prime}).$
(9)

The third term describes the interactions of ions with external fields with
potential energies, $u_{\pm}(\mathbb{r})$:

$U_{n_{+}n_{-}}^{(ext)}=\sum\limits_{\alpha}\int
d\mathbb{r}\hat{c}_{\alpha}(\mathbb{r})u_{\alpha}(\mathbb{r})=\sum\limits_{\alpha}(\hat{c}_{\alpha}u_{\alpha}).$
(10)

The last term in eq. (4) is the self-energy of the ions determined by the
following equation
$U_{self}=\left({n_{+}q_{+}^{2}}+{n_{-}q_{-}^{2}}\right)G(0)/2+({n_{+}U_{++}(0)}+{n_{-}U_{--}(0)})/2$.

Futher, using the standard Hubbard-Stratonovich transformations

$e^{-\frac{\beta}{2}(\hat{\rho}G\hat{\rho})}=\int\frac{\mathcal{D}\varphi}{C}e^{-\frac{\beta}{2}(\varphi
G^{-1}\varphi)+i\beta(\hat{\rho}\varphi)},$ (11)
$e^{-\frac{\beta}{2}\sum\limits_{\alpha\gamma}(\hat{c}_{\alpha}U_{\alpha\gamma}\hat{c}_{\gamma})}=\int\frac{\mathcal{D}\Phi}{\mathcal{N}_{U}}e^{-\frac{\beta}{2}\sum\limits_{\alpha\gamma}(\Phi_{\alpha}U^{-1}_{\alpha\gamma}\Phi_{\gamma})+i\beta\sum\limits_{\alpha}(\hat{c}_{\alpha}\Phi_{\alpha})},$
(12)

after some algebra, we obtain the following functional integral representation
of the grand partition function

$\Xi=\int\frac{\mathcal{D}\varphi}{C}e^{-\frac{\beta}{2}\left(\varphi
G^{-1}\varphi\right)+i\beta(\rho_{ext}\varphi)}\Xi_{R}[\varphi],$ (13)

where $C=\int\mathcal{D}\varphi e^{-\frac{\beta}{2}\left(\varphi
G^{-1}\varphi\right)}$ and
$\mathcal{N}_{U}=\int{\mathcal{D}\Phi}e^{-\frac{\beta}{2}\sum\limits_{\alpha\gamma}(\Phi_{\alpha}U^{-1}_{\alpha\gamma}\Phi_{\gamma})}$
are the normalization constants of the Gaussian measures. The inverse
electrostatic Green function has the following form

$G^{-1}(\mathbb{r}-\mathbb{r}^{\prime})=-\varepsilon_{b}\nabla^{2}\delta(\mathbb{r}-\mathbb{r}^{\prime}).$
(14)

The grand partition function, $\Xi_{R}$, of the reference fluid with pure
repulsive interactions between the particles can be written as the following
functional integral

$\Xi_{R}=\int\frac{\mathcal{D}\Phi}{\mathcal{N}_{U}}\exp\left[-\frac{\beta}{2}\sum\limits_{\alpha\gamma}(\Phi_{\alpha}U^{-1}_{\alpha\gamma}\Phi_{\gamma})+\sum\limits_{\alpha}\bar{z}_{\alpha}\int
d\mathbb{r}e^{i\beta(q_{\alpha}\varphi+\Phi_{\alpha}+iu_{\alpha})}\right],$
(15)

where $\bar{z}_{\pm}=e^{\frac{\beta{U}_{\pm\pm}(0)}{2}+\frac{\beta
q_{\pm}^{2}G(0)}{2}}z_{\pm}\theta_{\pm}$ are the renormalized fugacities. The
inverse operator, $U^{-1}$, is determined by the integral relation

$\int
d\mathbb{r}^{\prime\prime}\sum\limits_{\lambda}U_{\alpha\lambda}^{-1}(\mathbb{r}-\mathbb{r}^{\prime\prime})U_{\lambda\gamma}(\mathbb{r}^{\prime\prime}-\mathbb{r}^{\prime})=\delta_{\alpha\gamma}\delta(\mathbb{r}-\mathbb{r}^{\prime}).$
(16)

Also, using the cumulant expansion, we obtain

$\Xi_{R}=\exp\left[\sum\limits_{n=1}^{\infty}\frac{\left<W^{n}[\Phi;\chi]\right>_{c}}{n!}\right],$
(17)

where we have introduced the functional

$W[\Phi;\chi]=\sum\limits_{\alpha}\bar{z}_{\alpha}\int
d\mathbb{r}e^{i\beta\Phi_{\alpha}+i\beta\chi_{\alpha}}$ (18)

with the auxiliary variables $\chi_{\alpha}=q_{\alpha}\varphi+iu_{\alpha}$;
and the symbol $\left<(\cdot)\right>_{c}$ denotes the cumulant average Kubo
(1962), so that the following relations take place

$\left<W[\Phi;\chi]\right>_{c}=\left<W[\Phi;\chi]\right>,$ (19)
$\left<W^{2}[\Phi;\chi]\right>_{c}=\left<W^{2}[\Phi;\chi]\right>-\left<W[\Phi;\chi]\right>^{2},$
(20)
$\left<W^{3}[\Phi;\chi]\right>_{c}=\left<W^{3}[\Phi;\chi]\right>-3\left<W^{2}[\Phi;\chi]\right>\left<W[\Phi;\chi]\right>+2\left<W[\Phi;\chi]\right>^{3},...$
(21)

where we have introduced the notation for averaging over the Gaussian measure,
i.e.

$\left<(\cdot)\right>=\int\frac{\mathcal{D}\Phi}{\mathcal{N}_{U}}\exp\left[-\frac{\beta}{2}\sum\limits_{\alpha\gamma}(\Phi_{\alpha}U^{-1}_{\alpha\gamma}\Phi_{\gamma})\right](\cdot).$
(22)

Therefore, using the binomial formula, we obtain

$\sum\limits_{n=1}^{\infty}\frac{\left<W^{n}[\Phi;\chi]\right>_{c}}{n!}=\sum\limits_{n=1}^{\infty}\sum\limits_{k=0}^{n}\frac{\bar{z}_{+}^{k}\bar{z}_{-}^{n-k}}{k!(n-k)!}\int
d\mathbb{r}_{1}^{(+)}...\int d\mathbb{r}_{k}^{(+)}\int
d\mathbb{r}_{1}^{(-)}...\int d\mathbb{r}_{n-k}^{(-)}$
$\times\left<e^{i\beta\sum\limits_{s=1}^{k}\Phi_{+}\left(\mathbb{r}_{s}^{(+)}\right)+i\beta\sum\limits_{l=1}^{n-k}\Phi_{-}\left(\mathbb{r}_{s}^{(-)}\right)}\right>_{c}e^{i\beta\sum\limits_{s=1}^{k}\chi_{+}\left(\mathbb{r}_{s}^{(+)}\right)+i\beta\sum\limits_{l=1}^{n-k}\chi_{-}\left(\mathbb{r}_{l}^{(-)}\right)}.$
(23)

Then, after the calculations of the Gaussian functional integrals Zinn-Justin
(2002), we get

$\left<e^{i\beta\sum\limits_{s=1}^{k}\Phi_{+}(\mathbb{r}_{s}^{(+)})+i\beta\sum\limits_{l=1}^{n-k}\Phi_{-}(\mathbb{r}_{s}^{(-)})}\right>_{c}=e^{\frac{\beta
kU_{++}(0)}{2}+\frac{\beta(n-k)U_{--}(0)}{2}}B^{nk}_{R}\left(\mathbb{r}_{1}^{(+)},...,\mathbb{r}_{k}^{(+)};\mathbb{r}_{1}^{(-)},...,\mathbb{r}_{n-k}^{(-)}\right),$
(24)

where we have introduced the cluster functions:

$B_{R}^{11}\left(\mathbb{r}^{(+)}\right)=B_{R}^{10}\left(\mathbb{r}^{(-)}\right)=1,$
(25)
$B_{R}^{22}\left(\mathbb{r}_{1}^{(+)},\mathbb{r}_{2}^{(+)}\right)=e^{-\beta
U_{++}\left(\mathbb{r}_{1}^{(+)}-\mathbb{r}_{2}^{(+)}\right)}-1,\leavevmode\nobreak\
B_{R}^{21}\left(\mathbb{r}_{1}^{(+)};\mathbb{r}_{1}^{(-)}\right)=e^{-\beta
U_{+-}\left(\mathbb{r}_{1}^{(+)}-\mathbb{r}_{1}^{(-)}\right)}-1,$
$B_{R}^{20}\left(\mathbb{r}_{1}^{(-)},\mathbb{r}_{2}^{(-)}\right)=e^{-\beta
U_{--}\left(\mathbb{r}_{1}^{(-)}-\mathbb{r}_{2}^{(-)}\right)}-1,$ (26)
$B_{R}^{33}\left(\mathbb{r}_{1}^{(+)},\mathbb{r}_{2}^{(+)},\mathbb{r}_{3}^{(+)}\right)=e^{-\beta\left(U_{++}\left(\mathbb{r}_{1}^{(+)}-\mathbb{r}_{2}^{(+)}\right)+U_{++}\left(\mathbb{r}_{1}^{(+)}-\mathbb{r}_{3}^{(+)}\right)+U_{++}\left(\mathbb{r}_{2}^{(+)}-\mathbb{r}_{3}^{(+)}\right)\right)}$
$-e^{-\beta
U_{++}\left(\mathbb{r}_{1}^{(+)}-\mathbb{r}_{2}^{(+)}\right)}-e^{-\beta
U_{++}\left(\mathbb{r}_{1}^{(+)}-\mathbb{r}_{3}^{(+)}\right)}-e^{-\beta
U_{++}\left(\mathbb{r}_{2}^{(+)}-\mathbb{r}_{3}^{(+)}\right)}+2,$
$B_{R}^{32}\left(\mathbb{r}_{1}^{(+)},\mathbb{r}_{2}^{(+)};\mathbb{r}_{1}^{(-)}\right)=e^{-\beta\left(U_{++}\left(\mathbb{r}_{1}^{(+)}-\mathbb{r}_{2}^{(+)}\right)+U_{+-}\left(\mathbb{r}_{1}^{(+)}-\mathbb{r}_{1}^{(-)}\right)+U_{+-}\left(\mathbb{r}_{2}^{(+)}-\mathbb{r}_{1}^{(-)}\right)\right)}$
$-e^{-\beta
U_{++}\left(\mathbb{r}_{1}^{(+)}-\mathbb{r}_{2}^{(+)}\right)}-e^{-\beta
U_{+-}\left(\mathbb{r}_{1}^{(+)}-\mathbb{r}_{1}^{(-)}\right)}-e^{-\beta
U_{+-}\left(\mathbb{r}_{2}^{(+)}-\mathbb{r}_{1}^{(-)}\right)}+2,$
$B_{R}^{31}\left(\mathbb{r}_{1}^{(+)};\mathbb{r}_{1}^{(-)},\mathbb{r}_{2}^{(-)}\right)=e^{-\beta\left(U_{+-}\left(\mathbb{r}_{1}^{(+)}-\mathbb{r}_{1}^{(-)}\right)+U_{+-}\left(\mathbb{r}_{1}^{(+)}-\mathbb{r}_{2}^{(-)}\right)+U_{--}\left(\mathbb{r}_{1}^{(-)}-\mathbb{r}_{2}^{(-)}\right)\right)}$
$-e^{-\beta
U_{+-}\left(\mathbb{r}_{1}^{(+)}-\mathbb{r}_{1}^{(-)}\right)}-e^{-\beta
U_{+-}\left(\mathbb{r}_{1}^{(+)}-\mathbb{r}_{2}^{(-)}\right)}-e^{-\beta
U_{--}\left(\mathbb{r}_{1}^{(-)}-\mathbb{r}_{2}^{(-)}\right)}+2,$
$B_{R}^{30}\left(\mathbb{r}_{1}^{(-)},\mathbb{r}_{2}^{(-)},\mathbb{r}_{3}^{(-)}\right)=e^{-\beta\left(U_{--}\left(\mathbb{r}_{1}^{(-)}-\mathbb{r}_{2}^{(-)}\right)+U_{--}\left(\mathbb{r}_{2}^{(-)}-\mathbb{r}_{3}^{(-)}\right)+U_{--}\left(\mathbb{r}_{1}^{(-)}-\mathbb{r}_{3}^{(-)}\right)\right)}$
$-e^{-\beta
U_{--}\left(\mathbb{r}_{1}^{(-)}-\mathbb{r}_{2}^{(-)}\right)}-e^{-\beta
U_{--}\left(\mathbb{r}_{2}^{(-)}-\mathbb{r}_{3}^{(-)}\right)}-e^{-\beta
U_{--}\left(\mathbb{r}_{1}^{(-)}-\mathbb{r}_{3}^{(-)}\right)}+2,...$ (27)

where index $k$ denotes the number of positively charged ions among the group
of $n$ ions. At sufficiently small radii of the steric interactions, the
cluster functions, $B^{nk}_{R}$, significantly deviate from zero, only when
the positions of the particles are very close to each other. This allows us to
adopt the following approximation for cumulant expansion

$\sum\limits_{n=1}^{\infty}\frac{\left<W^{n}[\Phi;\chi]\right>_{c}}{n!}=\sum\limits_{n=1}^{\infty}\sum\limits_{k=0}^{n}\frac{\tilde{z}_{+}^{k}\tilde{z}_{-}^{n-k}}{k!(n-k)!}\int
d\mathbb{r}_{1}^{(+)}..\int d\mathbb{r}_{k}^{(+)}\int
d\mathbb{r}_{1}^{(-)}..\int d\mathbb{r}_{n-k}^{(-)}$ $\times
B^{nk}_{R}\left(\mathbb{r}_{1}^{(+)},...,\mathbb{r}_{k}^{(+)};\mathbb{r}_{1}^{(-)},...,\mathbb{r}_{n-k}^{(-)}\right)e^{i\beta\sum\limits_{s=1}^{k}\chi_{+}\left(\mathbb{r}_{s}^{(+)}\right)+i\beta\sum\limits_{l=1}^{n-k}\chi_{-}\left(\mathbb{r}_{l}^{(-)}\right)}$
$\approx\int
d\mathbb{r}\sum\limits_{n=1}^{\infty}\sum\limits_{k=0}^{n}\frac{\tilde{z}_{+}^{k}\tilde{z}_{-}^{n-k}}{k!(n-k)!}e^{i\beta
k\chi_{+}\left(\mathbb{r}\right)+i\beta(n-k)\chi_{-}\left(\mathbb{r}\right)}$
$\times\int d\mathbb{r}_{2}^{(+)}..\int d\mathbb{r}_{k}^{(+)}\int
d\mathbb{r}_{1}^{(-)}..\int
d\mathbb{r}_{n-k}^{(-)}B^{nk}_{R}\left(\mathbb{r}_{1}^{(+)},...,\mathbb{r}_{k}^{(+)};\mathbb{r}_{1}^{(-)},...,\mathbb{r}_{n-k}^{(-)}\right)$
$=\int
d\mathbb{r}\sum\limits_{n=1}^{\infty}\sum\limits_{k=0}^{n}b_{nk}{\tilde{z}_{+}^{k}\tilde{z}_{-}^{n-k}}e^{i\beta
k\chi_{+}\left(\mathbb{r}\right)+i\beta(n-k)\chi_{-}\left(\mathbb{r}\right)},$

where $\tilde{z}_{\pm}=z_{\pm}e^{\beta q_{\pm}^{2}G(0)/2}\theta_{\pm}$ and

$b_{nk}=\frac{1}{k!(n-k)!}\int d\mathbb{r}_{2}^{(+)}...\int
d\mathbb{r}_{k}^{(+)}\int d\mathbb{r}_{1}^{(-)}...\int
d\mathbb{r}_{n-k}^{(-)}B^{nk}_{R}\left(\mathbb{r}_{1}^{(+)},...,\mathbb{r}_{k}^{(+)};\mathbb{r}_{1}^{(-)},...,\mathbb{r}_{n-k}^{(-)}\right)$
(28)

are the cluster integrals for a two-component fluid Hill (1986). In other
words, the cluster functions in this ’local’ approximation can be approximated
as

$B^{nk}_{R}\left(\mathbb{r}_{1}^{(+)},...,\mathbb{r}_{k}^{(+)};\mathbb{r}_{1}^{(-)},...,\mathbb{r}_{n-k}^{(-)}\right)=k!(n-k)!b_{nk}\prod\limits_{j=2}^{k}\delta\left(\mathbb{r}_{j}^{(+)}-\mathbb{r}_{1}^{(+)}\right)\prod\limits_{l=1}^{n-k}\delta\left(\mathbb{r}_{l}^{(-)}-\mathbb{r}_{1}^{(+)}\right).$
(29)

Therefore, in the local approximation, the cumulant expansion is simply a
cluster expansion Hill (1986) of the bulk pressure as a function of shifted
chemical potentials,
$\bar{\mu}_{\pm}=\mu_{\pm}+i\chi_{\pm}+{q_{\pm}^{2}G(0)}/{2}+k_{B}T\ln\theta_{\pm}$,
i.e.

$\sum\limits_{n=1}^{\infty}\frac{\left<W^{n}[\Phi;\chi]\right>_{c}}{n!}\approx\beta\int
d\mathbb{r}{P(T,\bar{\mu}_{+},\bar{\mu}_{-})}.$ (30)

Redefining the chemical potentials of the species
$\bar{\mu}_{\pm}\rightarrow\bar{\mu}_{\pm}-{q_{\pm}^{2}G(0)}/{2}-k_{B}T\ln\theta_{\pm}$,
we obtain the following local approximation for the grand partition function
of the reference fluid

$\Xi_{R}=\exp\left[\beta\int
d\mathbb{r}P(T,\bar{\mu}_{+},\bar{\mu}_{-})\right],$ (31)

where $P\left(T,\bar{\mu}_{+},\bar{\mu}_{-}\right)$ is the reference fluid
pressure as a function of the intrinsic chemical potentials
$\bar{\mu}_{\pm}=\mu_{\pm}+iq_{\pm}\varphi-u_{\pm}$ of species and
temperature. We would like to note that this approximation is in accordance
with the thermodynamic relation $\Omega=-PV$ for the bulk systems.

Therefore, the grand partition function of the ionic fluid takes the form

$\Xi=\int\frac{\mathcal{D}\varphi}{C}\exp\left[-\beta\Omega_{ef}[\varphi]\right],$
(32)

where we have introduced the following grand thermodynamic potential

$\Omega_{ef}[\varphi]=\frac{1}{2}\left(\varphi
G^{-1}\varphi\right)-i(\rho_{ext}\varphi)-\int
d\mathbb{r}P\left(T,\bar{\mu}_{+},\bar{\mu}_{-}\right),$ (33)

as a functional of the fluctuating electrostatic potential, $\varphi$. We
would like to note that the functional (33) can be considered as the Landau-
Ginzburg-Wilson effective Hamiltonian for the ionic fluid
($\varphi(\mathbb{r})$ is an order parameter) similar to those obtained within
the similar approaches for simple fluids Ivanchenko and Lisyanskii (1984);
Brilliantov (1998); Brilliantov, Rubí, and Budkov (2020).

Then, using the mean-field approximation, i.e. the functional derivative of
the GTP to zero

$\frac{\delta\Omega_{ef}}{\delta\varphi(\mathbb{r})}=0$ (34)

and introducing the self-consistent field potential (saddle-point)
$\varphi^{(MF)}=i\psi$, we arrive at the following self-consistent field
equation

$\Delta\psi(\mathbb{r})=-\frac{1}{\varepsilon_{b}}\left(\rho_{ext}(\mathbb{r})+q_{+}\bar{c}_{+}(\mathbb{r})+q_{-}\bar{c}_{-}(\mathbb{r})\right)$
(35)

where the local concentrations, $\bar{c}_{\pm}(\mathbb{r})$, are determined by
the following thermodynamic relations

$\bar{c}_{\gamma}(\mathbb{r})=\frac{\delta\Omega_{ef}}{\delta{u}_{\gamma}(\mathbb{r})}=\frac{\partial{P\left(T,\bar{\mu}_{+}(\mathbb{r}),\bar{\mu}_{-}(\mathbb{r})\right)}}{\partial{\bar{\mu}_{\gamma}}}=c_{\gamma}\left(T,\bar{\mu}_{+}(\mathbb{r}),\bar{\mu}_{-}(\mathbb{r})\right),$
(36)

where
$c_{\gamma}\left(T,\bar{\mu}_{+}(\mathbb{r}),\bar{\mu}_{-}(\mathbb{r})\right)$
are the concentrations of the reference fluid species. Thus, within the mean-
field approximation, the ionic concentrations are determined as the molecular
concentrations of the reference fluid and depend on the position,
$\mathbb{r}$, via the intrinsic chemical potentials
$\bar{\mu}_{\pm}(\mathbb{r})=\mu_{\pm}-q_{\pm}\psi(\mathbb{r})-u_{\pm}(\mathbb{r})$.
Solving eq. (35) with the appropriate boundary conditions, we can calculate
the GTP in the mean-field approximation as

$\Omega=\Omega_{ef}[\varphi^{(MF)}]=\int
d\mathbb{r}\left(-\frac{\varepsilon_{b}(\nabla{\psi})^{2}}{2}+\rho_{ext}\psi-P(T,\bar{\mu}_{+},\bar{\mu}_{-})\right).$
(37)

However, when the ionic fluid adjoins with conductors possessing fixed
electric potentials or surface charges and there are no fixed volume charges,
we can use the following functional

$\Omega[\psi]=\int
d\mathbb{r}\left(-\frac{\varepsilon_{b}(\nabla{\psi})^{2}}{2}-P(T,\bar{\mu}_{+},\bar{\mu}_{-})\right),$
(38)

related to the following self-consistent field equation

$\Delta\psi(\mathbb{r})=-\frac{1}{\varepsilon_{b}}\left(q_{+}\bar{c}_{+}(\mathbb{r})+q_{-}\bar{c}_{-}(\mathbb{r})\right).$
(39)

Note that equation (35) can be obtained directly by varying functional (37)
with respect to the potential, $\psi$. In other words, functional (37) reaches
the extreme value at solutions to eq. (35). This idea was used in papers
Budkov  _et al._ (2018); Budkov, Kolesnikov, and Kiselev (2016); Maggs and
Podgornik (2016); Budkov  _et al._ (2020); Blossey, Maggs, and Podgornik
(2017) to obtain the self-consistent field equations by varying the GTP
rewritten as a functional of electrostatic potential Budkov, Kolesnikov, and
Kiselev (2016); Maggs and Podgornik (2016) and potentials of specific
interactions Budkov  _et al._ (2018); Blossey, Maggs, and Podgornik (2017) by
using the local Legendre transformation. In this paper, we justify this
heuristic approach, starting from the grand partition function of ionic fluid.

In order to apply this theory to inhomogeneous electrolyte solutions or room
temperature ionic liquids, it is necessary to specify the dependence of
osmotic pressure on chemical potentials, $P(T,\bar{\mu}_{+},\bar{\mu}_{-})$.
For the case of an ideal gas reference system, for which

$P(T,\bar{\mu}_{+},\bar{\mu}_{-})=\frac{k_{B}T}{\Lambda_{+}^{3}}e^{\beta\bar{\mu}_{+}}+\frac{k_{B}T}{\Lambda_{-}^{3}}e^{\beta\bar{\mu}_{-}},$
(40)

where $\Lambda_{\pm}$ are the thermal de Broglie wavelengths, we obtain the
classical Poisson-Boltzmann equation

$\Delta\psi(\mathbb{r})=-\frac{1}{\varepsilon_{b}}\left(q_{+}c_{+,b}e^{-\beta
q_{+}\psi(\mathbb{r})}+q_{-}c_{-,b}e^{-\beta q_{-}\psi(\mathbb{r})}\right),$
(41)

where we used the relations for bulk chemical potentials
$\mu_{\pm}=k_{B}T\ln\left(\Lambda_{\pm}^{3}c_{\pm,b}\right)$, with $c_{\pm,b}$
being the ionic concentrations in the bulk (where $\psi=0$). Another widely
used reference fluid system is a symmetric lattice gas, for which such
dependence has the form Hill (1986); Maggs and Podgornik (2016)

$P(T,\bar{\mu}_{+},\bar{\mu}_{-})=\frac{k_{B}T}{v}\ln\left(1+e^{\beta\bar{\mu}_{+}}+e^{\beta\bar{\mu}_{-}}\right),$
(42)

where $v$ is the lattice cell volume. In this case, when $q_{\pm}=\pm q$ and
$\mu_{\pm}=k_{B}T\ln{cv/(1-2cv)}$ the self-consistent field equation
transforms into the so-called Poisson-Fermi equation Kornyshev (2007); Bazant
_et al._ (2009); Borukhov, Andelman, and Orland (1997)

$\Delta\psi(\mathbb{r})=\frac{2cq}{\varepsilon_{b}}\frac{\sinh{\beta
q\psi(\mathbb{r})}}{1+2cv\left(\cosh\beta q\psi(\mathbb{r})-1\right)},$ (43)

where $c$ is the bulk ionic concentration. In order to take into account the
ionic size asymmetry, we can use the asymmetric lattice gas model as the
reference system, which was systematically considered in paper Maggs and
Podgornik (2016) in the context of the electric double layer theory. To take
into account the ionic size asymmetry, we can also use the hard sphere
mixture, for which the local pressure and chemical potentials can be obtained
from the following excess free energy density within the Percus-Yevick
approximation:

$\beta
f_{ex}=-n_{0}\ln(1-n_{3})+\frac{n_{1}n_{2}}{1-n_{3}}+\frac{n_{2}^{3}}{24\pi(1-n_{3})^{2}}$
(44)

or using the Carnahan-Starling approximation Roth (2010); Maggs and Podgornik
(2016). The auxiliary variables are: $n_{0}=\bar{c}_{+}+\bar{c}_{-}$,
$n_{1}=1/2\leavevmode\nobreak\ d_{+}\bar{c}_{+}+1/2\leavevmode\nobreak\
d_{-}\bar{c}_{-}$, $n_{2}=\pi\leavevmode\nobreak\
d_{+}^{2}\bar{c}_{+}+\pi\leavevmode\nobreak\ d_{-}^{2}\bar{c}_{-}$,
$n_{3}=\pi/6\leavevmode\nobreak\
d_{+}^{3}\bar{c}_{+}+\pi/6\leavevmode\nobreak\ d_{-}^{3}\bar{c}_{-}$. The
pressure and chemical potentials of species ($\alpha=\pm$) can be derived via
standard thermodynamic relations Roth  _et al._ (2002):

$\beta
P=n_{0}\frac{1+n_{3}+n_{3}^{2}}{(1-n_{3})^{3}}+\frac{1/12\pi\leavevmode\nobreak\
n_{2}^{3}+n_{1}n_{2}(1-n_{3})-3n_{0}n_{3}}{(1-n_{3})^{3}},$ (45)
$\beta\bar{\mu}_{\alpha}=\ln\frac{\bar{c}_{\alpha}\Lambda_{\alpha}^{3}}{1-n_{3}}+\frac{n_{2}d_{\alpha}}{2(1-n_{3})}+\pi
d_{\alpha}^{2}\left(\frac{n_{1}}{1-n_{3}}+\frac{n_{2}^{2}}{8\pi(1-n_{3})^{2}}\right)+\frac{\beta\pi
d_{\alpha}^{3}P}{6},$ (46)

where $d_{\alpha}$ are the ionic effective diameters.

###  II.2 Model II. Polarizable ionic fluid mean-field models

Now let us consider an ionic fluid, whose ions carry charges $q_{\pm}+e_{\pm}$
at the center and additional ’peripheral’ charges $-e_{\pm}$ at a fluctuating
distance from it, so that the total charges of the ions are $q_{\pm}$. The
fluctuating distance, $\mathbb{\xi}$, is described by the distribution
functions $g_{\pm}(\mathbb{\xi})$. This model can be used to describe room
temperature ionic liquids and molten salts with polarizable/polar ions. A
similar ’nonlocal’ description of the bound charge of ions and molecules was
also presented in papers Buyukdagli and Ala-Nissila (2013); Budkov (2019,
2020, 2018). The configuration integrals in the grand partition function (61)
take the form

$Q_{n_{+}n_{-}}=\int d\Gamma_{n_{+}n_{-}}\exp\left[-\beta
U_{n_{+}n_{-}}\right],$ (47)

where

$\int
d\Gamma_{n_{+}n_{-}}(\cdot)=\int..\int\prod\limits_{j=1}^{n_{+}}{d\mathbb{r}_{j}^{(+)}}d\mathbb{\xi}^{(+)}_{j}g_{+}\left(\mathbb{\xi}^{(+)}_{j}\right)\prod\limits_{k=1}^{n_{-}}{d\mathbb{r}_{k}^{(-)}}d\mathbb{\xi}^{(-)}_{k}g_{-}\left(\mathbb{\xi}^{(-)}_{k}\right)(\cdot).$
(48)

The total potential energy is determined by the same expression (4) with the
electrostatic energy expressed via the following local charge density

$\hat{\rho}(\mathbb{r})=\hat{\rho}_{i}(\mathbb{r})+{\rho}_{ext}(\mathbb{r}),$
(49)

where

$\hat{\rho}_{i}(\mathbb{r})=\sum\limits_{\alpha}q_{\alpha}\hat{c}_{\alpha}(\mathbb{r})+\sum\limits_{\alpha=\pm}\sum\limits_{j_{\alpha}=1}^{n_{\alpha}}e_{\alpha}\left(\delta\left(\mathbb{r}-\mathbb{r}^{(\alpha)}_{j_{\alpha}}\right)-\delta\left(\mathbb{r}-\mathbb{r}^{(\alpha)}_{j_{\alpha}}-\mathbb{\xi}_{j_{\alpha}}^{(\alpha)}\right)\right)$
(50)

is the microscopic charge density of the ions taking into account their
polarizability (the second term).

Performing the same transformations as in Subsection A, we arrive at

$\Xi=\int\frac{\mathcal{D}\varphi}{C}e^{-\frac{\beta}{2}\left(\varphi
G^{-1}\varphi\right)+i\beta(\rho_{ext}\varphi)}\Xi_{R},$ (51)

where we have introduced the reference grand partition function

$\Xi_{R}=\int\frac{\mathcal{D}\Phi}{\mathcal{N}_{U}}\exp\left[-\frac{\beta}{2}\sum\limits_{\alpha\gamma}(\Phi_{\alpha}U^{-1}_{\alpha\gamma}\Phi_{\gamma})+W[\Phi;\chi]\right]$
(52)

with the functional $W[\Phi;\chi]$ determined by eq. (18) with

$\chi_{\alpha}(\mathbb{r})=q_{\alpha}\varphi(\mathbb{r})+iu_{\alpha}(\mathbb{r})+k_{B}T\ln\left[\int
d\mathbb{\xi}g_{\alpha}(\mathbb{\xi})e^{i\beta
e_{\alpha}\left(\varphi(\mathbb{r})-\varphi(\mathbb{r}+\mathbb{\xi})\right)}\right].$
(53)

Then, performing the same calculations based on the cumulant expansion and
adopting the same approximations as in Subsection A, we obtain the GTP in the
absence of external volume electric charges

$\Omega[\psi]=\int
d\mathbb{r}\left(-\frac{\varepsilon_{b}(\nabla{\psi})^{2}}{2}-P(T,\bar{\mu}_{+},\bar{\mu}_{-})\right),$
(54)

where

$\bar{\mu}_{\alpha}(\mathbb{r})=\mu_{\alpha}-q_{\alpha}\psi(\mathbb{r})+k_{B}T\ln\left[\int
d\mathbb{\xi}g_{\alpha}(\mathbb{\xi})e^{-\beta
e_{\alpha}\left(\psi(\mathbb{r})-\psi(\mathbb{r}+\mathbb{\xi})\right)}\right].$
(55)

Now we would like to specify the probability distribution functions
$g_{\alpha}(\mathbb{\xi})$. We consider the Gaussian distributions

$g_{\alpha}(\mathbb{\xi})=(2\pi\sigma_{\alpha}^{2})^{-3/2}\int\frac{d\mathbb{n}_{\alpha}}{4\pi}\exp\left[-\frac{(\mathbb{\xi}-\mathbb{s}_{\alpha})^{2}}{2\sigma_{\alpha}^{2}}\right],$
(56)

averaged over orientations of unit vectors
$\mathbb{n}_{\alpha}=\mathbb{s}_{\alpha}/|\mathbb{s}_{\alpha}|$, where
$\mathbb{s}_{\alpha}$ is the vector of the average distance between charges
$\pm e_{\alpha}$. In this case, we obtain

$\int d\mathbb{\xi}g_{\alpha}(\mathbb{\xi})e^{-\beta
e_{\alpha}\left(\psi(\mathbb{r})-\psi(\mathbb{r}+\mathbb{\xi})\right)}=\int\frac{d\mathbb{n}_{\alpha}}{4\pi}\int\frac{d\mathbb{\xi}}{(2\pi\sigma_{\alpha}^{2})^{3/2}}e^{-\frac{(\mathbb{\xi}-\mathbb{s}_{\alpha})^{2}}{2\sigma_{\alpha}^{2}}}e^{-\beta
e_{\alpha}\left(\psi(\mathbb{r})-\psi(\mathbb{r}+\mathbb{\xi})\right)}$
$\approx\int\frac{d\mathbb{n}_{\alpha}}{4\pi}\int\frac{d\mathbb{\xi}}{(2\pi\sigma_{\alpha}^{2})^{3/2}}e^{-\frac{(\mathbb{\xi}-\mathbb{s}_{\alpha})^{2}}{2\sigma_{\alpha}^{2}}}e^{-\beta
e_{\alpha}\mathbb{\xi}\nabla\psi(\mathbb{r})}=\frac{\sinh(\beta
p_{\alpha}|\nabla\psi(\mathbb{r})|)}{\beta
p_{\alpha}|\nabla\psi(\mathbb{r})|}e^{\frac{\gamma_{\alpha}(\nabla\psi(\mathbb{r}))^{2}}{2k_{B}T}},$
(57)

where we have used the approximation
$\psi(\mathbb{r}+\mathbb{\xi})\approx\psi(\mathbb{r})+\mathbb{\xi}\nabla\psi(\mathbb{r})$
and introduced the permanent dipole moments,
$p_{\alpha}=e_{\alpha}|\mathbb{s}_{\alpha}|$, and static polarizabilities,
$\gamma_{\alpha}=\beta\sigma_{\alpha}^{2}e_{\alpha}^{2}$, of the ions.
Therefore, we obtain the following expressions for the intrinsic chemical
potentials of the species

$\bar{\mu}_{\alpha}(\mathbb{r})=\mu_{\alpha}-q_{\alpha}\psi(\mathbb{r})+\frac{\gamma_{\alpha}(\nabla\psi(\mathbb{r}))^{2}}{2}+k_{B}T\ln\left[\frac{\sinh(\beta
p_{\alpha}|\nabla\psi(\mathbb{r})|)}{\beta
p_{\alpha}|\nabla\psi(\mathbb{r})|}\right].$ (58)

Variation of the functional (54) with respect to the electrostatic potential,
$\psi(\mathbb{r})$, with an account for eq. (58) leads to the following self-
consistent field equation

$-\nabla\cdot(\epsilon(\mathbb{r})\nabla\psi(\mathbb{r}))=q_{+}\bar{c}_{+}(\mathbb{r})+q_{-}\bar{c}_{-}(\mathbb{r}),$
(59)

where

$\epsilon(\mathbb{r})=\varepsilon_{b}+\sum\limits_{\alpha=\pm}\left(\gamma_{\alpha}+\frac{p_{\alpha}^{2}}{k_{B}T}\frac{L(\beta
p_{\alpha}\mathcal{E}(\mathbb{r}))}{\beta
p_{\alpha}\mathcal{E}(\mathbb{r})}\right)\bar{c}_{\alpha}(\mathbb{r})$ (60)

is the effective permittivity of the medium,
$\mathcal{E}(\mathbb{r})=|\nabla\psi(\mathbb{r})|$ is the electric field
absolute value, $L(x)=\coth x-1/x$ is the Langevin function. We would like to
mention that eqs. (59-60), in case of the symmetric lattice gas model, were
used in ref. Budkov, Zavarzin, and Kolesnikov (2021) to model the ionic
liquids with a polarizable/polar cation on a charged electrode. A similar
theory was also used to model two-component electrolyte solution with strongly
polar cations on charged electrode Budkov  _et al._ (2020).

###  II.3 Model III. Mean-field ion-dipole mixture models (solvent in salt
models)

Now let us consider the case of an ionic fluid with an admixture of
electrically neutral but polarizable molecules. This model can be used to
describe room temperature ionic fluids and molten salts with an additive of
polar solvents (solvent in salt) or electrolyte solutions with an admixture of
highly polar molecules, such as amino acids. As above, we assume that the ions
carry charges $q_{\pm}+e_{\pm}$ at the center and additional charges
$-e_{\pm}$ at a fluctuating distance from it. The neutral molecules also have
charge $e_{0}$ at the center and charge $-e_{0}$ placed at a fluctuating
distance from the center. The probability distribution functions
$g_{\alpha}(\mathbb{\xi})$ ($\alpha=\pm,0$) are the same averaged Gaussian
functions (56). The grand partition function, thus, can be written as follows

$\Xi=\sum\limits_{n_{+}=0}^{\infty}\sum\limits_{n_{-}=0}^{\infty}\sum\limits_{n_{0}=0}^{\infty}\frac{z_{+}^{n_{+}}z_{-}^{n_{-}}z_{0}^{n_{0}}}{n_{+}!n_{-}!n_{0}!}Q_{n_{+}n_{-}n_{0}}$
(61)

with the fugacities,
$z_{\alpha}=\Lambda_{\alpha}^{-3}e^{\beta\mu_{\alpha}}\theta_{\alpha}$,
($\alpha=\pm,0$) of the ionic species and the configuration integrals

$Q_{n_{+}n_{-}n_{0}}=\int d\Gamma_{n_{+}n_{-}n_{0}}\exp\left[-\beta
U_{n_{+}n_{-}n_{0}}\right],$ (62) $\int
d\Gamma_{n_{+}n_{-}n_{0}}(\cdot)=\int..\int\prod\limits_{j=1}^{n_{+}}{d\mathbb{r}_{j}^{(+)}}d\mathbb{\xi}^{(+)}_{j}g_{+}\left(\mathbb{\xi}^{(+)}_{j}\right)\\\
\times\prod\limits_{k=1}^{n_{-}}{d\mathbb{r}_{k}^{(-)}}d\mathbb{\xi}^{(-)}_{k}g_{-}\left(\mathbb{\xi}^{(-)}_{k}\right)\prod\limits_{l=1}^{n_{0}}{d\mathbb{r}_{l}^{(0)}}d\mathbb{\xi}^{(0)}_{l}g_{0}\left(\mathbb{\xi}^{(0)}_{l}\right)(\cdot),$
(63)

where the total potential energy is determined by the same eq. (4). The total
local charge density, in this case, can be written in the form

$\hat{\rho}(\mathbb{r})=\hat{\rho}_{i}(\mathbb{r})+\hat{\rho}_{d}(\mathbb{r})+{\rho}_{ext}(\mathbb{r}),$
(64)

where $\hat{\rho}_{i}(\mathbb{r})$ is already defined above by eq. (50) and

$\hat{\rho}_{d}(\mathbb{r})=\sum\limits_{j=1}^{N_{0}}e_{0}\left(\delta\left(\mathbb{r}-\mathbb{r}^{(0)}_{j}\right)-\delta\left(\mathbb{r}-\mathbb{r}^{(0)}_{j}-\mathbb{\xi}_{j}^{(0)}\right)\right)$
(65)

is the microscopic density of the bound charge of the molecules. Then,
performing the same transformations as in Subsection A, we obtain

$\Xi=\int\frac{\mathcal{D}\varphi}{C}e^{-\frac{\beta}{2}\left(\varphi
G^{-1}\varphi\right)+i\beta(\rho_{ext}\varphi)}\Xi_{R}[\varphi],$ (66)

where we have introduced the reference grand partition function

$\Xi_{R}=\int\frac{\mathcal{D}\Phi}{\mathcal{N}_{U}}\exp\left[-\frac{\beta}{2}\sum\limits_{\alpha\gamma}(\Phi_{\alpha}U^{-1}_{\alpha\gamma}\Phi_{\gamma})+W[\Phi;\chi]\right]$
(67)

with the functional

$W[\Phi;\chi]=\sum\limits_{\alpha=\pm,0}\bar{z}_{\alpha}\int
d\mathbb{r}e^{i\beta\Phi_{\alpha}+i\beta\chi_{\alpha}},$ (68)

where

$\chi_{\alpha}(\mathbb{r})=q_{\alpha}\varphi(\mathbb{r})+iu_{\alpha}(\mathbb{r})+k_{B}T\ln\left[\int
d\mathbb{\xi}g_{\alpha}(\mathbb{\xi})e^{i\beta
e_{\alpha}\left(\varphi(\mathbb{r})-\varphi(\mathbb{r}+\mathbb{\xi})\right)}\right],\leavevmode\nobreak\
\alpha=\pm$ (69) $\chi_{0}(\mathbb{r})=iu_{0}(\mathbb{r})+k_{B}T\ln\left[\int
d\mathbb{\xi}g_{0}(\mathbb{\xi})e^{i\beta
e_{0}\left(\varphi(\mathbb{r})-\varphi(\mathbb{r}+\mathbb{\xi})\right)}\right].$
(70)

After the same approximation, as in Subsection B, in the absence of external
volume charge, we arrive at the following GTP

$\Omega[\psi]=\int
d\mathbb{r}\left(-\frac{\varepsilon_{b}(\nabla{\psi})^{2}}{2}-P(T,\bar{\mu}_{+},\bar{\mu}_{-},\bar{\mu}_{0})\right)$
(71)

and self-consistent field equation

$-\nabla\cdot(\epsilon(\mathbb{r})\nabla\psi(\mathbb{r}))=q_{+}\bar{c}_{+}(\mathbb{r})+q_{-}\bar{c}_{-}(\mathbb{r}),$
(72)

where

$\epsilon(\mathbb{r})=\varepsilon_{b}+\sum\limits_{\alpha=\pm,0}\left(\gamma_{\alpha}+\frac{p_{\alpha}^{2}}{k_{B}T}\frac{L(\beta
p_{\alpha}\mathcal{E}(\mathbb{r}))}{\beta
p_{\alpha}\mathcal{E}(\mathbb{r})}\right)\bar{c}_{\alpha}(\mathbb{r})$ (73)

is the effective permittivity taking into account polarizability of the ions
and molecules and

$\bar{\mu}_{\alpha}(\mathbb{r})=\mu_{\alpha}-q_{\alpha}\psi(\mathbb{r})+\frac{\gamma_{\alpha}(\nabla\psi(\mathbb{r}))^{2}}{2}+k_{B}T\ln\left[\frac{\sinh(\beta
p_{\alpha}|\nabla\psi(\mathbb{r})|)}{\beta
p_{\alpha}|\nabla\psi(\mathbb{r})|}\right],\leavevmode\nobreak\ \alpha=\pm,$
(74)
$\bar{\mu}_{0}(\mathbb{r})=\mu_{0}+\frac{\gamma_{0}(\nabla\psi(\mathbb{r}))^{2}}{2}+k_{B}T\ln\left[\frac{\sinh(\beta
p_{0}|\nabla\psi(\mathbb{r})|)}{\beta p_{0}|\nabla\psi(\mathbb{r})|}\right]$
(75)

are the intrinsic chemical potentials, where $p_{0}=e_{0}|\mathbb{s}_{0}|$ and
$\gamma_{0}=\beta\sigma_{0}^{2}e_{0}^{2}$ are the permanent dipole moment and
static polarizability of the molecules, respectively. The local concentrations
of the species are

$\bar{c}_{\alpha}(\mathbb{r})=\frac{\partial{P}}{\partial{\bar{\mu}_{\alpha}(\mathbb{r})}}=c_{\alpha}\left(T,\bar{\mu}_{+}(\mathbb{r}),\bar{\mu}_{-}(\mathbb{r}),\bar{\mu}_{0}(\mathbb{r})\right),\leavevmode\nobreak\
\alpha=\pm,0.$ (76)

We would like to note that in case when $\gamma_{\pm}=\gamma_{0}=p_{\pm}=0$
and $p_{0}=p\neq 0$ and

$P(T,\bar{\mu}_{+},\bar{\mu}_{-},\bar{\mu}_{0})=\frac{k_{B}T}{\Lambda_{+}^{3}}e^{\beta\bar{\mu}_{+}}+\frac{k_{B}T}{\Lambda_{-}^{3}}e^{\beta\bar{\mu}_{-}}+\frac{k_{B}T}{\Lambda_{0}^{3}}e^{\beta\bar{\mu}_{0}},$
(77)

we obtain the well-known Poisson-Boltzmann-Langevin (PBL) equation Abrashkin,
Andelman, and Orland (2007); Coalson and Duncan (1996)

$-\nabla\cdot(\epsilon(\mathbb{r})\nabla\psi(\mathbb{r}))=q_{+}c_{+,b}e^{-\beta
q_{+}\psi(\mathbb{r})}+q_{-}c_{-,b}e^{-\beta q_{-}\psi(\mathbb{r})},$ (78)

where

$\epsilon(\mathbb{r})=\varepsilon_{b}+\frac{c_{0,b}p^{2}}{k_{B}T}\mathcal{G}(\beta
p\mathcal{E}(\mathbb{r})),$ (79)

$\mathcal{G}(x)=\sinh xL(x)/x^{2}$, $c_{0,b}$ is the bulk concentration of the
molecules. Note that we took into consideration the relations for bulk
chemical potentials $\mu_{\pm}=k_{B}T\ln(c_{\pm,b}\Lambda_{\pm}^{3})$,
$\mu_{0}=k_{B}T\ln(c_{0,b}\Lambda_{0}^{3})$. To take into account the excluded
volume of the ions and molecules, we can use a three-component hard-sphere
system as the reference system in the Percus-Yevick approximation determined
by eqs. (44-46), where $n_{0}=\bar{c}_{+}+\bar{c}_{-}+\bar{c}_{0}$,
$n_{1}=1/2\leavevmode\nobreak\ d_{+}\bar{c}_{+}+1/2\leavevmode\nobreak\
d_{-}\bar{c}_{-}+1/2\leavevmode\nobreak\ d_{0}\bar{c}_{0}$,
$n_{2}=\pi\leavevmode\nobreak\ d_{+}^{2}\bar{c}_{+}+\pi\leavevmode\nobreak\
d_{-}^{2}\bar{c}_{-}+\pi\leavevmode\nobreak\ d_{0}^{2}\bar{c}_{0}$,
$n_{3}=\pi/6\leavevmode\nobreak\
d_{+}^{3}\bar{c}_{+}+\pi/6\leavevmode\nobreak\
d_{-}^{3}\bar{c}_{-}+\pi/6\leavevmode\nobreak\ d_{0}^{3}\bar{c}_{0}$ ($d_{0}$
is the effective diameter of molecules).

##  III Stress tensor and macroscopic forces in ionic fluids

As it has been shown above, the GTP of an ionic fluid for all considered
models in the mean-field approximation can be written as follows

$\Omega=\int d\mathbb{r}\omega(\psi,\nabla\psi),$ (80)

where we have introduced the GTP density

$\omega=-\frac{\varepsilon_{b}(\nabla{\psi})^{2}}{2}-P.$ (81)

We assume that the dependence of the GTP density on coordinates stems only
from the electrostatic potential $\psi(\mathbb{r})$, so that for simplicity we
consider the case $u_{\pm}(\mathbb{r})=0$. The case $u_{\pm}(\mathbb{r})\neq
0$ will be considered below. Variation of the GTP leads to self-consistent
field equation (39), which can be written as follows

$\frac{\partial}{\partial
x_{i}}\left(\frac{\partial\omega}{\partial\psi_{,i}}\right)-\frac{\partial\omega}{\partial\psi}=0,$
(82)

where $x_{i}$ are the Cartesian coordinates and
$\psi_{,i}=\partial{\psi}/\partial{x}_{i}$. Note that we use the Einstein
rule, implying the summation over repeated indices. Then, with using eq. (82),
we get

$\frac{\partial\omega}{\partial x_{i}}=\frac{\partial}{\partial
x_{k}}\left(\frac{\partial\omega}{\partial\psi_{,k}}\right)\psi_{,i}+\left(\frac{\partial\omega}{\partial\psi_{,k}}\right)\frac{\partial\psi_{,k}}{\partial
x_{i}}=\frac{\partial}{\partial
x_{k}}\left(\psi_{,i}\frac{\partial{\omega}}{\partial{\psi_{,k}}}\right)$ (83)

or

$\frac{\partial\sigma_{ik}}{\partial{x_{k}}}=0,$ (84)

where

$\sigma_{ik}=\delta_{ik}\omega-\psi_{,i}\frac{\partial{\omega}}{\partial{\psi_{,k}}}.$
(85)

Eq. (84) determines the local mechanical equilibrium condition for the ionic
fluid with stress tensor, $\sigma_{ik}$. For model I, eq.(85) yields

$\sigma_{ik}=\varepsilon_{b}\left(\mathcal{E}_{i}\mathcal{E}_{k}-\frac{1}{2}\mathcal{E}^{2}\delta_{ik}\right)-\delta_{ik}P_{0},$
(86)

where the first term is the standard Maxwell stress tensor Landau, Lifshitz,
and Pitaevskii (1984) in continuous dielectric media with the electrostatic
field components $\mathcal{E}_{i}=-\psi_{,i}$, whereas the second term is the
stress tensor of an ionic fluid with the local osmotic pressure
$P_{0}(\mathbb{r})=P(T,\bar{\mu}_{+}(\mathbb{r}),\bar{\mu}_{-}(\mathbb{r}))$,
$\delta_{ik}$ is the unit tensor (Kronecker symbol). For model II, based on
the GTP (54) and intrinsic chemical potentials (58), we obtain the stress
tensor with an account for the polarizabilities and permanent dipole moments
of the ions. Using eq. (85), we obtain the expression for the stress tensor

$\sigma_{ik}=\epsilon\mathcal{E}_{i}\mathcal{E}_{k}-\delta_{ik}\left(P_{0}+\frac{\varepsilon_{b}\mathcal{E}^{2}}{2}\right).$
(87)

Note that in the case of unpolarizable ions ($p_{\pm}=\gamma_{\pm}=0$,
$\epsilon=\varepsilon_{b}$), eq. (87) transforms into eq. (86). For model III,
using eq. (85) and taking into account eqs. (71), (74), and (75), we arrive at
the same expression (87) for the stress tensor with
$P_{0}(\mathbb{r})=P(T,\bar{\mu}_{+}(\mathbb{r}),\bar{\mu}_{-}(\mathbb{r}),\bar{\mu}_{0}(\mathbb{r}))$.

We would like to note that the same approach was used in Landau and Lifshitz
(2013) to express the energy-momentum tensor of an arbitrary physical system
through its Lagrangian (Noether’s theorem Bogolyubov and Shirkov (1973)). Note
that even if we cannot obtain the GTP density from the first principles, we
still can utilise eq. (85) to formulate the stress tensor. However, it would
require ad hoc assumptions regarding the GTP functional dependence on the
electrostatic potential and field.

Now let us consider the case when ions are under the influence of external
potential forces with total potential energies $u_{\pm}(\mathbb{r})$. For
example, these can be van der Waals or specific interactions between the ions
and the charged surfaces (electrodes or membranes, for instance). In this
case, the GTP density depends explicitly on the coordinates via the
aforementioned potential energies. Analogous calculations result in the
following equation

$\frac{\partial\sigma_{ik}}{\partial{x_{k}}}-\sum\limits_{\alpha}\bar{c}_{\alpha}\frac{\partial
u_{\alpha}}{\partial x_{i}}=0,$ (88)

which determines the mechanical equilibrium condition of the ionic fluid in
the presence of external volume forces.

Knowledge of the stress tensor allows us to calculate the macroscopic force
acting on a conductive (electrode) or dielectric (a membrane or a colloid
particle) body immersed in an ionic fluid. It can be calculated as the
following integral over the surface of an immersed body

$F_{i}=\oint\limits_{S}\sigma_{ik}n_{k}dS,$ (89)

where $n_{k}$ are the components of the external normal and $dS$ is the
elementary area. Note that in the absence of volume nonelectrostatic forces
($u_{\pm}(\mathbb{r})=0$), the integration in eq. (89) can be performed over
any closed surface around a macroscopic body because the ionic fluid is in
mechanical equilibrium Landau, Lifshitz, and Pitaevskii (1984). In the
presence of nonelectrostatic volume forces, as it follows from the mechanical
equilibrium, the macroscopic force can be calculated as follows

$F_{i}=\oint\limits_{S}\sigma_{ik}n_{k}dS+\int\limits_{V}f_{i}dV,$ (90)

where the first integral is taken over an arbitrary closed surface $S$ around
the macroscopic body, whereas the second integral – over the volume between
surface $S$ and the surface of the macroscopic body;
$f_{i}=-\sum_{\alpha}\bar{c}_{\alpha}{\partial u_{\alpha}}/{\partial x_{i}}$
are the components of the total nonelectrostatic volume density force. In
order to calculate the macroscopic force by eq. (89) or (90), it is necessary
to solve self-consistent field equation (35) with appropriate boundary
conditions 111The boundary conditions are determined by the nature of immersed
body. The boundary conditions for dielectric or conductive body are standard
Landau, Lifshitz, and Pitaevskii (1984). The boundary conditions for surfaces
with regulated charge are proposed were Podgornik (2018). to obtain the local
osmotic pressure, $P_{0}(\mathbb{r})$, and local ionic concentrations,
$\bar{c}_{\alpha}(\mathbb{r})$, of the ions and the electrostatic field on the
body surface and around it. We would like to note that similar theories based
on ad hoc assumptions were applied in Neu (1999); Trizac and Raimbault (1999);
Filippov  _et al._ (2006) to describe the electrostatic interactions between
charged particles in ionic fluid within the ideal gas and symmetric lattice
gas models. We applied our theory to describe the electrosorption-iduced
deformation of electrified porous carbon with nonconvex pores in electrolyte
solutions Kolesnikov, Mazur, and Budkov (2022).

Now it is instructive to discuss how the general mean-field approximation of
disjoining pressure for an ionic liquid in a slit nanopore can be derived.
Placing the origin of $z$ axis on one charged wall and another one at $z=H$,
we can calculate the disjoining pressure Derjaguin, Churaev, and Muller (1987)
as follows

$\Pi=-\frac{\partial\Omega}{\partial H}-P_{b},$ (91)

where $P_{b}=P(T,\mu_{+},\mu_{-},\mu_{0})$ is the ionic fluid pressure in the
bulk. Then, we have

$\frac{\partial\Omega}{\partial
H}=\int\limits_{0}^{H}dz\frac{\delta\Omega}{\delta\psi(z)}\frac{\partial\psi(z)}{\partial
H}+\omega-\psi^{\prime}\frac{\partial\omega}{\partial\psi^{\prime}}\bigg{|}_{z=H}-\int\limits_{0}^{H}dz\sum\limits_{\alpha}\frac{\partial
P}{\partial u_{\alpha}}\frac{\partial u_{\alpha}}{\partial H}$
$=\int\limits_{0}^{H}dz\sum\limits_{\alpha}\bar{c}_{\alpha}(z)\frac{\partial
u_{\alpha}(z)}{\partial H}+\sigma_{zz}(H),$ (92)

where we took into account that ${\partial P}/{\partial u_{\alpha}}=-{\partial
P}/{\partial\bar{\mu}_{\alpha}}=-\bar{c}_{\alpha}$. Therefore, we eventually
obtain the general expression for disjoining pressure

$\Pi=-\sum\limits_{\alpha}\int\limits_{0}^{H}dz\bar{c}_{\alpha}(z)\frac{\partial
u_{\alpha}(z)}{\partial H}-\sigma_{zz}(H)-P_{b},$ (93)

where
$\sigma_{zz}(H)=\omega-\psi^{\prime}{\partial\omega}/{\partial\psi^{\prime}}|_{z=H}={\varepsilon_{b}\mathcal{E}^{2}(H)}/{2}-P_{0}(H)$
is the normal stress at $z=H$. Now let us assume that external potentials are
created by identical walls, i.e.

$u_{\alpha}(z)=\phi_{\alpha}(z)+\phi_{\alpha}(H-z),$ (94)

where $\phi_{\alpha}$ is the single wall potential, and taking into account
that $\bar{c}_{\alpha}(z)=\bar{c}_{\alpha}(H-z)$, we arrive at

$\Pi=-\sum\limits_{\alpha}\int\limits_{0}^{H}dz\bar{c}_{\alpha}(z)\phi_{\alpha}^{\prime}(z)-\sigma_{zz}(H)-P_{b}.$
(95)

Eq. (95) can be rewritten in a form that is more suitable for applications.
Using the mechanical equilibrium equation

$\frac{d\sigma_{zz}}{dz}-\sum\limits_{\alpha}\bar{c}_{\alpha}(z)u_{\alpha}^{\prime}(z)=0,$
(96)

after the integration from $z=H/2$ to $z=H$ we obtain

$\sigma_{zz}(H)=\sigma_{zz}\left(\frac{H}{2}\right)+\sum\limits_{\alpha}\int\limits_{H/2}^{H}dz\bar{c}_{\alpha}(z)\phi_{\alpha}^{\prime}(z)-\sum\limits_{\alpha}\int\limits_{0}^{H/2}dz\bar{c}_{\alpha}(z)\phi_{\alpha}^{\prime}(z),$
(97)

where we have taken into account eq. (94) and condition
$\bar{c}_{\alpha}(z)=\bar{c}_{\alpha}(H-z)$. Substituting expression (97) for
$\sigma_{zz}(H)$ in eq.(95), after some algebra, we obtain

$\Pi=P_{m}-P_{b}-2\sum\limits_{\alpha}\int\limits_{H/2}^{H}dz\bar{c}_{\alpha}(z)\phi_{\alpha}^{\prime}(z),$
(98)

where $P_{m}=P_{0}\left({H}/{2}\right)$ is the pressure in the middle of a
pore and we take into account that $\mathcal{E}(H/2)=0$. Note that a similar
expression for disjoining pressure for simple (nonionic) fluids was obtained
by Henderson in Ref.Henderson (1986) within a different approach. For
negligible volume forces ($\phi_{\alpha}(z)=0$), we obtain the Langmuir
equation Derjaguin, Churaev, and Muller (1987), $\Pi=P_{m}-P_{b}$, for
disjoining pressure. Eq. (98) is a generalization of the classical DLVO
expression for disjoining pressure for the case of an arbitrary reference
fluid model in the presence of external volume forces.

##  IV Concluding remarks

In conclusion, we proposed a field-theoretical approach based on the
thermodynamic perturbation theory allowing us to derive the grand
thermodynamic potential of the inhomogeneous ionic fluid as a functional of
electrostatic potential for an arbitrary reference system. We obtained the
modified Poisson-Boltzmann equation as the Euler-Lagrange equation for the
grand thermodynamic potential. Using Noether’s theorem, we derived a general
mean-field expression for the stress tensor consistent with the respective
modified Poisson-Boltzmann equation. We derived a general mean-field
expression for the disjoining pressure of an ionic fluid in a slit pore. We
applied the developed formalism to describe three ionic liquid models of
practical importance: nonpolarizable models (including the Poisson-Boltzmann
and Poisson-Fermi equations), polarizable models (the ions have nonzero
permanent dipole moment or static polarizability), and ion-dipole mixtures
(mixture of polar/polarizable ions and polar/polarizable molecules). For all
these models, we derived modified Poisson-Boltzmann equations and respective
stress tensors, which can be used in different applications, where it is
necessary to estimate the electrostatic potential and ionic concentrations
around dielectric or conductive bodies (electrodes, colloids, membranes, dust
particles, etc.) immersed in an ionic fluid, as well as the macroscopic forces
acting on them.

## Acknowledgments

The development of the self-consistent field theories presented in section 2
is supported by the Russian Science Foundation (No. 21-11-00031). Theory
presented in section 3 is supported by grant of the President of the Russian
Federation (project No. MD-341.2021.1.3).

## References

  * Kralj-Iglič and Iglič (1996) V. Kralj-Iglič and A. Iglič, ‘‘A simple statistical mechanical approach to the free energy of the electric double layer including the excluded volume effect,’’ Journal de Physique II 6, 477–491 (1996).
  * Borukhov, Andelman, and Orland (1997) I. Borukhov, D. Andelman, and H. Orland, ‘‘Steric effects in electrolytes: A modified poisson-boltzmann equation,’’ Physical review letters 79, 435 (1997).
  * Maggs and Podgornik (2016) A. Maggs and R. Podgornik, ‘‘General theory of asymmetric steric interactions in electrostatic double layers,’’ Soft matter 12, 1219–1229 (2016).
  * Kornyshev (2007) A. A. Kornyshev, ‘‘Double-layer in ionic liquids: paradigm change?’’ The Journal of Physical Chemistry B 111, 5545–5557 (2007).
  * Goodwin, Feng, and Kornyshev (2017) Z. A. Goodwin, G. Feng, and A. A. Kornyshev, ‘‘Mean-field theory of electrical double layer in ionic liquids with account of short-range correlations,’’ Electrochimica Acta 225, 190–197 (2017).
  * Blossey, Maggs, and Podgornik (2017) R. Blossey, A. Maggs, and R. Podgornik, ‘‘Structural interactions in ionic liquids linked to higher-order poisson-boltzmann equations,’’ Physical Review E 95, 060602 (2017).
  * Frydel (2011) D. Frydel, ‘‘Polarizable poisson–boltzmann equation: The study of polarizability effects on the structure of a double layer,’’ The Journal of chemical physics 134, 234704 (2011).
  * Budkov  _et al._ (2020) Y. A. Budkov, A. V. Sergeev, S. V. Zavarzin, and A. L. Kolesnikov, ‘‘Two-component electrolyte solutions with dipolar cations on a charged electrode: Theory and computer simulations,’’ The Journal of Physical Chemistry C 124, 16308–16314 (2020).
  * Budkov, Zavarzin, and Kolesnikov (2021) Y. A. Budkov, S. V. Zavarzin, and A. L. Kolesnikov, ‘‘Theory of ionic liquids with polarizable ions on a charged electrode,’’ The Journal of Physical Chemistry C 125, 21151–21159 (2021).
  * Budkov, Kolesnikov, and Kiselev (2015) Y. A. Budkov, A. Kolesnikov, and M. Kiselev, ‘‘A modified poisson-boltzmann theory: Effects of co-solvent polarizability,’’ EPL (Europhysics Letters) 111, 28002 (2015).
  * Budkov, Kolesnikov, and Kiselev (2016) Y. A. Budkov, A. Kolesnikov, and M. Kiselev, ‘‘On the theory of electric double layer with explicit account of a polarizable co-solvent,’’ The Journal of Chemical Physics 144, 184703 (2016).
  * Budkov  _et al._ (2018) Y. A. Budkov, A. L. Kolesnikov, Z. A. Goodwin, M. G. Kiselev, and A. A. Kornyshev, ‘‘Theory of electrosorption of water from ionic liquids,’’ Electrochimica Acta 284, 346–354 (2018).
  * Coalson and Duncan (1996) R. D. Coalson and A. Duncan, ‘‘Statistical mechanics of a multipolar gas: a lattice field theory approach,’’ The Journal of Physical Chemistry 100, 2612–2620 (1996).
  * Abrashkin, Andelman, and Orland (2007) A. Abrashkin, D. Andelman, and H. Orland, ‘‘Dipolar poisson-boltzmann equation: ions and dipoles close to charge interfaces,’’ Physical review letters 99, 077801 (2007).
  * Iglič, Gongadze, and Bohinc (2010) A. Iglič, E. Gongadze, and K. Bohinc, ‘‘Excluded volume effect and orientational ordering near charged surface in solution of ions and langevin dipoles,’’ Bioelectrochemistry 79, 223–227 (2010).
  * Buyukdagli and Ala-Nissila (2013) S. Buyukdagli and T. Ala-Nissila, ‘‘Microscopic formulation of nonlocal electrostatics in polar liquids embedding polarizable ions,’’ Physical Review E 87, 063201 (2013).
  * Slavchov (2014) R. I. Slavchov, ‘‘Quadrupole terms in the maxwell equations: Debye-hückel theory in quadrupolarizable solvent and self-salting-out of electrolytes,’’ The Journal of Chemical Physics 140, 074503 (2014).
  * Podgornik (2018) R. Podgornik, ‘‘General theory of charge regulation and surface differential capacitance,’’ The Journal of chemical physics 149, 104701 (2018).
  * Uematsu, Netz, and Bonthuis (2018) Y. Uematsu, R. R. Netz, and D. J. Bonthuis, ‘‘The effects of ion adsorption on the potential of zero charge and the differential capacitance of charged aqueous interfaces,’’ Journal of Physics: Condensed Matter 30, 064002 (2018).
  * Buyukdagli (2020) S. Buyukdagli, ‘‘Schwinger-dyson equations for composite electrolytes governed by mixed electrostatic couplings strengths,’’ The Journal of chemical physics 152, 014902 (2020).
  * de Souza and Bazant (2020) J. P. de Souza and M. Z. Bazant, ‘‘Continuum theory of electrostatic correlations at charged surfaces,’’ The Journal of Physical Chemistry C 124, 11414–11421 (2020).
  * Bazant, Storey, and Kornyshev (2011) M. Z. Bazant, B. D. Storey, and A. A. Kornyshev, ‘‘Double layer in ionic liquids: Overscreening versus crowding,’’ Physical review letters 106, 046102 (2011).
  * Kubo (1962) R. Kubo, ‘‘Generalized cumulant expansion method,’’ Journal of the Physical Society of Japan 17, 1100–1120 (1962).
  * Zinn-Justin (2002) J. Zinn-Justin,  _Quantum field theory and critical phenomena_ , Vol. 113 (Clarendon Press, Oxford, 2002).
  * Hill (1986) T. L. Hill,  _An introduction to statistical thermodynamics_ (Courier Corporation, 1986).
  * Ivanchenko and Lisyanskii (1984) Y. M. Ivanchenko and A. Lisyanskii, ‘‘Ginzburg-landau functional for liquid-vapor phase transition,’’ Theoretical and Mathematical Physics 58, 97–103 (1984).
  * Brilliantov (1998) N. V. Brilliantov, ‘‘Effective magnetic hamiltonian and ginzburg criterion for fluids,’’ Physical Review E 58, 2628 (1998).
  * Brilliantov, Rubí, and Budkov (2020) N. V. Brilliantov, J. M. Rubí, and Y. A. Budkov, ‘‘Molecular fields and statistical field theory of fluids: Application to interface phenomena,’’ Physical Review E 101, 042135 (2020).
  * Bazant  _et al._ (2009) M. Z. Bazant, M. S. Kilic, B. D. Storey, and A. Ajdari, ‘‘Towards an understanding of induced-charge electrokinetics at large applied voltages in concentrated solutions,’’ Advances in colloid and interface science 152, 48–88 (2009).
  * Roth (2010) R. Roth, ‘‘Fundamental measure theory for hard-sphere mixtures: a review,’’ Journal of Physics: Condensed Matter 22, 063102 (2010).
  * Roth  _et al._ (2002) R. Roth, R. Evans, A. Lang, and G. Kahl, ‘‘Fundamental measure theory for hard-sphere mixtures revisited: the white bear version,’’ Journal of Physics: Condensed Matter 14, 12063 (2002).
  * Budkov (2019) Y. A. Budkov, ‘‘A statistical field theory of salt solutions of "hairy" dielectric particles,’’ Journal of Physics: Condensed Matter 32, 055101 (2019).
  * Budkov (2020) Y. A. Budkov, ‘‘Statistical field theory of ion–molecular solutions,’’ Physical Chemistry Chemical Physics 22, 14756–14772 (2020).
  * Budkov (2018) Y. A. Budkov, ‘‘Nonlocal statistical field theory of dipolar particles in electrolyte solutions,’’ Journal of Physics: Condensed Matter 30, 344001 (2018).
  * Landau, Lifshitz, and Pitaevskii (1984) L. D. Landau, E. M. Lifshitz, and L. P. Pitaevskii,  _Course of Theoretical Physics VIII: Electrodynamics of continuous media_ (Pergamon Press, 1984).
  * Landau and Lifshitz (2013) L. D. Landau and E. M. Lifshitz,  _The classical theory of fields_ (Elsevier, 2013).
  * Bogolyubov and Shirkov (1973) N. Bogolyubov and D. Shirkov, ‘‘Introduction to quantum fields theory,’’ Nauka Eds (1973).
  * Note (1) The boundary conditions are determined by the nature of immersed body. The boundary conditions for dielectric or conductive body are standard Landau, Lifshitz, and Pitaevskii (1984). The boundary conditions for surfaces with regulated charge are proposed were Podgornik (2018).
  * Neu (1999) J. C. Neu, ‘‘Wall-mediated forces between like-charged bodies in an electrolyte,’’ Physical review letters 82, 1072 (1999).
  * Trizac and Raimbault (1999) E. Trizac and J.-L. Raimbault, ‘‘Long-range electrostatic interactions between like-charged colloids: Steric and confinement effects,’’ Physical Review E 60, 6530 (1999).
  * Filippov  _et al._ (2006) A. V. Filippov, A. Pal, A. N. Starostin, and A. Ivanov, ‘‘Electrostatic interaction between two macroparticles in the poisson-boltzmann model,’’ JETP letters 83, 546–552 (2006).
  * Kolesnikov, Mazur, and Budkov (2022) A. Kolesnikov, D. Mazur, and Y. Budkov, ‘‘Electrosorption-induced deformation of porous electrode with non-convex pore geometry in electrolyte solutions: theoretical study,’’ to be published (2022).
  * Derjaguin, Churaev, and Muller (1987) B. Derjaguin, N. Churaev, and V. Muller, ‘‘The derjaguin-landau-verwey-overbeek (dlvo) theory of stability of lyophobic colloids,’’ in  _Surface Forces_ (Springer, 1987) pp. 293–310.
  * Henderson (1986) J. Henderson, ‘‘Compressibility route to solvation structure,’’ Molecular Physics 59, 89–96 (1986).