# Statistical field theory of mechanical stresses in Coulomb fluids. General
covariant approach vs Noether’s theorem

Petr E. Brandyshev  Laboratory of Computational Physics, HSE University,
Tallinskaya st. 34, 123458 Moscow, Russia  Yury A. Budkov<EMAIL_ADDRESS>Laboratory of Computational Physics, HSE University, Tallinskaya st. 34,
123458 Moscow, Russia  Laboratory of Multiscale Modeling of Molecular Systems,
G.A. Krestov Institute of Solution Chemistry of the Russian Academy of
Sciences, 153045, Akademicheskaya st. 1, Ivanovo, Russia

###### Abstract

In this paper, we introduce a statistical field theory that describes the
macroscopic forces in inhomogeneous Coulomb fluids. Our approach employs a
general covariant methodology, in conjunction with a functional Legendre
transformation method, to calculate the stress tensor for Coulomb fluids. This
tensor encompasses the mean field stress tensor and the fluctuation
corrections derived through the one-loop approximation. The correction for
fluctuations includes a term that accounts for the thermal fluctuations of the
local electrostatic potential and field in the vicinity of the mean field
configuration. This correlation stress tensor determines how electrostatic
correlation affects local stresses in a nonuniform Coulomb fluid. We also
formulate the generalization of Noether’s first theorem for the case of
fluctuating order parameter and derive within it the same total stress tensor.
We would like to emphasize that our general approach is applicable not only to
Coulomb fluids but also to nonionic simple or complex fluids, for which the
field-theoretic Hamiltonian is known as a functional of the relevant scalar
order parameters.

##  I Introduction

Coulomb fluids, such as plasma, electrolyte solutions, molten salts, and room
temperature ionic liquids, have become a popular topic among chemical
engineers and researchers due to their use in various applications [1]. These
applications include lipid and ion-exchange membranes, biomacromolecules,
colloids, batteries, fuel cells, and supercapacitors, all of which involve
Coulomb fluids interacting with charged surfaces or being confined in charged
nanopores. However, the inhomogeneity of the ionic fluid violates its local
electrical neutrality, which requires the use of numerical solutions of the
self-consistent field equation (Poisson-Boltzmann equation or its
modifications) for the electrostatic potential with appropriate boundary
conditions [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12].

In practical applications of Coulomb fluids confined in nanosized pores of
varying geometries, it is necessary to calculate the mechanical stresses
described by the stress tensor in addition to the concentration and
electrostatic potential profiles. This knowledge of the local stress tensor,
which is consistent with a certain self-consistent field equation, enables the
calculation of physical quantities such as solvation pressure and shear
stresses, which are useful for estimating the deformation of pore materials
(in batteries and supercapacitors, for instance) [13, 14]. Furthermore, the
stress tensor can be used to calculate the macroscopic force acting on the
charged macroscopic conductor or dielectric immersed in the Coulomb fluid [13,
15, 16, 17]. Therefore, a first-principle approach to derive the stress tensor
of inhomogeneous Coulomb fluids would be beneficial for practical application
purposes.

Recent progress has been made in this field, as demonstrated by two recent
articles [3, 18]. In the first [3], Budkov and Kolesnikov applied Noether’s
(first) theorem [19, 20] to the grand thermodynamic potential of a Coulomb
fluid as a functional of the electrostatic potential. They established a
conservation law for the symmetric stress tensor, $\sigma_{ik}$, which
represents the local mechanical equilibrium condition. This tensor is composed
of two terms: the electrostatic Maxwell stress tensor, which is related to the
local electric field, and the hydrostatic isotropic stress tensor, which is
determined by the local osmotic pressure of the ions. The authors extended
this equilibrium condition to include cases where external potential forces
act on the ions. They also derived a general analytical expression for the
electrostatic disjoining pressure of a Coulomb fluid confined in a charged
nanopore slit, which goes beyond the conventional DLVO expression to include
different reference models of fluid. In another study [18], Budkov and Kalikin
presented a self-consistent field theory of macroscopic forces in
inhomogeneous flexible chain polyelectrolyte solutions. The authors derived an
analytical expression for a stress tensor by subjecting the system to a small
dilation and considering the extremum of the grand thermodynamic potential.
This stress tensor includes the previously mentioned hydrostatic and Maxwell
stress tensors, as well as a conformational stress tensor generated by the
conformational Lifshitz entropy of flexible polymer chains. The authors
applied their theory to the investigation of polyelectrolyte solution
constrained in a conducting slit nanopore and observed anomalies in the
disjoining pressure and electric differential capacitance at small pore
thicknesses.

Despite the progress made in deriving the local stress tensor within the mean
field theory, it is still unclear how to compute it beyond the mean field
approximation. This occurs when we consider the fluctuation corrections for
the thermodynamic potential. The fluctuation corrections for the mean field
approximation are always non-local functionals [21, 22], even on a one-loop
correction level, which makes the use of Noether’s theorem, formulated for
local functionals [20], quite problematic. Moreover, it would be valuable to
have a generalization of the Noether’s theorem formalism [20] for the case of
fluctuating order parameters when dealing with the formulation in terms of a
functional integral over the fluctuating order parameters. On the other hand,
an alternative approach to Noether’s theorem for this case may be the general
covariant approach recently proposed by us in paper [23]. Our approach is
based on Noether’s second theorem, which enables us to derive the symmetric
stress tensor for an arbitrary model of an inhomogeneous liquid as a
functional derivative of a grand thermodynamic potential with respect to the
metric tensor. We have applied this approach to several phenomenological
models of inhomogeneous Coulomb fluids – Cahn-Hilliard-like model [24],
Bazant-Storey-Kornyshev model [11], and Maggs-Podgornik-Blossey model [24]. It
should be noted that this method is similar to the one used by D. Hilbert in
the general theory of relativity to derive the energy-momentum tensor from the
action functional [25, 26].

The present paper aims to demonstrate the application of the aforementioned
general covariant approach, combined with the functional Legendre transform,
to obtain fluctuation corrections to the mean field approximation of the
stress tensor for a Coulomb fluid. Additionally, the paper proposes the
application of Noether’s first theorem to the grand partition function, which
is presented in functional integral form.

##  II Functional Legendre transform

We start from the grand partition function as the following functional
integral [3] over the fluctuating electrostatic potential $\phi$ with the
auxiliary volume charge density, $\rho(\mathbb{r})$,

$\Xi[\rho]=\int\frac{\mathcal{D}\phi}{C_{0}}\exp\bigg{(}-\frac{\mathcal{H}[\phi]+(\rho\phi)}{k_{B}T}\bigg{)}$
(1)

where

$\mathcal{H}[\phi]=\int
d\mathbb{r}\bigg{[}\frac{\varepsilon}{2}(\nabla\phi)^{2}-P(\\{\mu_{\alpha}-u_{\alpha}+iq_{\alpha}\phi\\})\bigg{]},$
(2)

is the effective field-theoretic Hamiltonian obtained in [3] and

$C_{0}=\exp\bigg{[}-\frac{1}{2}tr\ln\bigg{(}-\frac{\varepsilon}{k_{B}T}\Delta\bigg{)}\bigg{]},$
(3)

is the normalizing multiplayer of the Gaussian measure;
$P(\\{\bar{\mu}_{\alpha}\\})$ is the pressure of the reference fluid system
(see ref. [3]) dependent on the "shifted" chemical potentials
$\bar{\mu}_{\alpha}=\mu_{\alpha}+iq_{\alpha}\phi-u_{\alpha}$; $k_{B}$ is the
Boltzmann constant, $T$ is the temperature, $u_{\alpha}$ are the external
potentials, $q_{\alpha}$ electric charge of the ion of $\alpha$th kind;
$\Delta=\nabla^{2}$ is the Laplace operator. We have also introduced the
short-hand notation

$(\rho\phi)=\int d\mathbb{r}\rho(\mathbb{r})\phi(\mathbb{r}).$ (4)

We would like to emphasize that the present study only deals with the simplest
model of the Coulomb fluid (Model I, as classified in [3]). In other words,
this model does not account for the static polarizabilities and permanent
dipole moments of the ions and does not explicitly consider the polar solvent.
Nevertheless, a generalization of the theory for the latter cases (Models II
and III) can be performed directly.

The grand partition function can be rewritten as follows

$\Xi[\rho]=\exp\bigg{(}-\frac{W[\rho]}{k_{B}T}\bigg{)},$ (5)

where $W[\rho]$ is some functional of the charge density. Thus, the
expectation value is

$\left<\phi(\mathbb{r})\right>=i\varphi(\mathbb{r})=\frac{\delta
W[\rho]}{\delta\rho(\mathbb{r})}.$ (6)

The effective grand thermodynamic potential (GTP) can be derived from the
following functional Legendre transform [27, 21, 22]

$\Omega[\phi]=W[\rho]-(\rho\phi).$ (7)

Thereby

$\frac{\delta\Omega[\phi]}{\delta\phi(\mathbb{r})}=\int
d\mathbb{r}^{\prime}\frac{\delta
W[\rho]}{\delta\rho(\mathbb{r}^{\prime})}\frac{\delta\rho(\mathbb{r}^{\prime})}{\delta\phi(\mathbb{r})}-\int
d\mathbb{r}^{\prime}\frac{\delta\rho(\mathbb{r}^{\prime})}{\delta\phi(\mathbb{r})}\phi(\mathbb{r}^{\prime})-\rho(\mathbb{r})$
(8)

and with account of eq. (6), we get

$\frac{\delta\Omega[i\varphi]}{\delta\varphi(\mathbb{r})}=-i\rho(\mathbb{r}).$
(9)

The mean field potential (or saddle-point), $\psi(\mathbb{r})$, can be
obtained from the following Euler-Lagrange (EL) equation

$\frac{\delta\mathcal{H}[i\psi]}{\delta\psi(\mathbb{r})}=-i\rho(\mathbb{r}).$
(10)

Introducing the fluctuation, $\eta(\mathbb{r})$, near the mean field potential
by

$\phi=i\psi+\eta,$ (11)

we can expand $\mathcal{H}$ in a functional series in $\eta$

$\frac{\mathcal{H}[i\psi+\eta]}{k_{B}T}=\frac{\mathcal{H}[i\psi]}{k_{B}T}+\frac{\varepsilon}{2k_{B}T}\int
d\mathbb{r}\eta(\mathbb{r})\bigg{(}-\Delta+\varkappa^{2}(\psi)\bigg{)}\eta(\mathbb{r}),$
(12)

where

$\varkappa^{2}(\psi)=\frac{1}{\varepsilon}\sum\limits_{\alpha\gamma}{q_{\alpha}q_{\gamma}}\frac{\partial^{2}P}{\partial\mu_{\alpha}\partial\mu_{\gamma}},$
(13)

and calculating the Gaussian integral (1) over the fluctuations $\eta$, we
arrive at

$W[\rho]=\mathcal{H}[i\psi]+\frac{\Omega_{1}[\psi]}{2\beta}+i(\rho\psi)+O\bigg{(}\frac{1}{\beta^{2}}\bigg{)}.$
(14)

where the functional dependence of $\psi$ on $\rho$ is determined by eq. (10),
$\beta=1/k_{B}T$ and the following notation

$\Omega_{1}[\psi]=tr\ln\bigg{(}\frac{-\Delta+\varkappa^{2}}{-\Delta}\bigg{)}$
(15)

is introduced. We would like to note that $1/\beta$ plays the role of the
Planck’s constant $\hbar$ in quantum field theory [27]. The occurrence of the
multiplayer $1/\beta=k_{B}T$ in the second term on the right-hand side of
equation (14) reflects the fact that this term describes the contribution of
thermal fluctuations of the electrostatic potential near its mean value,
$\varphi$. Note that in the present work we hold only the first-order terms on
$1/\beta$.

Using eqs. (6) and (14), we can get

$i\varphi(\mathbb{r})=\int
d\mathbb{r}^{\prime}\bigg{(}\frac{\delta\mathcal{H}[i\psi]}{\delta\psi(\mathbb{r}^{\prime})}+i\rho(\mathbb{r}^{\prime})\bigg{)}\frac{\delta\psi(\mathbb{r}^{\prime})}{\delta\rho(\mathbb{r})}+i\psi(\mathbb{r})+\frac{i\chi(\mathbb{r})}{2\beta},$
(16)

where

$i\chi(\mathbb{r})=\int
d\mathbb{r}^{\prime}\frac{\delta\Omega_{1}[\psi]}{\delta\psi(\mathbb{r}^{\prime})}\frac{\delta\psi(\mathbb{r}^{\prime})}{\delta\rho(\mathbb{r})}.$
(17)

Taking into account eq. (10) as it follows from eq. (16), we obtain

$\psi(\mathbb{r})=\varphi(\mathbb{r})-\frac{\chi(\mathbb{r})}{2\beta}.$ (18)

Substituting this into eq. (14) and expanding it in a series in $\chi$, we get

$\displaystyle W[\rho]=\mathcal{H}[i\varphi]+i(\rho\varphi)-\int
d\mathbb{r}\bigg{(}\frac{\delta\mathcal{H}[i\psi]}{\delta\psi(\mathbb{r})}+i\rho(\mathbb{r})\bigg{)}\frac{\chi(\mathbb{r})}{2\beta}+\frac{\Omega_{1}[\varphi]}{2\beta}.$
(19)

and using eq. (10) again, we arrive at

$\displaystyle
W[\rho]=\mathcal{H}[i\varphi]+\frac{\Omega_{1}[\varphi]}{2\beta}+i(\rho\varphi).$
(20)

Then performing the Legendre transform (7), we eventually obtain

$\Omega[i\varphi]=\mathcal{H}[i\varphi]+\frac{\Omega_{1}[\varphi]}{2\beta}.$
(21)

Eq. (10) at $\rho=0$ can be written as

$\varepsilon\Delta\psi=-\sum_{\alpha}q_{\alpha}\bar{c}_{\alpha},$ (22)

where we took into account that ${\partial
P(\psi)}/{\partial\psi}=-\sum_{\alpha}q_{\alpha}\bar{c}_{\alpha}$ and
$\partial{P}/\partial{\mu_{\alpha}}=\bar{c}_{\alpha}$. Eq. (22) is nothing
more than the modified Poisson-Boltzmann equation [3].

Variation of $\Omega_{1}$ is derived from the expression

$\frac{\delta\Omega_{1}[\varphi]}{\delta\varphi(\mathbb{r})}=\beta\varepsilon\frac{\partial
U(\varphi)}{\partial\varphi(\mathbb{r})}G(\mathbb{r},\mathbb{r}|\varphi),\quad
U(\varphi)=\varkappa^{2}(\varphi),$ (23)

where the Green function $G(\mathbb{r},\mathbb{r}^{\prime}|\varphi)$ is
determined by equation

$\varepsilon\beta\bigg{(}-\Delta+U(\varphi)\bigg{)}G(\mathbb{r},\mathbb{r}^{\prime}|\varphi)=\delta(\mathbb{r}-\mathbb{r}^{\prime}).$
(24)

Thus the EL equations (9) at $\rho=0$ has the form

$-\frac{\partial
P(\varphi)}{\partial\varphi}+\varepsilon\Delta\varphi+\frac{\varepsilon}{2}\frac{\partial
U(\varphi)}{\partial\varphi(\mathbb{r})}G(\mathbb{r},\mathbb{r}|\varphi)=0.$
(25)

Then using the expansion (16) and holding first-order terms on $1/\beta$, as
mentioned above, we can write the equation. (25) with account of eq. (22) as
follows

$\bigg{(}-\Delta+U(\psi)\bigg{)}\frac{\chi(\mathbb{r})}{\beta}-\frac{\partial
U(\psi)}{\partial\psi(\mathbb{r})}G(\mathbb{r},\mathbb{r}|\psi)=0.$ (26)

Now, using eqs. (10) and (24), we can show that

$\frac{\delta\psi(\mathbb{r}^{\prime})}{\delta\rho(\mathbb{r})}=i\beta
G(\mathbb{r},\mathbb{r}^{\prime}|\psi).$ (27)

Thus, substituting eqs. (23) and (27) in eq. (17), we get

$\chi(\mathbb{r})=\varepsilon\beta^{2}\int
d\mathbb{r}^{\prime}G(\mathbb{r},\mathbb{r}^{\prime}|\psi)G(\mathbb{r}^{\prime},\mathbb{r}^{\prime}|\psi)\frac{\partial
U(\psi)}{\partial\psi(\mathbb{r}^{\prime})}.$ (28)

By examining eq. (24), we can conclude that
$G(\mathbb{r},\mathbb{r}^{\prime}|\psi)$ is a first-order value on $1/\beta$.
As a result, it is expected that $\chi$ would also be a zero-order value on
$1/\beta$, corroborating the expansion given in eq. (16). Furthermore, we can
verify that eq. (28) is indeed a solution to the EL equation (26).

##  III Stress tensor. General covariant approach

Now, we turn to discussion how to calculate the stress tensor from the derived
above GTP (21). To achieve this goal, we will adopt the approach presented in
our recent work [23]. In this approach, the stress tensor can be obtained
using the following expression

$\sigma_{ik}=\frac{2}{\sqrt{g}}\frac{\delta\Omega}{\delta
g_{ik}}\bigg{|}_{g_{ik}=\delta_{ik}},$ (29)

where $g_{ij}$ is the metric tensor, and $g=\det{g_{ij}}$ – its determinant.

To apply eq. (29), we have to express the GTP in general covariant form. Thus,
we have

$\Omega[i\varphi]=\mathcal{H}[i\varphi]+\mathcal{H}^{\prime}[\varphi],\quad\mathcal{H}^{\prime}[\varphi]=\frac{\Omega_{1}[\varphi]}{2\beta},$
(30)

where

$\mathcal{H}=\int
d\mathbb{r}\sqrt{g}\bigg{[}-\frac{\varepsilon}{2}g^{ij}\partial_{i}\varphi\partial_{j}\varphi-P(\\{\bar{\mu}_{\alpha}\\})\bigg{]},$
(31)

is the "mean field" functional with the aforementioned shifted chemical
potentials $\bar{\mu}_{\alpha}=\mu_{\alpha}-q_{\alpha}\varphi$. Note that we
have implied the summation over repeated coordinate indices. We also assumed
here for simplicity $u_{\alpha}=0$.

The one-loop correction is

$\mathcal{H}^{\prime}[\varphi]=\frac{k_{B}T}{2}tr\ln\bigg{(}\frac{-\Delta+U(\varphi)}{-\Delta}\bigg{)},$
(32)

where Laplacian can be written in general covariant form [25, 28]

$\Delta
f=\frac{1}{\sqrt{g}}\partial_{i}\bigg{(}\sqrt{g}g^{ij}\partial_{j}f\bigg{)}.$
(33)

General coordinate transformations lead to

$d\mathbb{r}\rightarrow Jd\mathbb{r},\quad\sqrt{g}\rightarrow J^{-1}\sqrt{g},$
(34)

where $J$ is the Jacobian determinant. Thus, the invariant delta-function can
be defined by

$\int d\mathbb{r}\sqrt{g(\mathbb{r})}\delta(\mathbb{r})=1.$ (35)

Let us consider infinitesimal transformation of the metric tensor

$g_{ij}(\mathbb{r})\rightarrow g_{ij}(\mathbb{r})+\delta g_{ij}(\mathbb{r}).$
(36)

Then, we have [25]

$\delta(\sqrt{g})=\frac{1}{2}\sqrt{g}g^{ij}\delta g_{ij},\quad\delta
g^{ij}=-g^{im}g^{jn}\delta g_{mn}.$ (37)

As it was obtained in ref. [3], the mean field approximation for the stress
tensor is

$\sigma^{(MF)}_{ik}=\varepsilon\partial_{i}\psi\partial_{k}\psi-\frac{\varepsilon}{2}\delta_{ik}\partial_{j}\psi\partial_{j}\psi-P\delta_{ik}.$
(38)

In what follows, utilizing eq. (29), we will calculate the fluctuation
corrections to the mean field approximation (38). The zeroth order term is

$\sigma^{(0)}_{ik}=\frac{2}{\sqrt{g}}\frac{\delta\mathcal{H}[i\varphi]}{\delta
g_{ik}}\bigg{|}_{g_{ik}=\delta_{ik}},$ (39)

that yields the same functional form that realized in the mean field
approximation, i.e.

$\sigma^{(0)}_{ik}=\varepsilon\partial_{i}\varphi\partial_{k}\varphi-\frac{\varepsilon}{2}\delta_{ik}\partial_{j}\varphi\partial_{j}\varphi-P\delta_{ik}.$
(40)

However, it should not be confused with the mean field approximation (38),
since the electrostatic potential $\varphi$ satisfies the equation of EL (25)
taking into account the electrostatic correlations. It can be shown that the
divergence of this tensor is

$\partial_{i}\sigma^{(0)}_{ik}=\bigg{(}-\frac{\partial
P(\varphi)}{\partial\varphi}+\varepsilon\Delta\varphi\bigg{)}\partial_{k}\varphi.$
(41)

Expanding $\sigma^{(0)}_{ik}$ in series on $\chi$ and holding first order
terms in $k_{B}T$, and expanding the chemical potentials
$\mu_{\alpha}=\mu_{\alpha}^{(0)}+\mu_{\alpha}^{(1)}$ up to the first order in
$k_{B}T$, we get

$\sigma^{(0)}_{ik}=\sigma^{(MF)}_{ik}-\delta_{ik}\sum\limits_{\lambda}\mu_{\lambda}^{(1)}\bar{c}_{\lambda}+\frac{\varepsilon}{2\beta}\bigg{(}\partial_{i}\chi\partial_{k}\psi+\partial_{i}\psi\partial_{k}\chi-\delta_{ik}\partial_{j}\chi\partial_{j}\psi\bigg{)}-\delta_{ik}\frac{\chi}{2\beta}\frac{\partial
P}{\partial\psi},$ (42)

where we took into account that
$P(\\{\mu_{\alpha}-q_{\alpha}\psi\\})=P(\\{\mu_{\alpha}^{(0)}-q_{\alpha}\psi\\})+\sum_{\lambda}\mu_{\lambda}^{(1)}\bar{c}_{\lambda}$.
Note that the fluctuation corrections $\mu_{\alpha}^{(1)}$ for the chemical
potentials of the mean field should be calculated for each individual case
[21, 22]. The third and forth terms in the right hand side of eq. (42)
determine a contribution to the total stress tensor arising from the mismatch
between the electrostatic potential of the mean field, $\psi$, and its
expectation value, $\varphi$.

Thus, the stress tensor is

$\sigma_{ik}=\sigma^{(0)}_{ik}+\sigma^{(1)}_{ik},$ (43)

where the first-order term can be obtained in the same way

$\sigma^{(1)}_{ik}=\frac{2}{\sqrt{g}}\frac{\delta\mathcal{H}^{\prime}}{\delta
g_{ik}}\bigg{|}_{g_{ik}=\delta_{ik}}.$ (44)

In order to perform the calculation in eq. (44), we have to introduce a
general covariant definition of the trace, i.e.

$tr(A)=\int d\mathbb{r}\sqrt{g(\mathbb{r})}A(\mathbb{r},\mathbb{r}),$ (45)

For composition of two integral operators,

$C=AB,$ (46)

we can introduce the rule

$C(\mathbb{r},\mathbb{r}^{\prime})=\int
d\mathbb{r}^{\prime\prime}\sqrt{g(\mathbb{r}^{\prime\prime})}A(\mathbb{r},\mathbb{r}^{\prime\prime})B(\mathbb{r}^{\prime\prime},\mathbb{r}^{\prime}).$
(47)

Action of the operator $A$ on a function $f(\mathbb{r})$ is determined by

$Af(\mathbb{r})=\int
d\mathbb{r}^{\prime}\sqrt{g(\mathbb{r}^{\prime})}A(\mathbb{r},\mathbb{r}^{\prime})f(\mathbb{r}^{\prime}).$
(48)

The trace variation is

$\delta tr(A)=\int d\mathbb{r}\sqrt{g(\mathbb{r})}\delta
A(\mathbb{r},\mathbb{r})+\frac{1}{2}\int
d\mathbb{r}\sqrt{g(\mathbb{r})}g^{ij}(\mathbb{r})\delta
g_{ij}(\mathbb{r})A(\mathbb{r},\mathbb{r}),$ (49)

which can be rewritten as

$\delta tr(A)=tr(\bar{\delta}A),$ (50)

where we have introduced the infinitesimal operator $\bar{\delta}A$ that has
the kernel determined by the identity

$\bar{\delta}A(\mathbb{r},\mathbb{r}^{\prime})=\delta
A(\mathbb{r},\mathbb{r}^{\prime})+\frac{1}{2}A(\mathbb{r},\mathbb{r}^{\prime})g^{ij}(\mathbb{r}^{\prime})\delta
g_{ij}(\mathbb{r}^{\prime}).$ (51)

Thus, we have

$\delta tr(A)=\int
d\mathbb{r}\sqrt{g(\mathbb{r})}\bar{\delta}A(\mathbb{r},\mathbb{r}).$ (52)

Using eq. (37), we can show that the operator $\Delta$ is transformed as

$\Delta\rightarrow\Delta+\Delta^{\prime},$ (53)

where

$\displaystyle\Delta^{\prime}f=\frac{1}{2\sqrt{g}}\partial_{i}\bigg{(}\sqrt{g}g^{mn}\delta
g_{mn}g^{ij}\partial_{j}f\bigg{)}-\frac{1}{2\sqrt{g}}g^{mn}\delta
g_{mn}\partial_{i}\bigg{(}\sqrt{g}g^{ij}\partial_{j}f\bigg{)}$ (54)
$\displaystyle-\frac{1}{\sqrt{g}}\partial_{i}\bigg{(}\sqrt{g}\delta
g_{mn}g^{im}g^{jn}\partial_{j}f\bigg{)},$

Functional (32) can be rewritten as

$\mathcal{H}^{\prime}=\frac{k_{B}T}{2}tr\ln\bigg{(}I+\frac{\varepsilon
G_{0}U}{k_{B}T}\bigg{)},$ (55)

where operator $G_{0}$ is determined by equation

$-\frac{\varepsilon}{k_{B}T}\Delta G_{0}=I,$ (56)

where $I$ is the unity operator. Varying both sides of eq. (56) on metric
tensor, we arrive at (see Appendix A)

$-\frac{\varepsilon}{k_{B}T}\bigg{(}\Delta^{\prime}G_{0}+\Delta\bar{\delta}G_{0}\bigg{)}=0.$
(57)

that in turn yields

$\bar{\delta}G_{0}=\frac{\varepsilon}{k_{B}T}G_{0}\Delta^{\prime}G_{0}.$ (58)

Thus, the functional variation of the fluctuation correction is

$\delta\mathcal{H}^{\prime}=\frac{\varepsilon}{2}tr\bigg{[}\bigg{(}I+\frac{\varepsilon
G_{0}U}{k_{B}T}\bigg{)}^{-1}\bar{\delta}G_{0}U\bigg{]}.$ (59)

Therefore, we obtain

$\delta\mathcal{H}^{\prime}=\frac{\varepsilon}{2}tr\bigg{[}\bigg{(}I+\frac{\varepsilon
G_{0}U}{k_{B}T}\bigg{)}^{-1}\frac{\varepsilon}{k_{B}T}G_{0}\Delta^{\prime}G_{0}U\bigg{]}.$
(60)

Let us introduce the operator $G$, determined by the identity (24) that can be
rewritten as follows

$\frac{\varepsilon}{k_{B}T}\bigg{(}-\Delta+U\bigg{)}G=I.$ (61)

Then using (56) and (61), we can get

$\bigg{(}I+\frac{\varepsilon G_{0}U}{k_{B}T}\bigg{)}^{-1}=GG_{0}^{-1}.$ (62)

Thus, the functional variation is

$\delta\mathcal{H}^{\prime}=\frac{\varepsilon^{2}}{2k_{B}T}tr(G\Delta^{\prime}G_{0}U)$
(63)

Using eq. (62) taking into account (56), we arrive at

$\frac{\varepsilon G_{0}U}{k_{B}T}=G_{0}G^{-1}-I.$ (64)

Thus, by using the invariance of the trace with respect to cyclic permutations
of operators in Equation (63), we can obtain the following.

$\displaystyle\delta\mathcal{H}^{\prime}=-\frac{\varepsilon}{2}tr(\Delta^{\prime}\bar{G}).$
(65)

where

$\bar{G}=G-G_{0}.$ (66)

Expression (54) can be rewritten as

$\displaystyle\Delta^{\prime}f=Df-\frac{1}{2}g^{mn}\delta g_{mn}\Delta f,$
(67)

where we have introduced the following auxiliary differential operator

$\displaystyle Df=\frac{1}{2\sqrt{g}}\partial_{i}\bigg{(}\sqrt{g}g^{mn}\delta
g_{mn}g^{ij}\partial_{j}f\bigg{)}-\frac{1}{\sqrt{g}}\partial_{i}\bigg{(}\sqrt{g}\delta
g_{mn}g^{im}g^{jn}\partial_{j}f\bigg{)}.$ (68)

Then, we have

$\displaystyle\delta\mathcal{H}^{\prime}=-\frac{\varepsilon}{2}tr(D\bar{G})+\frac{\varepsilon}{4}tr\bigg{(}g^{mn}\delta
g_{mn}\Delta\bar{G}\bigg{)},$ (69)

where

$\displaystyle tr(DG)=\int d\mathbb{r}\int
d\mathbb{r}^{\prime}\sqrt{g(\mathbb{r}^{\prime})}G(\mathbb{r}^{\prime},\mathbb{r})\sqrt{g(\mathbb{r})}D_{\mathbb{r}}\delta(\mathbb{r}-\mathbb{r}^{\prime})$
(70)

and for brevity we use the notation
$G(\mathbb{r},\mathbb{r}^{\prime})=G(\mathbb{r},\mathbb{r}^{\prime}|\varphi)$.

Integrating it by parts with account of $\delta g_{mn}=0$ on the boundary of
integration, we get

$tr(DG)=-\int d\mathbb{r}^{\prime}\int
d\mathbb{r}\sqrt{g(\mathbb{r}^{\prime})}\partial_{i}G(\mathbb{r}^{\prime},\mathbb{r})\sqrt{g(\mathbb{r})}\delta
g_{mn}(\mathbb{r})\partial_{j}\delta(\mathbb{r}-\mathbb{r}^{\prime})\gamma^{ijmn}(\mathbb{r}),$
(71)

where we have introduced the following tensor function

$\gamma^{ijmn}(\mathbb{r})=\frac{1}{2}g^{mn}(\mathbb{r})g^{ij}(\mathbb{r})-g^{im}(\mathbb{r})g^{jn}(\mathbb{r}).$
(72)

Thus, we have

$\displaystyle\int
d\mathbb{r}^{\prime}\sqrt{g(\mathbb{r}^{\prime})}\partial_{i}G(\mathbb{r}^{\prime},\mathbb{r})\partial_{j}\delta(\mathbb{r}-\mathbb{r}^{\prime})=\lim_{\mathbb{r}^{\prime}\rightarrow\mathbb{r}}\partial_{j}\int
d\mathbb{r}^{\prime\prime}\sqrt{g(\mathbb{r}^{\prime\prime})}\partial_{i}^{\prime}G(\mathbb{r}^{\prime\prime},\mathbb{r}^{\prime})\delta(\mathbb{r}-\mathbb{r}^{\prime\prime})$
(73)

so that

$\displaystyle\int
d\mathbb{r}^{\prime}\sqrt{g(\mathbb{r}^{\prime})}\partial_{i}G(\mathbb{r}^{\prime},\mathbb{r})\partial_{j}\delta(\mathbb{r}-\mathbb{r}^{\prime})=\lim_{\mathbb{r}^{\prime}\rightarrow\mathbb{r}}\partial_{j}\partial_{i}^{\prime}G(\mathbb{r},\mathbb{r}^{\prime}).$
(74)

Therefore, we obtain

$tr(DG)=-\int d\mathbb{r}\sqrt{g(\mathbb{r})}\delta
g_{mn}(\mathbb{r})\gamma^{ijmn}(\mathbb{r})\lim_{\mathbb{r}^{\prime}\rightarrow\mathbb{r}}\partial_{j}\partial_{i}^{\prime}G(\mathbb{r},\mathbb{r}^{\prime}).$
(75)

The second term in eq. (69) is

$tr\bigg{(}g^{mn}\delta g_{mn}\Delta G\bigg{)}=\int
d\mathbb{r}^{\prime}g^{mn}(\mathbb{r}^{\prime})\delta
g_{mn}(\mathbb{r}^{\prime})\sqrt{g(\mathbb{r}^{\prime})}\int
d\mathbb{r}\sqrt{g(\mathbb{r})}G(\mathbb{r},\mathbb{r}^{\prime})\Delta_{\mathbb{r}^{\prime}}\delta(\mathbb{r}^{\prime}-\mathbb{r}).$
(76)

Using the expression

$\int
d\mathbb{r}\sqrt{g(\mathbb{r})}G(\mathbb{r},\mathbb{r}^{\prime})\Delta_{\mathbb{r}^{\prime}}\delta(\mathbb{r}^{\prime}-\mathbb{r})=\lim_{\mathbb{r}^{\prime\prime}\rightarrow\mathbb{r}^{\prime}}\Delta_{\mathbb{r}^{\prime}}\int
d\mathbb{r}\delta(\mathbb{r}^{\prime}-\mathbb{r})\sqrt{g(\mathbb{r})}G(\mathbb{r},\mathbb{r}^{\prime\prime}),$
(77)

we can write

$tr\bigg{(}g^{mn}\delta g_{mn}\Delta G\bigg{)}=\\\ \int
d\mathbb{r}^{\prime}g^{mn}(\mathbb{r}^{\prime})\delta
g_{mn}(\mathbb{r}^{\prime})\sqrt{g(\mathbb{r}^{\prime})}\lim_{\mathbb{r}^{\prime\prime}\rightarrow\mathbb{r}^{\prime}}\bigg{(}\Delta_{\mathbb{r}^{\prime}}G(\mathbb{r}^{\prime},\mathbb{r}^{\prime\prime})\bigg{)}$
(78)

Then, using eq. (44), we can get

$\sigma^{(1)}_{ij}(\mathbb{r})=\frac{\varepsilon}{2}\lim_{\mathbb{r}^{\prime}\rightarrow\mathbb{r}}\bigg{(}\delta_{ij}\partial_{k}\partial_{k}\bar{G}(\mathbb{r},\mathbb{r}^{\prime})+\delta_{ij}\partial_{k}\partial_{k}^{\prime}\bar{G}(\mathbb{r},\mathbb{r}^{\prime})-\partial_{i}\partial_{j}^{\prime}\bar{G}(\mathbb{r},\mathbb{r}^{\prime})-\partial_{j}\partial_{i}^{\prime}\bar{G}(\mathbb{r},\mathbb{r}^{\prime})\bigg{)}.$
(79)

The latter can be expressed in more compact form

$\displaystyle\sigma^{(1)}_{ij}(\mathbb{r})=\frac{\varepsilon}{2}\lim_{\mathbb{r}^{\prime}\rightarrow\mathbb{r}}\hat{D}_{ij}\bar{G}(\mathbb{r},\mathbb{r}^{\prime}),$
(80)

where the following differential operator

$\hat{D}_{ij}=\delta_{ij}\partial_{k}\partial_{k}+\delta_{ij}\partial_{k}\partial_{k}^{\prime}\\\
-\partial_{i}\partial_{j}^{\prime}-\partial_{j}\partial_{i}^{\prime}$ (81)

has been introduced.

The divergence of this tensor is

$\displaystyle\partial_{i}\sigma^{(1)}_{ij}(\mathbb{r})=\frac{\varepsilon}{2}\lim_{\mathbb{r}^{\prime}\rightarrow\mathbb{r}}\bigg{[}\bigg{(}\partial_{i}+\partial_{i}^{\prime}\bigg{)}\bigg{(}\delta_{ij}\partial_{k}\partial_{k}\bar{G}(\mathbb{r},\mathbb{r}^{\prime})+\delta_{ij}\partial_{k}\partial_{k}^{\prime}\bar{G}(\mathbb{r},\mathbb{r}^{\prime})$
(82)
$\displaystyle-\partial_{i}\partial_{j}^{\prime}\bar{G}(\mathbb{r},\mathbb{r}^{\prime})-\partial_{j}\partial_{i}^{\prime}\bar{G}(\mathbb{r},\mathbb{r}^{\prime})\bigg{)}\bigg{]}.$

Therefore, after some algebra, we can get

$\displaystyle\partial_{i}\sigma^{(1)}_{ik}(\mathbb{r})=\frac{\varepsilon}{2}\lim_{\mathbb{r}^{\prime}\rightarrow\mathbb{r}}\bigg{[}\partial_{k}\partial_{i}\partial_{i}\bar{G}(\mathbb{r},\mathbb{r}^{\prime})-\partial_{i}^{\prime}\partial_{i}^{\prime}\partial_{k}\bar{G}(\mathbb{r},\mathbb{r}^{\prime})\bigg{]},$
(83)

and using eqs. (56) and (61), we rewrite eq. (83) in the form

$\displaystyle\partial_{i}\sigma^{(1)}_{ik}(\mathbb{r})=\frac{\varepsilon}{2}G(\mathbb{r},\mathbb{r})\partial_{k}U(\varphi),$
(84)

where

$\partial_{k}U(\varphi)=\frac{\partial
U(\varphi)}{\partial\varphi}\partial_{k}\varphi.$ (85)

Thus, using eqs. (25), (41) and (84), we obtain

$\partial_{i}\sigma_{ik}(\mathbb{r})=\frac{\delta\Omega[\varphi]}{\delta\varphi(\mathbb{r})}\partial_{k}\varphi(\mathbb{r}).$
(86)

where

$\frac{\delta\Omega[\varphi]}{\delta\varphi(\mathbb{r})}=-\frac{\partial
P(\varphi)}{\partial\varphi}+\varepsilon\Delta\varphi+\frac{\varepsilon}{2}\frac{\partial
U(\varphi)}{\partial\varphi(\mathbb{r})}G(\mathbb{r},\mathbb{r}).$ (87)

Thereby, the local mechanical equilibrium condition [3, 19, 29]

$\partial_{i}\sigma_{ik}=0$ (88)

is satisfied if $\varphi$ is a solution of EL equation (25). Thus, we have
demonstrated the exactness of eq. (86) with respect to the terms of any order
in $1/\beta$, for the functional (30). However, as the functional (30) was
only calculated with precision up to the first-order terms, solving equation
(25) with higher-order accuracy would not be meaningful.

Further, we can show that the tensor (79) can be rewritten in the form

$\displaystyle\sigma^{(1)}_{ij}(\mathbb{r})=\sigma^{(cor)}_{ij}(\mathbb{r})=\frac{\varepsilon}{2}\bigg{(}U(\psi)Q(\mathbb{r})+\mathcal{D}_{kk}(\mathbb{r})\bigg{)}\delta_{ij}-\varepsilon\mathcal{D}_{ij}(\mathbb{r}),$
(89)

where we took into account terms of first order on $1/\beta$ and introduce the
following short-hand notations

$Q(\mathbb{r})=\bigg{\langle}\eta^{2}(\mathbb{r})\bigg{\rangle}=\lim_{\mathbb{r}^{\prime}\rightarrow\mathbb{r}}G(\mathbb{r},\mathbb{r}^{\prime}|\psi),$
(90)
$\mathcal{D}_{ij}(\mathbb{r})=\lim_{\mathbb{r}^{\prime}\rightarrow\mathbb{r}}\partial_{i}\partial_{j}^{\prime}\bar{G}(\mathbb{r},\mathbb{r}^{\prime}|\psi).$
(91)

where the expectation values are calculated in the Gaussian approximation
derived from the expansion near the mean field configuration. Eq. (89) being
the main result of this paper determines the fluctuation contribution to the
total stress tensor rising from the thermal fluctuations of the local
electrostatic potential and field near their mean field configuration. In
other words, this tensor describes the effect of electrostatic correlation on
local stresses in the inhomogeneous Coulomb fluid. We call it correlation
stress tensor.

We would like to point out that the formulation of our theory bears formal
similarities to that of the theory of van der Waals forces by E.M. Lifshitz
[30]. Nevertheless, our theory deals with pure classical electrostatic
fluctuations, as opposed to the electromagnetic quantum fluctuations presented
in the Lifshitz theory.

In order to calculate the macroscopic force acting on the dielectric or
conducting body immersed in Coulomb fluid, it is necessary to solve the EL
equation (25) for $\varphi(\mathbb{r})$ with appropriate boundary conditions
and equation for the Green function (24) and then calculate the following
surface integral over the body surface [3, 18]

$F_{i}=\oint\limits_{\mathcal{A}}\sigma_{ik}n_{k}d\mathcal{A},$ (92)

where $n_{k}$ is the external normal and $d\mathcal{A}$ is the elementary
area.

##  IV Stress tensor. Noether’s theorem

In this section, we would like to formulate one more approach to derive the
stress tensor of Coulomb fluids taking into account electrostatic field
thermal fluctuations. This approach is based on a generalization of Noether’s
first theorem for the grand partition function presented in the functional
integral form, i.e.

$\Xi=\int\frac{\mathcal{D}\phi}{C_{0}}\exp\bigg{(}-\frac{\mathcal{H}[\phi])}{k_{B}T}\bigg{)},$
(93)

where the functional $\mathcal{H}$ is determined by eq. (2). Let us perform
the global infinitesimal shift transformation of the coordinates

$x^{\prime}_{k}=x_{k}+h_{k},$ (94)

under which the grand partition function have to be invariant, i.e.

$\delta\ln\Xi=-\frac{1}{k_{B}T}\bigg{\langle}\delta\mathcal{H}[\phi]\bigg{\rangle}=0.$
(95)

Then, after some algebra [19, 23], we obtain

$\partial_{i}\bigg{\langle}\hat{\sigma}_{ik}h_{k}\bigg{\rangle}=\bigg{\langle}\frac{\delta\mathcal{H}[\phi]}{\delta\phi}\partial_{k}\phi
h_{k}\bigg{\rangle},$ (96)

where we have introduced the stress tensor determined on the fluctuating
random order parameter (electrostatic potential)

$\hat{\sigma}_{ik}=\frac{\varepsilon}{2}\delta_{ik}\partial_{j}\phi\partial_{j}\phi-\varepsilon\partial_{i}\phi\partial_{k}\phi-P\delta_{ik}.$
(97)

Since $h_{k}$ is the arbitrary constant infinitesimal vector, we get

$\partial_{i}\bigg{\langle}\hat{\sigma}_{ik}\bigg{\rangle}=\bigg{\langle}\frac{\delta\mathcal{H}[\phi]}{\delta\phi}\partial_{k}\phi\bigg{\rangle},$
(98)

where the expectation value is

$\sigma_{ik}=\bigg{\langle}\hat{\sigma}_{ik}\bigg{\rangle}=\frac{1}{\Xi}\int\frac{\mathcal{D}\phi}{C_{0}}\hat{\sigma}_{ik}\exp\bigg{(}-\frac{\mathcal{H}[\phi])}{k_{B}T}\bigg{)}.$
(99)

Further, taking into account that [31]

$\bigg{\langle}\frac{\delta\mathcal{H}[\phi]}{\delta\phi(\mathbb{r})}\partial_{k}\phi(\mathbb{r})\bigg{\rangle}=\partial_{k}\frac{\delta\phi(\mathbb{r})}{\delta\phi(\mathbb{r})}=\partial_{k}\delta(\mathbb{0})=\lim_{\mathbb{r}^{\prime}\rightarrow\mathbb{r}}\bigg{(}\partial_{k}+\partial_{k}^{\prime}\bigg{)}\delta(\mathbb{r}-\mathbb{r}^{\prime})=0,$
(100)

we arrive at Noether’s first theorem, generalized for the case of fluctuating
order parameter

$\partial_{i}\bigg{\langle}\hat{\sigma}_{ik}\bigg{\rangle}=0.$ (101)

The latter expression represents another important result of this paper.

Now, let us estimate the average stress tensor,
$\sigma_{ik}=\left<\hat{\sigma}_{ik}\right>$. For this purpose, we introduce
the fluctuation, $\eta(\mathbb{r})$, near the mean field potential by

$\phi=i\psi+\eta.$ (102)

Substituting (102) into (97), we can show that the fluctuating stress tensor
has the form with accuracy up to terms of order $k_{B}T$

$\displaystyle\hat{\sigma}_{ik}=\sigma^{(MF)}_{ik}-\delta_{ik}\sum\limits_{\lambda}\mu_{\lambda}^{(1)}\bar{c}_{\lambda}+i\varepsilon\bigg{(}\delta_{ik}\partial_{j}\eta\partial_{j}\psi-\partial_{i}\eta\partial_{k}\psi-\partial_{i}\psi\partial_{k}\eta\bigg{)}+i\eta\delta_{ik}\frac{\partial
P}{\partial\psi}$ (103)
$\displaystyle+\frac{\varepsilon}{2}\varkappa^{2}(\psi)\eta^{2}\delta_{ik}+\frac{\varepsilon}{2}\delta_{ik}\partial_{j}\eta\partial_{j}\eta-\varepsilon\partial_{i}\eta\partial_{k}\eta.$

Fluctuation $\eta$ is the value of the order $(k_{B}T)^{\frac{1}{2}}$. Thus,
to take into account corrections to the stress tensor up to the first order on
$k_{B}T$, we have to restrict ourselves to an expansion of $\mathcal{H}$ up to
third order on $\eta$. Then, using EL equation

$\frac{\delta\mathcal{H}[i\psi]}{\delta\psi}=0,$ (104)

we can expand $\mathcal{H}$ in a functional series in $\eta$ as follows

$\displaystyle\frac{\mathcal{H}[i\psi+\eta]}{k_{B}T}=\frac{\mathcal{H}[i\psi]}{k_{B}T}+\frac{\varepsilon}{2k_{B}T}\int
d\mathbb{r}\eta(\mathbb{r})\bigg{(}-\Delta+\varkappa^{2}(\psi)\bigg{)}\eta(\mathbb{r})$
(105) $\displaystyle-\frac{i\varepsilon}{6k_{B}T}\int
d\mathbb{r}\frac{\partial}{\partial\psi}\bigg{(}\varkappa^{2}(\psi)\bigg{)}\eta^{3}(\mathbb{r})+\frac{O(\eta^{4})}{k_{B}T}.$

Then we have

$\displaystyle
e^{-\beta\mathcal{H}[i\psi+\eta]}=e^{-\beta\mathcal{H}[i\psi]}e^{-\frac{1}{2}(\eta
G^{-1}\eta)}\bigg{(}1+\frac{i\varepsilon\beta}{6}\int
d\mathbb{r}\frac{\partial}{\partial\psi}\bigg{(}\varkappa^{2}(\psi)\bigg{)}\eta^{3}(\mathbb{r})+O\bigg{(}\frac{1}{\beta}\bigg{)}\bigg{)},$
(106)

where we have introduced the following short-hand notation

$(\eta G^{-1}\eta)=\int d\mathbb{r}\int
d\mathbb{r}^{\prime}\eta(\mathbb{r})G^{-1}(\mathbb{r},\mathbb{r}^{\prime}|\psi)\eta(\mathbb{r}^{\prime}),$
(107)

Substituting (103) and (106) into (99) and integrating on $\eta$ we get the
stress tensor with accuracy up to terms of order $1/\beta$

$\displaystyle\sigma_{ik}=\sigma^{(0)}_{ik}+\frac{\varepsilon}{2\beta}\bigg{(}\partial_{i}\chi\partial_{k}\psi+\partial_{i}\psi\partial_{k}\chi-\delta_{ik}\partial_{j}\chi\partial_{j}\psi\bigg{)}-\delta_{ik}\frac{\chi}{2\beta}\frac{\partial
P}{\partial\psi}+\sigma^{(cor)}_{ij},$ (108)

where

$\sigma^{(0)}_{ik}=\sigma^{(MF)}_{ik}-\delta_{ik}\sum\limits_{\lambda}\mu_{\lambda}^{(1)}\bar{c}_{\lambda},$
(109)
$\displaystyle\sigma^{(cor)}_{ij}(\mathbb{r})=\frac{\varepsilon}{2}U(\psi)G(\mathbb{r},\mathbb{r}|\psi)+\frac{\varepsilon}{2}\mathcal{D}_{kk}(\mathbb{r})\delta_{ij}-\varepsilon\mathcal{D}_{ij}(\mathbb{r}),$
(110)

and following short-hand notations

$\mathcal{D}_{ij}(\mathbb{r})=\lim_{\mathbb{r}^{\prime}\rightarrow\mathbb{r}}\partial_{i}\partial_{j}^{\prime}G(\mathbb{r},\mathbb{r}^{\prime}|\psi),$
(111) $\chi(\mathbb{r})=\frac{\varepsilon\beta^{2}}{3}\int
d\mathbb{r}^{\prime}\frac{\partial{U(\psi(\mathbb{r}^{\prime}))}}{\partial\psi(\mathbb{r}^{\prime})}\bigg{\langle}\eta(\mathbb{r})\eta^{3}(\mathbb{r}^{\prime})\bigg{\rangle}_{0}$
(112)

have been introduced.

Using determination of the average over the Gaussian measure

$\left<(\cdot)\right>_{0}=\int\frac{\mathcal{D}\eta}{C_{1}}e^{-\frac{1}{2}(\eta
G^{-1}\eta)}(\cdot),$ (113) $C_{1}=\int\mathcal{D}\eta e^{-\frac{1}{2}(\eta
G^{-1}\eta)}.$ (114)

and well-known Wick’s theorem [32], eq. (112) yields

$\chi(\mathbb{r})=\varepsilon\beta^{2}\int
d\mathbb{r}^{\prime}\frac{\partial{U(\psi(\mathbb{r}^{\prime}))}}{\partial\psi(\mathbb{r}^{\prime})}G(\mathbb{r},\mathbb{r}^{\prime})G(\mathbb{r}^{\prime},\mathbb{r}^{\prime}).$
(115)

Subtracting from the obtaned stress tensor the divergenceless tensor

$\frac{\varepsilon}{2}\mathcal{D}^{(free)}_{kk}(\mathbb{r})\delta_{ij}-\varepsilon\mathcal{D}^{(free)}_{ij}(\mathbb{r}),$
(116)

where
$\mathcal{D}^{(free)}_{ij}(\mathbb{r})=\lim_{\mathbb{r}^{\prime}\rightarrow\mathbb{r}}\partial_{i}\partial_{j}^{\prime}G_{0}(\mathbb{r},\mathbb{r}^{\prime})$,
we arrive at the same expression for the correlation stress tensor (89), which
has been obtained within the general covariant approach.

##  V 1D case

In this section, we would like to specify the correlation stress tensor
obtained above for a practically important one-dimensional case, i.e. when the
Coulomb fluid is confined in a slit-like pore or in close proximity to a flat
electrified surface. Using the cylindrical coordinates, we can write

$G(\mathbb{r},\mathbb{r}^{\prime})=G(\bm{\rho}-\bm{\rho}^{\prime}|z,z^{\prime}),$
(117)

where $\bm{\rho}$ is the two-dimensional vector lying in the plane of the pore
and the $z$-axis is perpendicular to the pore plane. Considering the Fourier
transform

$G(\bm{\rho}-\bm{\rho}^{\prime}|z,z^{\prime})=\int\frac{d^{2}\mathbb{q}}{(2\pi)^{2}}e^{-i\mathbb{q}(\bm{\rho}-\bm{\rho}^{\prime})}G(q|z,z^{\prime}),$
(118)

where $G(q|z,z^{\prime})$ is the even function of $\mathbb{q}$ depending only
on the vector modulus $q=|\mathbb{q}|$. This system is one-dimensional in the
sense that the functions $\psi$ and $\mathcal{D}_{ij}$ depend only on $z$. Let
us first calculate the cross elements

$\mathcal{D}_{xy}(z)=\lim_{z^{\prime}\rightarrow
z}\int\frac{d^{2}\mathbb{q}}{(2\pi)^{2}}q_{x}q_{y}G(q|z,z^{\prime}),$ (119)
$\mathcal{D}_{xz}(z)=\lim_{z^{\prime}\rightarrow
z}\int\frac{d^{2}\mathbb{q}}{(2\pi)^{2}}iq_{x}\partial_{z}G(q|z,z^{\prime}),$
(120) $\mathcal{D}_{yz}(z)=\lim_{z^{\prime}\rightarrow
z}\int\frac{d^{2}\mathbb{q}}{(2\pi)^{2}}iq_{y}\partial_{z}G(q|z,z^{\prime}),$
(121)

It is obvious that these elements are zero, since the integral of an odd
function is zero

$\mathcal{D}_{xy}(z)=\mathcal{D}_{xz}(z)=\mathcal{D}_{yz}(z)=0.$ (122)

Diagonal elements are

$\mathcal{D}_{xx}(z)=\int\frac{d^{2}\mathbb{q}}{(2\pi)^{2}}q_{x}^{2}Q(q,z),\leavevmode\nobreak\
\mathcal{D}_{yy}(z)=\int\frac{d^{2}\mathbb{q}}{(2\pi)^{2}}q_{y}^{2}Q(q,z),$
(123)
$\mathcal{D}_{zz}(z)=\int\frac{d^{2}\mathbb{q}}{(2\pi)^{2}}\mathcal{D}_{zz}(q,z),$
(124)

where

$\mathcal{D}_{zz}(q,z)=\lim_{z^{\prime}\rightarrow
z}\partial_{z}\partial_{z^{\prime}}G(q|z,z^{\prime}),$ (125)
$Q(q,z)=\lim_{z^{\prime}\rightarrow z}G(q|z,z^{\prime}).$ (126)

Thus the trace is

$\mathcal{D}_{kk}(z)=\int\frac{d^{2}\mathbb{q}}{(2\pi)^{2}}\bigg{(}q^{2}Q(q,z)+\mathcal{D}_{zz}(q,z)\bigg{)}.$
(127)

The Fourier image of the Green function can be found from the equation

$\varepsilon\beta\bigg{(}-\partial^{2}_{z}+q^{2}+\varkappa^{2}(z)\bigg{)}G(q|z,z^{\prime})=\delta(z-z^{\prime}).$
(128)

The correlation stress tensor elements are defined by the following
expressions

$\displaystyle\sigma^{(cor)}_{ij}(z)=\int\frac{d^{2}\mathbb{q}}{(2\pi)^{2}}\sigma^{(cor)}_{ij}(q,z),$
(129)
$\displaystyle\sigma^{(cor)}_{xx}(q,z)=\frac{\varepsilon}{2}\bigg{(}Q(q,z)\varkappa^{2}(z)+\mathcal{D}_{zz}(q,z)\bigg{)},$
(130) $\displaystyle\sigma^{(cor)}_{yy}(q,z)=\sigma^{(cor)}_{xx}(q,z),$ (131)
$\displaystyle\sigma^{(cor)}_{zz}(q,z)=\frac{\varepsilon}{2}\bigg{(}Q(q,z)\left(q^{2}+\varkappa^{2}(z)\right)-\mathcal{D}_{zz}(q,z)\bigg{)}.$
(132)

By analyzing the asymptotic behavior of the normal stresses for larger pore
thickness, it is possible to better understand the developed theoretical
background. Specifically, the analysis would focus on the behavior of the
normal stresses as the pore thickness $H\to\infty$. In this case,
$\varkappa(z)=\kappa=$const (inverse screening length of bulk Coulomb fluid
[33]) and

$G(q|z,z^{\prime})\simeq\frac{\exp\left(-\kappa_{q}|z-z^{\prime}|\right)}{2\beta\varepsilon\kappa_{q}},$
(133)

so that we have

$Q(q,z)\simeq\frac{1}{2\beta\varepsilon\kappa_{q}},\leavevmode\nobreak\
\mathcal{D}_{zz}(q,z)\simeq-\frac{\kappa_{q}}{2\beta\varepsilon},$ (134)

where $\kappa_{q}=\sqrt{\kappa^{2}+q^{2}}$. Thus, we have for the normal
correlation stress

$\sigma_{zz}^{(cor)}=\frac{k_{B}T}{2}\int\frac{d^{2}\mathbb{q}}{(2\pi)^{2}}\left(\kappa_{q}-q-\frac{\kappa^{2}}{2q}\right)=-\frac{k_{B}T\kappa^{3}}{12\pi},$
(135)

where we have subtracted from the final result the infinite value. To obtain
the total normal stress at $H\to\infty$, we have to calculate (see eq. (42))

$\sigma_{zz}^{(0)}=-k_{B}T\sum\limits_{\alpha}c_{\alpha}-\sum\limits_{\alpha}c_{\alpha}\mu_{\alpha}^{(1)},$
(136)

where $c_{\alpha}$ are the bulk ionic concentrations. This contribution could
be easily obtained for the case, when the reference system is the mixture of
the ideal gases, for which

$P(\\{\bar{\mu}_{\alpha}\\})=k_{B}T\sum\limits_{\alpha}\Lambda_{\alpha}^{-3}e^{\beta\bar{\mu}_{\alpha}},$
(137)

where $\Lambda_{\alpha}$ are the de Broglie thermal wavelengths. In this case,
we have

$\displaystyle\bar{c}_{\alpha}(\mathbb{r})=\frac{\delta\Omega}{\delta
u_{\alpha}(\mathbb{r})}=\Lambda_{\alpha}^{-3}e^{\beta\bar{\mu}_{\alpha}}+\frac{\varepsilon}{2}\frac{\partial
U(\varphi)}{\partial u_{\alpha}(\mathbb{r})}G(\mathbb{r},\mathbb{r}|\varphi)$
$\displaystyle=\Lambda_{\alpha}^{-3}e^{\beta\bar{\mu}_{\alpha}}\left(1-\frac{q_{\alpha}^{2}}{2(k_{B}T)^{2}}G(\mathbb{r},\mathbb{r}|\varphi)\right)$
(138)

that yields

$\bar{\mu}_{\alpha}=k_{B}T\ln\left(\Lambda_{\alpha}^{3}\bar{c}_{\alpha}\right)+\frac{q_{\alpha}^{2}}{2k_{B}T}G(\mathbb{r},\mathbb{r}|\varphi).$
(139)

Thus, for the bulk, where $u_{\alpha}=0$ and $\varphi=0$, we have

$\mu_{\alpha}=k_{B}T\ln\left(\Lambda_{\alpha}^{3}{c}_{\alpha}\right)+\frac{q_{\alpha}^{2}}{2k_{B}T}\left(G(0)-G_{0}(0)\right),$
(140)

We subtracted the infinite constant $q_{\alpha}^{2}G_{0}(0)/{2k_{B}T}$ from
the final result to avoid infinite terms. This can be justified by recognizing
that the chemical potentials are determined up to an arbitrary constant. We
also considered that in the bulk the Green function is translation invariant,
i.e. $G(\mathbb{r},\mathbb{r}^{\prime}|0)=G(\mathbb{r}-\mathbb{r}^{\prime})$.
The second term in eq. (140) is the fluctuation correction to the ideal gas
bulk chemical potential which can be rewritten as follows

$\mu_{\alpha}^{(1)}=\frac{q_{\alpha}^{2}}{2k_{B}T}\int\frac{d^{2}\mathbb{q}}{(2\pi)^{2}}\left(G(q|z,z)-G_{0}(q|z,z)\right),$
(141)

where
$G_{0}(q|z,z^{\prime})={\exp\left(-q|z-z^{\prime}|\right)}/{2q\beta\varepsilon}$.
Further, calculating the integral (141) and taking into account that for ideal
gas reference system
$\kappa=\left(\sum_{\alpha}q_{\alpha}^{2}c_{\alpha}/\varepsilon
k_{B}T\right)^{1/2}$, after some algebra, we arrive at the classical Debye-
Hückel expression

$\sigma_{zz}=-P_{o}=-k_{B}T\sum\limits_{\alpha}c_{\alpha}+\frac{k_{B}T\kappa^{3}}{24\pi},$
(142)

where $P_{o}$ is the osmotic pressure of the ions in the bulk. Performing the
same calculations for other nonzero components of the correlation stress
tensor, we obtain

$\sigma^{(cor)}_{xx}=\sigma^{(cor)}_{yy}=\frac{k_{B}T}{4}\int\frac{d^{2}\mathbb{q}}{(2\pi)^{2}}\left(q-\frac{q^{2}}{\sqrt{q^{2}+\kappa^{2}}}-\frac{\kappa^{2}}{2q}\right)=-\frac{k_{B}T\kappa^{3}}{12\pi}.$
(143)

Taking into account that
$\sigma_{xx}^{(0)}=\sigma_{yy}^{(0)}=\sigma_{zz}^{(0)}$, we obtain, as should
be, the isotropic stress tensor in the bulk fluid, that is,
$\sigma_{ij}=-P_{o}\delta_{ij}$.

##  VI Concluding remarks

In conclusion, we have formulated a first-principles statistical field theory
of the macroscopic forces in spatially inhomogeneous Coulomb fluids. Utilizing
the previously formulated general covariant approach, in combination with the
method based on the functional Legendre transform, we derived the total stress
tensor for the Coulomb fluid consisting of the previously derived mean field
stress tensor and fluctuation corrections on the one-loop correction level.
The fluctuation corrections include a term that accounts for the thermal
fluctuations of the local electrostatic potential and field in the vicinity of
the mean field configuration. This correlation stress tensor determines how
electrostatic correlation affects the local stresses in a non-uniform Coulomb
fluid. In addition, we have formulated a generalization of Noether’s first
theorem for the case of a fluctuating order parameter and derived the same
total stress tensor using this formulation.

Finally, we would like to note that proposed approaches is applicable not only
to the Coulomb fluids, but also to nonionic simple or complex fluids
(including quantum ones) for which we know the field-theoretic Hamiltonian as
a functional of the appropriate fluctuating order parameters.

Acknowledgements. This work is an output of a research project implemented as
part of the Basic Research Program at the National Research University Higher
School of Economics (HSE University).

## Appendix A

Let us introduce the function $\tilde{\delta}(\mathbb{r}-\mathbb{r}^{\prime})$
determined by expression

$\int
d\mathbb{r}\sqrt{\tilde{g}(\mathbb{r})}\tilde{\delta}(\mathbb{r}-\mathbb{r}^{\prime})f(\mathbb{r})=f(\mathbb{r}^{\prime}).$
(144)

where we consider the metric variation as the type function variation (it is
not variation under coordinate transformations)

$\tilde{g}_{ij}(\mathbb{r})=g_{ij}(\mathbb{r})+\delta g_{ij}(\mathbb{r}),$
(145)

and then we can define the variation

$\delta^{\prime}(\mathbb{r}-\mathbb{r}^{\prime})=\tilde{\delta}(\mathbb{r}-\mathbb{r}^{\prime})-\delta(\mathbb{r}-\mathbb{r}^{\prime}).$
(146)

Then, using eq. (35), we can show that the equality (144) is fulfilled if the
variation of delta-function is determined by the following identity

$\delta^{\prime}(\mathbb{r}-\mathbb{r}^{\prime})=-\frac{1}{2}g^{ij}(\mathbb{r}^{\prime})\delta
g_{ij}(\mathbb{r}^{\prime})\delta(\mathbb{r}-\mathbb{r}^{\prime}).$ (147)

Let us note that using known delta-function properties we can write

$g^{ij}(\mathbb{r}^{\prime})\delta
g_{ij}(\mathbb{r}^{\prime})\delta(\mathbb{r}-\mathbb{r}^{\prime})=g^{ij}(\mathbb{r})\delta
g_{ij}(\mathbb{r})\delta(\mathbb{r}-\mathbb{r}^{\prime}).$ (148)

It is obviously that in this case

$\delta\bigg{(}\sqrt{g(\mathbb{r}^{\prime})}\delta(\mathbb{r}-\mathbb{r}^{\prime})\bigg{)}=0.$
(149)

Under transformations (145) the operator $G_{0}$ is infinitesimally
transformed as follows

$\tilde{G_{0}}=G_{0}+\delta G_{0}.$ (150)

where the transformed operator $\tilde{G_{0}}$ is determined by the equation

$-\frac{\varepsilon}{k_{B}T}\bigg{(}\Delta_{\mathbb{r}}+\Delta^{\prime}_{\mathbb{r}}\bigg{)}\tilde{G_{0}}(\mathbb{r},\mathbb{r}^{\prime})=\tilde{\delta}(\mathbb{r}-\mathbb{r}^{\prime}).$
(151)

The equation (102) can be written in the form

$-\frac{\varepsilon}{k_{B}T}\Delta_{\mathbb{r}}G_{0}(\mathbb{r},\mathbb{r}^{\prime})=\delta(\mathbb{r}-\mathbb{r}^{\prime}).$
(152)

Then subtracting (152) from expression (151) with account of eq. (147) and
holding only first-order terms, we can get

$-\frac{\varepsilon}{k_{B}T}\bigg{(}\Delta^{\prime}_{\mathbb{r}}G_{0}(\mathbb{r},\mathbb{r}^{\prime})+\Delta_{\mathbb{r}}\delta
G_{0}(\mathbb{r},\mathbb{r}^{\prime})\bigg{)}=-\frac{1}{2}g^{ij}(\mathbb{r}^{\prime})\delta
g_{ij}(\mathbb{r}^{\prime})\delta(\mathbb{r}-\mathbb{r}^{\prime}).$ (153)

Using eq. (152) again, we can rewrite this equality in the form

$-\frac{\varepsilon}{k_{B}T}\bigg{(}\Delta^{\prime}_{\mathbb{r}}G_{0}(\mathbb{r},\mathbb{r}^{\prime})+\Delta_{\mathbb{r}}\delta
G_{0}(\mathbb{r},\mathbb{r}^{\prime})+\frac{1}{2}g^{ij}(\mathbb{r}^{\prime})\delta
g_{ij}(\mathbb{r}^{\prime})\Delta_{\mathbb{r}}G_{0}(\mathbb{r},\mathbb{r}^{\prime})\bigg{)}=0.$
(154)

Using definition (51), eq. (154) can be reduced to eq. (57).

Let us consider the following expressed diffeomorphic transformations

$\displaystyle\tilde{x}^{k}=x^{k}+\xi^{k}(\mathbb{r}),$ (155)
$\displaystyle\tilde{g}_{ij}(\tilde{\mathbb{r}})=g_{ij}(\mathbb{r})+\tilde{\delta}g_{ij}(\mathbb{r}),$

where $\xi^{k}$ and $\tilde{\delta}g_{ij}$ are generally independent
transformation parameters, $\tilde{x}^{k}$ and $x^{k}$ are $k$th components of
vectors $\tilde{\mathbb{r}}$ and $\mathbb{r}$, respectively. Then, using the
identity

$\int
d\tilde{\mathbb{r}}_{1}\sqrt{\tilde{g}(\tilde{\mathbb{r}}_{1})}\tilde{\delta}(\tilde{\mathbb{r}}_{1}-\tilde{\mathbb{r}}_{2})f(\tilde{\mathbb{r}}_{1})=f(\tilde{\mathbb{r}}_{2}).$
(156)

we can show that the delta-function is transformed under (155), in accordance
with the expression

$\tilde{\delta}(\tilde{\mathbb{r}}_{1}-\tilde{\mathbb{r}}_{2})=\bigg{(}1-\partial_{k}\xi^{k}(\mathbb{r}_{1})-\frac{1}{2}g^{ij}(\mathbb{r}_{1})\tilde{\delta}g_{ij}(\mathbb{r}_{1})\bigg{)}\delta(\mathbb{r}_{1}-\mathbb{r}_{2}).$
(157)

and if we put that

$\tilde{\delta}g^{ik}(\mathbb{r})=\xi^{i,k}(\mathbb{r})+\xi^{k,i}(\mathbb{r}),$
(158)

where

$\xi^{i,k}(\mathbb{r})=g^{kj}\partial_{j}\xi^{i}(\mathbb{r}),$ (159)

then in this special case

$\tilde{\delta}(\tilde{\mathbb{r}}_{1}-\tilde{\mathbb{r}}_{2})=\delta(\mathbb{r}_{1}-\mathbb{r}_{2}).$
(160)

The delta function thus defined is invariant under a subgroup of
transformations (155) that correspond to eq. (158) (under standard
diffeomorphisms). Our approach requires us to consider another subgroup of
(155) that corresponds to the condition that $\xi^{k}=0$, and where
$\tilde{\delta}g_{ij}$ can be an arbitrary infinitesimal function. The delta
function is transformed under this subgroup according to the law (147).

## References

  * Naji  _et al._ [2013] A. Naji, M. Kanduč, J. Forsman, and R. Podgornik, The Journal of chemical physics 139, 150901 (2013).
  * Blossey [2023] R. Blossey, in  _The Poisson-Boltzmann Equation: An Introduction_ (Springer, 2023) pp. 53–96.
  * Budkov and Kolesnikov [2022] Y. A. Budkov and A. L. Kolesnikov, Journal of Statistical Mechanics: Theory and Experiment 2022, 053205 (2022).
  * Budkov and Kolesnikov [2021] Y. A. Budkov and A. L. Kolesnikov, Current Opinion in Electrochemistry 33, 100931 (2021).
  * Borukhov, Andelman, and Orland [1997] I. Borukhov, D. Andelman, and H. Orland, Physical review letters 79, 435 (1997).
  * Abrashkin, Andelman, and Orland [2007] A. Abrashkin, D. Andelman, and H. Orland, Physical review letters 99, 077801 (2007).
  * Avni, Adar, and Andelman [2020] Y. Avni, R. M. Adar, and D. Andelman, Physical Review E 101, 010601 (2020).
  * Podgornik [2018] R. Podgornik, The Journal of chemical physics 149, 104701 (2018).
  * Kornyshev [2007] A. A. Kornyshev, The Journal of Physical Chemistry B 111, 5545 (2007).
  * Iglič, Gongadze, and Bohinc [2010] A. Iglič, E. Gongadze, and K. Bohinc, Bioelectrochemistry 79, 223 (2010).
  * Bazant, Storey, and Kornyshev [2011] M. Z. Bazant, B. D. Storey, and A. A. Kornyshev, Physical review letters 106, 046102 (2011).
  * Slavchov [2014] R. I. Slavchov, The Journal of Chemical Physics 140, 074503 (2014).
  * Kolesnikov, Mazur, and Budkov [2022] A. L. Kolesnikov, D. A. Mazur, and Y. A. Budkov, Europhysics Letters 140, 16001 (2022).
  * Kolesnikov, Budkov, and Gor [2021] A. Kolesnikov, Y. A. Budkov, and G. Gor, Journal of Physics: Condensed Matter 34, 063002 (2021).
  * Neu [1999] J. C. Neu, Physical review letters 82, 1072 (1999).
  * Trizac and Raimbault [1999] E. Trizac and J.-L. Raimbault, Physical Review E 60, 6530 (1999).
  * de Souza and Bazant [2020] J. P. de Souza and M. Z. Bazant, The Journal of Physical Chemistry C 124, 11414 (2020).
  * Budkov and Kalikin [2023] Y. A. Budkov and N. N. Kalikin, Physical Review E 107, 024503 (2023).
  * Hermann and Schmidt [2022] S. Hermann and M. Schmidt, Journal of Physics: Condensed Matter 34, 213001 (2022).
  * Noether [1971] E. Noether, Transport theory and statistical physics 1, 186 (1971).
  * Netz and Orland [2000] R. R. Netz and H. Orland, The European Physical Journal E 1, 203 (2000).
  * Lau [2008] A. Lau, Physical Review E 77, 011502 (2008).
  * Brandyshev and Budkov [2023] P. E. Brandyshev and Y. A. Budkov, The Journal of chemical physics 158 (2023).
  * Blossey, Maggs, and Podgornik [2017] R. Blossey, A. Maggs, and R. Podgornik, Physical Review E 95, 060602 (2017).
  * Landau and Lifshitz [1971] L. Landau and E. Lifshitz, Pergamon Press, Oxford (1971).
  * Earman and Glymour [1978] J. Earman and C. Glymour, Archive for History of Exact Sciences , 291 (1978).
  * Weinberg [1995] S. Weinberg,  _The quantum theory of fields_ (Cambridge university press, 1995).
  * Weinberg [1972] S. Weinberg,  _Gravitation and cosmology: principles and applications of the general theory of relativity_ (1972).
  * Sprik [2021] M. Sprik, The Journal of Chemical Physics 155 (2021).
  * Lifshitz and Pitaevskii [2013] E. M. Lifshitz and L. P. Pitaevskii,  _Statistical physics: theory of the condensed state_ , Vol. 9 (Elsevier, 2013).
  * Buyukdagli [2020] S. Buyukdagli, The Journal of chemical physics 152, 014902 (2020).
  * Zinn-Justin [2002] J. Zinn-Justin,  _Quantum field theory and critical phenomena_ , Vol. 113 (Clarendon Press, Oxford, 2002).
  * Maggs and Podgornik [2016] A. Maggs and R. Podgornik, Soft matter 12, 1219 (2016).