Statistical field theory of mechanical stresses in Coulomb fluids: general covariant approach vs Noether's theorem

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–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–18]. Therefore, a first principles approach to derive the stress tensor of inhomogeneous Coulomb fluids would be beneficial for practical application purposes.

Progress has been made in this field, as demonstrated in two recent studies [3, 19]. In the first study [3], Budkov and Kolesnikov applied Noether's (first) theorem [20, 21] to the grand thermodynamic potential (GTP) of a Coulomb fluid as a functional of the electrostatic potential. They established a conservation law for the symmetric stress tensor, σ_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 ions. The authors extended this equilibrium condition to include cases in which 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 went beyond the conventional Derjaguin–Landau–Verwey–Overbeek expression to include different reference models of the fluid. In another study [19], 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 GTP. 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 investigate polyelectrolyte solutions constrained in a conducting slit nanopore and observed anomalies in the disjoining pressure and electric differential capacitance at small pore thicknesses.

It is important to mention several earlier studies in which the authors explored the derivation and explicit evaluation of electrostatic normal stress in Coulomb fluids confined in slit-like pores at both weak-coupling and strong-coupling limits [22–25]. Moreira and Netz utilized Monte Carlo simulations and statistical field theory, specifically on the level of one-loop approximation and strong coupling theory, to investigate the behavior of highly charged plates in the presence of multivalent counterions. Their study revealed a novel unbinding transition at equilibrium plate separation, where the pressure changes from attractive to repulsive. In another study [23], Dean and Horgan proposed a contact value theorem for Coulomb fluids in planar or filmlike geometries using a Hamiltonian field theoretical representation of the system. Furthermore, Jho et al examined the strong-coupling electrostatic interaction between two like-charged nanoparticles using the strong coupling theory and Monte Carlo simulations. In a recent study by Buyukdagli [25], a contact-value identity was derived by considering the field-theoretic partition function of an electrolyte confined between two anionic membranes. This identity holds true for a wide range of intramolecular solute structures and electrostatic coupling strengths.

Despite the progress made in deriving the local stress tensor within the mean-field theory, as well as normal electrostatic stress in slit-like pores even beyond the mean-field theory, there is currently a lack of clarity on how to compute all components of the stress tensor beyond the mean-field approximation.

Fluctuation corrections to the mean-field stress tensor should occur when considering the fluctuation corrections for the thermodynamic potential. The fluctuation corrections for the mean-field approximation of the thermodynamic potential are always nonlocal functionals [26, 27], even on a one-loop correction level, which makes the use of Noether's theorem, formulated for local functionals [21], quite problematic. Thus, it would be valuable to have a generalization of the Noether's theorem formalism [21] 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 could be the general covariant approach that we recently proposed in our paper [28]. Our approach was based on Noether's second theorem, which allowed us to derive the symmetric stress tensor for any model of an inhomogeneous liquid as a functional derivative of a GTP with respect to the metric tensor. It is important to note that the general covariant approach does not offer any advantages over Noether's first theorem when it comes to local functionals. However, it can be more advantageous for nonlocal functionals [29]. We have applied this approach to several phenomenological nonlocal models of inhomogeneous Coulomb fluids, such as the Cahn–Hilliard-like model [30, 31], the Bazant–Storey–Kornyshev model [11], and the Maggs–Podgornik–Blossey model [30], and have obtained the corresponding phenomenological stress tensors. It should be noted that this method is similar to the one used by Hilbert in the general theory of relativity to derive the energy–momentum tensor from the action functional [32, 33].

This study proposes two methods for deriving the stress tensor of inhomogeneous Coulomb fluids. The first method is based on the application of Noether's first theorem to the grand partition function, which is presented in functional integral form. The second method is based on the aforementioned general covariant approach in combination with the functional Legendre transform to obtain fluctuation corrections to the mean-field approximation of the stress tensor for a Coulomb fluid.

We start from the generating functional as the following functional integral [3] over the fluctuating electrostatic potential φ with the auxiliary function $\rho(\textbf{r})$:

$$\begin{equation} \Xi\left[\rho\right] = \int \frac{\mathcal{D}\phi}{C_0} \exp\left(-\frac{\mathcal{H}\left[\phi\right]+\left(\rho\phi\right)}{k_\mathrm{B}T}\right) \end{equation} \tag{ 1 }$$

where

$$\begin{equation} \mathcal{H}\left[\phi\right] = \int \textrm{d}\textbf{r} \left[\frac{\varepsilon}{2}\left(\nabla\phi\right)^2-P\left(\left\{\bar{\mu}_{\alpha}\right\}\right)\right], \end{equation} \tag{ 2 }$$

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

$$\begin{equation} C_0 = \exp\left[-\frac{1}{2}\mathrm{tr} \ln\left(-\frac{\varepsilon}{k_\mathrm{B}T}\Delta\right)\right], \end{equation} \tag{ 3 }$$

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

$$\begin{equation} \left(\rho\phi\right) = \int \textrm{d}\textbf{r} \rho\left(\textbf{r}\right)\phi\left(\textbf{r}\right). \end{equation} \tag{ 4 }$$

At $\rho(\textbf{r}) = 0$, the generating functional (1) transforms into the grand partition function of the Coulomb fluid obtained in [3]. We emphasize that this study deals only 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 polar solvents. Nevertheless, a generalization of the theory for the latter cases (models II and III) can be performed directly. We also emphasize that our study focuses on inhomogeneous Coulomb fluids, which involve the interaction between ions and surface external charges. These charges, although not included in the field-theoretic Hamiltonian, can be taken into account by incorporating them into the boundary conditions for the self-consistent field equations. Furthermore, the spatial inhomogeneity of the Coulomb fluid can also be attributed to the presence of external potential fields characterized by potentials $u_{\alpha}(\textbf{r})$.

The generating functional can be rewritten as follows:

$$\begin{equation} \Xi\left[\rho\right] = \exp\left(-\frac{W\left[\rho\right]}{k_\mathrm{B}T}\right), \end{equation} \tag{ 5 }$$

where $W[\rho]$ is a functional of the auxiliary function $\rho(\textbf{r})$. Thus, the expectation value can be obtained from the equation

$$\begin{equation} i\varphi\left(\textbf{r}\right) = \frac{\delta W\left[\rho\right]}{\delta \rho\left(\textbf{r}\right)}. \end{equation} \tag{ 6 }$$

Note that at ρ = 0, we have $\varphi(\textbf{r}) = -i\langle\phi(\textbf{r})\rangle$, where $\langle(..)\rangle$ means average over the Gibbs statistics of the Coulomb fluid with the field-theoretical Hamiltonian (2). The effective GTP can be derived from the following functional Legendre transform [26, 27, 34]:

$$\begin{equation} \Omega\left[\phi\right] = W\left[\rho\right]-\left(\rho\phi\right). \end{equation} \tag{ 7 }$$

Therefore,

$$\begin{equation} \frac{\delta\Omega\left[\phi\right]}{\delta\phi\left(\textbf{r}\right)} = \int \textrm{d}\textbf{r}^{\prime} \frac{\delta W\left[\rho\right]}{\delta\rho\left(\textbf{r}^{\prime}\right)}\frac{\delta\rho\left(\textbf{r}^{\prime}\right)}{\delta\phi\left(\textbf{r}\right)}-\int \textrm{d}\textbf{r}^{\prime}\frac{\delta\rho\left(\textbf{r}^{\prime}\right)}{\delta\phi\left(\textbf{r}\right)}\phi\left(\textbf{r}^{\prime}\right)-\rho\left(\textbf{r}\right) \end{equation} \tag{ 8 }$$

and with account of equation (6), we get

$$\begin{equation} \frac{\delta\Omega\left[i\varphi\right]}{\delta\varphi\left(\textbf{r}\right)} = -i\rho\left(\textbf{r}\right). \end{equation} \tag{ 9 }$$

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

$$\begin{equation} \frac{\delta\mathcal{H}\left[i\psi\right]}{\delta\psi\left(\textbf{r}\right)} = -i\rho\left(\textbf{r}\right). \end{equation} \tag{ 10 }$$

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

$$\begin{equation} \phi = i\psi+\eta, \end{equation} \tag{ 11 }$$

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

$$\begin{equation} \frac{\mathcal{H}\left[i\psi+\eta\right]}{k_\mathrm{B}T} = \frac{\mathcal{H}\left[i\psi\right]}{k_\mathrm{B}T}+ \frac{\varepsilon}{2k_\mathrm{B}T}\int \textrm{d}\textbf{r} \eta\left(\textbf{r}\right)\left(-\Delta+\varkappa^2\left(\textbf{r}\right)\right)\eta\left(\textbf{r}\right), \end{equation} \tag{ 12 }$$

where

$$\begin{equation} \varkappa^2\left(\textbf{r}\right) = U\left(\psi\left(\textbf{r}\right)\right) = \frac{1}{\varepsilon}\sum\limits_{\alpha\gamma}{q_\alpha q_\gamma}\frac{\partial^{2}P}{\partial\bar{\mu}_\alpha \partial\bar{\mu}_\gamma}, \end{equation} \tag{ 13 }$$

and calculating the Gaussian integral (1) over the fluctuations η, we arrive at

$$\begin{equation} W\left[\rho\right] = \mathcal{H}\left[i\psi\right]+\frac{\Omega_1\left[\psi\right]}{2\beta}+i\left(\rho\psi\right)+O\left(\frac{1}{\beta^2}\right) \end{equation} \tag{ 14 }$$

where the functional dependence of ψ on ρ is determined by equation (10), $\beta = 1/k_\mathrm{B}T$, and the following notation:

$$\begin{equation} \Omega_1\left[\psi\right] = \mathrm{tr}\ln\left(\frac{-\Delta+U}{-\Delta}\right) \end{equation} \tag{ 15 }$$

is introduced. Symbol $\mathrm{tr}(..)$ denotes the trace of operator. We would like to note that $1/\beta$ plays the role of Planck's constant $\hbar$ in the quantum field theory [34]. The occurrence of the multiplayer $1/\beta = k_\mathrm{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, ϕ. Note that in this study, we hold only the first-order terms on $1/\beta$.

Using equations (6) and (14), we can get

$$\begin{equation} i\varphi\left(\textbf{r}\right) = \int \textrm{d}\textbf{r}^{\prime}\left( \frac{\delta\mathcal{H}\left[i\psi\right]}{\delta \psi\left(\textbf{r}^{\prime}\right)}+i\rho\left(\textbf{r}^{\prime}\right)\right)\frac{\delta\psi\left(\textbf{r}^{\prime}\right)}{\delta \rho\left(\textbf{r}\right)}+i\psi\left(\textbf{r}\right)+\frac{i\chi\left(\textbf{r}\right)}{2\beta}, \end{equation} \tag{ 16 }$$

where

$$\begin{equation} i\chi\left(\textbf{r}\right) = \int \textrm{d}\textbf{r}^{\prime} \frac{\delta\Omega_1\left[\psi\right]}{\delta\psi\left(\textbf{r}^{\prime}\right)}\frac{\delta\psi\left(\textbf{r}^{\prime}\right)}{\delta \rho\left(\textbf{r}\right)}. \end{equation} \tag{ 17 }$$

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

$$\begin{equation} \psi\left(\textbf{r}\right) = \varphi\left(\textbf{r}\right)-\frac{\chi\left(\textbf{r}\right)}{2\beta}. \end{equation} \tag{ 18 }$$

Substituting this into equation (14) and expanding it in a series in χ, we get

$$\begin{equation} \begin{aligned} W\left[\rho\right] = \mathcal{H}\left[i\varphi\right]+ i\left(\rho\varphi\right)-\int \textrm{d}\textbf{r}\left( \frac{\delta\mathcal{H}\left[i\psi\right]}{\delta \psi\left(\textbf{r}\right)}+i\rho\left(\textbf{r}\right)\right)\frac{\chi\left(\textbf{r}\right)}{2\beta} +\frac{\Omega_1\left[\varphi\right]}{2\beta} \end{aligned} \end{equation} \tag{ 19 }$$

and using equation (10) again, we arrive at

$$\begin{equation} \begin{aligned} W\left[\rho\right] = \mathcal{H}\left[i\varphi\right]+ \frac{\Omega_1\left[\varphi\right]}{2\beta}+i\left(\rho\varphi\right). \end{aligned} \end{equation} \tag{ 20 }$$

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

$$\begin{equation} \Omega\left[i\varphi\right] = \mathcal{H}\left[i\varphi\right] +\frac{\Omega_1\left[\varphi\right]}{2\beta}. \end{equation} \tag{ 21 }$$

Equation (10) at ρ = 0 can be written as

$$\begin{equation} \varepsilon\Delta\psi = -\sum_{\alpha}q_{\alpha}\bar{c}_{\alpha}, \end{equation} \tag{ 22 }$$

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

A variation of Ω₁ is derived from the expression

$$\begin{equation} \frac{\delta \Omega_1\left[\varphi\right]}{\delta\varphi\left(\textbf{r}\right)} = \beta\varepsilon\frac{\partial U\left(\varphi\right)}{\partial\varphi\left(\textbf{r}\right)}G\left(\textbf{r},\textbf{r}|\varphi\right), \end{equation} \tag{ 23 }$$

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

$$\begin{equation} \varepsilon\beta\left(-\Delta+U\left(\varphi\right)\right)G\left(\textbf{r},\textbf{r}^{\prime}|\varphi\right) = \delta\left(\textbf{r}-\textbf{r}^{\prime}\right). \end{equation} \tag{ 24 }$$

Thus, the EL equation (9) at ρ = 0 has the form

$$\begin{equation} -\frac{\partial P}{\partial\varphi}+\varepsilon\Delta\varphi+\frac{\varepsilon}{2}\frac{\partial U\left(\varphi\right)}{\partial\varphi\left(\textbf{r}\right)}G\left(\textbf{r},\textbf{r}|\varphi\right) = 0. \end{equation} \tag{ 25 }$$

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

$$\begin{equation} \left(-\Delta+U\left(\psi\right)\right)\frac{\chi\left(\textbf{r}\right)}{\beta}-\frac{\partial U\left(\psi\right)}{\partial\psi\left(\textbf{r}\right)}G\left(\textbf{r},\textbf{r}|\psi\right) = 0. \end{equation} \tag{ 26 }$$

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

$$\begin{equation} \frac{\delta\psi\left(\textbf{r}^{\prime}\right)}{\delta\rho\left(\textbf{r}\right)} = i\beta G\left(\textbf{r},\textbf{r}^{\prime}|\psi\right). \end{equation} \tag{ 27 }$$

Thus, substituting equations (23) and (27) into equation (17), we get

$$\begin{equation} \chi\left(\textbf{r}\right) = \varepsilon\beta^2\int \textrm{d}\textbf{r}^{\prime} G\left(\textbf{r},\textbf{r}^{\prime}|\psi\right) G\left(\textbf{r}^{\prime},\textbf{r}^{\prime}|\psi\right)\frac{\partial U\left(\psi\right)}{\partial\psi\left(\textbf{r}^{\prime}\right)}. \end{equation} \tag{ 28 }$$

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

In this section, we formulate an approach based on a generalization of Noether's first theorem [3, 20, 21] to derive the stress tensor of Coulomb fluids from the grand partition function presented in the functional integral form, taking into account electrostatic field thermal fluctuations. The grand partition function is

$$\begin{equation} \Xi = \int \frac{\mathcal{D}\phi}{C_0} \exp\bigg(-\frac{\mathcal{H}\left[\phi\right])}{k_\mathrm{B}T}\bigg), \end{equation} \tag{ 29 }$$

where the functional $\mathcal{H}$ is determined by equation (2). First, let us consider the case of $u_{\alpha}(\textbf{r}) = 0$.

Let us perform the global infinitesimal shift transformation of the coordinates

$$\begin{equation} x^\prime_{k} = x_{k}+h_{k}, \end{equation} \tag{ 30 }$$

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

$$\begin{equation} \delta\ln\Xi = -\frac{1}{k_\mathrm{B}T}\bigg\langle\delta\mathcal{H}\left[\phi\right]\bigg\rangle = 0. \end{equation} \tag{ 31 }$$

Then, after some algebra [20, 28], we obtain

$$\begin{equation} \partial_{i}\bigg\langle\hat{\sigma}_{ik}h_{k}\bigg\rangle = \bigg\langle \frac{\delta\mathcal{H}\left[\phi\right]}{\delta\phi}h_{k}\partial_{k}\phi \bigg\rangle, \end{equation} \tag{ 32 }$$

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

$$\begin{equation} \hat{\sigma}_{ik} = \frac{\varepsilon}{2}\delta_{ik}\partial_j\phi\partial_j\phi -\varepsilon\partial_i\phi\partial_k\phi-P\delta_{ik}. \end{equation} \tag{ 33 }$$

Because h_k is the arbitrary constant infinitesimal vector, we can get

$$\begin{equation} \partial_{i}\langle\hat{\sigma}_{ik}\rangle = \bigg\langle\frac{\delta\mathcal{H}\left[\phi\right]}{\delta\phi}\partial_{k}\phi\bigg\rangle, \end{equation} \tag{ 34 }$$

where the expectation value is

$$\begin{equation} \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}\left[\phi\right])}{k_\mathrm{B}T}\bigg).\end{equation} \tag{ 35 }$$

Further, taking into account that [35]

$$\begin{equation} \bigg\langle\frac{\delta\mathcal{H}\left[\phi\right]}{\delta\phi\left(\textbf{r}\right)}\partial_{k}\phi\left(\textbf{r}\right)\bigg\rangle = \partial_{k}\frac{\delta\phi\left(\textbf{r}\right)}{\delta\phi\left(\textbf{r}\right)} = \partial_{k}\delta\left(\textbf0\right) = \lim_{\textbf{r}^{\prime}\rightarrow \textbf{r}}\left(\partial_{k}+\partial_{k}^{\prime}\right)\delta\left(\textbf{r}-\textbf{r}^{\prime}\right) = 0, \end{equation} \tag{ 36 }$$

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

$$\begin{equation} \partial_{i}\langle\hat{\sigma}_{ik}\rangle = 0. \end{equation} \tag{ 37 }$$

The latter expression, which is the local mechanical equilibrium condition with the average stress tensor $\langle\hat{\sigma}_{ik}\rangle$, represents another important result of this study.

In the case that external fields occur ($u_{\alpha}(\textbf{r})\neq 0$), the same calculations lead to the following mechanical equilibrium condition:

$$\begin{equation} \partial_{i}\langle\hat{\sigma}_{ik}\rangle-\sum\limits_{\alpha}\langle\hat{c}_{\alpha}\rangle\partial_{k}u_{\alpha} = 0, \end{equation} \tag{ 38 }$$

where $\hat{c}_{\alpha}(\textbf{r}) = \partial P/\partial\bar{\mu}_{\alpha}$ is the microscopic fluctuating ionic concentration.

Now, let us estimate the average stress tensor, $\sigma_{ik} = \langle\hat{\sigma}_{ik}\rangle$. For this purpose, we introduce the fluctuation, $\eta(\textbf{r})$, near the mean-field potential by

$$\begin{equation} \phi = i\psi+\eta. \end{equation} \tag{ 39 }$$

Taking into account that fluctuation η is the value of order $(k_\mathrm{B}T)^{\frac{1}{2}}$, substituting equation (39) into equation (33), we can show that the fluctuating stress tensor has the form with accuracy up to terms of order $k_\mathrm{B}T$

$$\begin{equation} \begin{aligned} \hat{\sigma}_{ik} & = \sigma^\mathrm{(MF)}_{ik}-\delta_{ik}\sum\limits_{\lambda}\mu_{\lambda}^{\left(1\right)}\bar{c}_{\lambda}+i\varepsilon\left(\delta_{ik}\partial_j\eta\partial_j\psi-\partial_i\eta\partial_k\psi-\partial_i\psi\partial_k\eta\right)+i\eta\delta_{ik} \frac{\partial P}{\partial \psi}\\ &\quad+\frac{\varepsilon}{2}U\left(\psi\right)\eta^2\delta_{ik}+\frac{\varepsilon}{2}\delta_{ik}\partial_j\eta\partial_j\eta-\varepsilon\partial_i\eta\partial_k\eta, \end{aligned} \end{equation} \tag{ 40 }$$

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 [22, 26, 27].

As obtained in [3], the mean-field approximation for the stress tensor is

$$\begin{equation} \sigma^\mathrm{(MF)}_{ik} = \varepsilon\partial_i\psi\partial_k\psi- \frac{\varepsilon}{2}\delta_{ik}\partial_j\psi\partial_j\psi-P\delta_{ik}. \end{equation} \tag{ 41 }$$

In what follows, utilizing equation (37), we will calculate the fluctuation corrections to the mean-field approximation (41).

Thus, to take into account corrections to the stress tensor up to the first order on $k_\mathrm{B}T$, we have to restrict ourselves to an expansion of $\mathcal{H}$ up to the third order on η. Then, using the EL equation

$$\begin{equation} \frac{\delta\mathcal{H}\left[i\psi\right]}{\delta\psi} = 0, \end{equation} \tag{ 42 }$$

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

$$\begin{align} \frac{\mathcal{H}\left[i\psi+\eta\right]}{k_\mathrm{B}T} & = \frac{\mathcal{H}\left[i\psi\right]}{k_\mathrm{B}T}+ \frac{\varepsilon}{2k_\mathrm{B}T}\int \textrm{d}\textbf{r} \eta\left(\textbf{r}\right)\left(-\Delta+U\left(\psi\right)\right)\eta\left(\textbf{r}\right)\nonumber\\ &\quad -\frac{i\varepsilon}{6k_\mathrm{B}T}\int \textrm{d}\textbf{r}\frac{\partial U\left(\psi\right)}{\partial\psi}\eta^3\left(\textbf{r}\right)+\frac{O\left(\eta^4\right)}{k_\mathrm{B}T}. \end{align} \tag{ 43 }$$

Then, we have

$$\begin{equation} \begin{aligned} \mathrm{e}^{-\beta\mathcal{H}\left[i\psi+\eta\right]} = \mathrm{e}^{-\beta\mathcal{H}\left[i\psi\right]}\mathrm{e}^{-\frac{1}{2}\left(\eta G^{-1}\eta\right)} \left(1+\frac{i\varepsilon\beta}{6}\int \textrm{d}\textbf{r}\frac{\partial U\left(\psi\right)}{\partial\psi}\eta^3\left(\textbf{r}\right)+O\left(\frac{1}{\beta}\right)\right), \end{aligned} \end{equation} \tag{ 44 }$$

where we have introduced the following short-hand notation:

$$\begin{equation} \left(\eta G^{-1}\eta\right) = \int \textrm{d}\textbf{r}\int \textrm{d}\textbf{r}^{\prime} \eta\left(\textbf{r}\right) G^{-1}\left(\textbf{r},\textbf{r}^{\prime}|\psi\right)\eta\left(\textbf{r}^{\prime}\right).\end{equation} \tag{ 45 }$$

Substituting equations (40) and (44) into equation (35) and integrating on η, we get the stress tensor with accuracy up to terms of order $1/\beta$

$$\begin{equation} \begin{aligned} \sigma_{ik} = \sigma^{\left(0\right)}_{ik}+\frac{\varepsilon}{2\beta}\left(\partial_i\chi\partial_k\psi+\partial_i\psi\partial_k\chi- \delta_{ik}\partial_j\chi\partial_j\psi\right)-\delta_{ik} \frac{\chi}{2\beta}\frac{\partial P}{\partial \psi}+\sigma^{\left(1\right)}_{ij}, \end{aligned} \end{equation} \tag{ 46 }$$

where

$$\begin{equation} \sigma^{\left(0\right)}_{ik} = \sigma^\mathrm{(MF)}_{ik}-\delta_{ik}\sum\limits_{\lambda}\mu_{\lambda}^{\left(1\right)}\bar{c}_{\lambda}, \end{equation} \tag{ 47 }$$

$$\begin{equation} \begin{aligned} \sigma^{\left(1\right)}\left(\textbf{r}\right) = \sigma^\mathrm{(cor)}_{ij}\left(\textbf{r}\right) = \frac{\varepsilon}{2}\left(U\left(\psi\right)G\left(\textbf{r},\textbf{r}|\psi\right) +\mathcal{D}_{kk}\left(\textbf{r}\right)\right)\delta_{ij}-\varepsilon \mathcal{D}_{ij}\left(\textbf{r}\right), \end{aligned} \end{equation} \tag{ 48 }$$

and following short-hand notations

$$\begin{equation} \mathcal{D}_{ij}\left(\textbf{r}\right) = \lim_{\textbf{r}^{\prime}\rightarrow \textbf{r}} \partial_i\partial_j^{\prime} G\left(\textbf{r},\textbf{r}^{\prime}|\psi\right), \end{equation} \tag{ 49 }$$

$$\begin{equation} \chi\left(\textbf{r}\right) = \frac{\varepsilon\beta^2}{3}\int \textrm{d}\textbf{r}^{\prime}\frac{\partial{U\left(\psi\left(\textbf{r}^{\prime}\right)\right)}}{\partial\psi\left(\textbf{r}^{\prime}\right)}\bigg\langle\eta\left(\textbf{r}\right)\eta^3\left(\textbf{r}^{\prime}\right) \bigg\rangle_{0} \end{equation} \tag{ 50 }$$

have been introduced. Note that we removed the terms that include the bare Green function, $G_{0}(\textbf{r},\textbf{r}^{\prime})$, because they pertain to identically divergenceless tensor and do not contribute to mechanical forces.

We can determine the average over the Gaussian measure using the expression

$$\begin{equation} \left<\left(\cdot\right)\right> _{0} = \int \frac{\mathcal{D}\eta}{C_1} \mathrm{e}^{-\frac{1}{2}\left(\eta G^{-1}\eta\right)} \left(\cdot\right), \end{equation} \tag{ 51 }$$

where C₁ is defined as

$$\begin{equation} C_1 = \int \mathcal{D}\eta \mathrm{e}^{-\frac{1}{2}\left(\eta G^{-1}\eta\right)}. \end{equation} \tag{ 52 }$$

By applying Wick's theorem [36] from equation (50), we can obtain the expression (28), as stated previously. Equation (48), which is the main result of this study, 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. This is known as the correlation stress tensor.

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

$$\begin{equation} F_{i} = \oint\limits_{\mathcal{A}}\sigma_{ik}n_{k} \mathrm{d}\mathcal{A}, \end{equation} \tag{ 53 }$$

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

Now, we discuss an alternative, more sophisticated approach to obtaining the total stress tensor from the GTP derived preceding equation (21). This approach is based on the general covariant methodology presented in our recent study [28]. As mentioned in Introduction, in this study, this approach does not have any advantages relative to the one based on Noether's theorem. However, we would like to consider it from a pedagogical standpoint.

In this approach, the stress tensor can be obtained using the following expression:

$$\begin{equation} \sigma_{ik} = \frac{2}{\sqrt{g}}\frac{\delta \Omega}{\delta g_{ik}}\bigg{|}_{g_{ik} = \delta_{ik}}, \end{equation} \tag{ 54 }$$

where g_ij is the metric tensor, $g = \det{g_{ij}}$—its determinant, and Ω is the GTP (21) obtained above within the functional Legendre transformation approach.

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

$$\begin{equation} \Omega\left[i\varphi\right] = \mathcal{H}\left[i\varphi\right]+\mathcal{H}^{^{\prime}}\left[\varphi\right],\quad \mathcal{H}^{^{\prime}}\left[\varphi\right] = \frac{\Omega_1\left[\varphi\right]}{2\beta}, \end{equation} \tag{ 55 }$$

where

$$\begin{equation} \mathcal{H} = \int \textrm{d}\textbf{r} \sqrt{g} \left[-\frac{\varepsilon}{2}g^{ij}\partial_{i}\varphi\partial_{j}\varphi-P\left(\left\{\bar{\mu}_{\alpha}\right\}\right)\right], \end{equation} \tag{ 56 }$$

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 that for simplicity, $u_{\alpha} = 0$.

The one-loop correction is

$$\begin{equation} \mathcal{H}^{^{\prime}}\left[\varphi\right] = \frac{k_\mathrm{B}T}{2}\mathrm{tr}\ln\left(\frac{-\Delta+U\left(\varphi\right)}{-\Delta}\right), \end{equation} \tag{ 57 }$$

where the Laplacian can be written in general covariant form [32, 37]

$$\begin{equation} \Delta f = \frac{1}{\sqrt{g}}\partial_{i}\left(\sqrt{g}g^{ij}\partial_{j}f\right). \end{equation} \tag{ 58 }$$

General coordinate transformations lead to

$$\begin{equation} \textrm{d}\textbf{r}\rightarrow J\textrm{d}\textbf{r},\quad \sqrt{g}\rightarrow J^{-1}\sqrt{g}, \end{equation} \tag{ 59 }$$

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

$$\begin{equation} \int \textrm{d}\textbf{r}\sqrt{g\left(\textbf{r}\right)}\delta\left(\textbf{r}\right) = 1. \end{equation} \tag{ 60 }$$

Let us consider an infinitesimal transformation of the metric tensor

$$\begin{equation} g_{ij}\left(\textbf{r}\right)\rightarrow g_{ij}\left(\textbf{r}\right)+\delta g_{ij}\left(\textbf{r}\right). \end{equation} \tag{ 61 }$$

Then, we have [32]

$$\begin{equation} \delta\left(\sqrt{g}\right) = \frac{1}{2}\sqrt{g}g^{ij}\delta g_{ij},\quad \delta g^{ij} = -g^{im}g^{jn}\delta g_{mn}.\end{equation} \tag{ 62 }$$

The zeroth-order term is

$$\begin{equation} \sigma^{\left(0\right)}_{ik} = \frac{2}{\sqrt{g}}\frac{\delta\mathcal{H}\left[i\varphi\right]}{\delta g_{ik}}\bigg|_{g_{ik} = \delta_{ik}}, \end{equation} \tag{ 63 }$$

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

$$\begin{equation} \sigma^{\left(0\right)}_{ik} = \varepsilon\partial_i\varphi\partial_k\varphi- \frac{\varepsilon}{2}\delta_{ik}\partial_j\varphi\partial_j\varphi-P\delta_{ik}. \end{equation} \tag{ 64 }$$

However, it should not be confused with the mean-field approximation (41) because the electrostatic potential ϕ satisfies the equation of EL (25), taking into account the electrostatic correlations. It can be shown that the divergence of this tensor is

$$\begin{equation} \partial_i\sigma^{\left(0\right)}_{ik} = \left(-\frac{\partial P\left(\varphi\right)}{\partial\varphi}+\varepsilon\Delta\varphi\right)\partial_k\varphi.\end{equation} \tag{ 65 }$$

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

$$\begin{equation} \sigma^{\left(0\right)}_{ik} = \sigma^\mathrm{(MF)}_{ik}-\delta_{ik}\sum\limits_{\lambda}\mu_{\lambda}^{\left(1\right)}\bar{c}_{\lambda}+\frac{\varepsilon}{2\beta}\left(\partial_i\chi\partial_k\psi+\partial_i\psi\partial_k\chi- \delta_{ik}\partial_j\chi\partial_j\psi\right)-\delta_{ik} \frac{\chi}{2\beta}\frac{\partial P}{\partial \psi}. \end{equation} \tag{ 66 }$$

The third and fourth terms on the right-hand side of equation (66) determine the contribution to the total stress tensor arising from the mismatch between the electrostatic potential of the mean field, ψ, and its expectation value, ϕ.

Thus, the stress tensor is

$$\begin{equation} \sigma_{ik} = \sigma^{\left(0\right)}_{ik}+\sigma^{\left(1\right)}_{ik}, \end{equation} \tag{ 67 }$$

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

$$\begin{equation} \sigma^{\left(1\right)}_{ik} = \frac{2}{\sqrt{g}}\frac{\delta\mathcal{H}^{^{\prime}}}{\delta g_{ik}}\bigg|_{g_{ik} = \delta_{ik}}. \end{equation} \tag{ 68 }$$

To perform the calculation in equation (68), we have to introduce a general covariant definition of the trace, i.e.

$$\begin{equation} \mathrm{tr}\left(A\right) = \int \mathrm{d}\textbf{r} \sqrt{g\left(\textbf{r}\right)}A\left(\textbf{r},\textbf{r}\right). \end{equation} \tag{ 69 }$$

For composition of two integral operators,

$$\begin{equation} C = AB, \end{equation} \tag{ 70 }$$

we can introduce the rule

$$\begin{equation} C\left(\textbf{r},\textbf{r}^{\prime}\right) = \int \textrm{d}\textbf{r}^{^{\prime\prime}} \sqrt{g\left(\textbf{r}^{^{\prime\prime}}\right)}A\left(\textbf{r},\textbf{r}^{^{\prime\prime}}\right)B\left(\textbf{r}^{^{\prime\prime}},\textbf{r}^{\prime}\right). \end{equation} \tag{ 71 }$$

The action of the operator A on a function $f(\textbf{r})$ is determined by

$$\begin{equation} Af\left(\textbf{r}\right) = \int \textrm{d}\textbf{r}^{\prime} \sqrt{g\left(\textbf{r}^{\prime}\right)} A\left(\textbf{r},\textbf{r}^{\prime}\right)f\left(\textbf{r}^{\prime}\right). \end{equation} \tag{ 72 }$$

The trace variation is

$$\begin{equation} \delta \mathrm{tr}\left(A\right) = \int \textrm{d}\textbf{r} \sqrt{g\left(\textbf{r}\right)}\delta A\left(\textbf{r},\textbf{r}\right)+\frac{1}{2}\int \textrm{d}\textbf{r} \sqrt{g\left(\textbf{r}\right)}g^{ij}\left(\textbf{r}\right)\delta g_{ij}\left(\textbf{r}\right) A\left(\textbf{r},\textbf{r}\right), \end{equation} \tag{ 73 }$$

which can be rewritten as

$$\begin{equation} \delta \mathrm{tr}\left(A\right) = \mathrm{tr}\left(\bar{\delta} A\right), \end{equation} \tag{ 74 }$$

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

$$\begin{equation} \bar{\delta} A\left(\textbf{r},\textbf{r}^{\prime}\right) = \delta A\left(\textbf{r},\textbf{r}^{\prime}\right)+\frac{1}{2} A\left(\textbf{r},\textbf{r}^{\prime}\right)g^{ij}\left(\textbf{r}^{\prime}\right)\delta g_{ij}\left(\textbf{r}^{\prime}\right). \end{equation} \tag{ 75 }$$

Thus, we have

$$\begin{equation} \delta \mathrm{tr}\left(A\right) = \int \mathrm{d}\textbf{r} \sqrt{g\left(\textbf{r}\right)}\bar{\delta} A\left(\textbf{r},\textbf{r}\right). \end{equation} \tag{ 76 }$$

Using equation (62), we can show that the operator Δ is transformed as

$$\begin{equation} \Delta\rightarrow \Delta+\Delta^{\prime}, \end{equation} \tag{ 77 }$$

where

$$\begin{align} \Delta^{\prime} f & = \frac{1}{2\sqrt{g}}\partial_{i}\left(\sqrt{g}g^{mn}\delta g_{mn}g^{ij}\partial_{j}f\right)- \frac{1}{2\sqrt{g}}g^{mn}\delta g_{mn}\partial_{i}\left(\sqrt{g}g^{ij}\partial_{j}f\right)\nonumber\\ &\quad-\frac{1}{\sqrt{g}}\partial_{i}\left(\sqrt{g}\delta g_{mn}g^{im}g^{jn}\partial_{j}f\right). \end{align} \tag{ 78 }$$

The functional (57) can be rewritten as

$$\begin{equation} \mathcal{H}^{^{\prime}} = \frac{k_\mathrm{B}T}{2}\mathrm{tr}\ln\left(I+\frac{\varepsilon G_{0}U}{k_\mathrm{B}T}\right), \end{equation} \tag{ 79 }$$

where operator G₀ is determined by equation

$$\begin{equation} -\frac{\varepsilon}{k_\mathrm{B}T}\Delta G_{0} = I, \end{equation} \tag{ 80 }$$

where I is the unity operator. Varying both sides of equation (80) on metric tensor, we arrive at (see appendix A)

$$\begin{equation} -\frac{\varepsilon}{k_\mathrm{B}T}\left(\Delta^{\prime} G_{0}+\Delta\bar{\delta} G_{0}\right) = 0 \end{equation} \tag{ 81 }$$

that in turn yields

$$\begin{equation} \bar{\delta} G_{0} = \frac{\varepsilon}{k_\mathrm{B}T}G_{0}\Delta^{\prime} G_{0}. \end{equation} \tag{ 82 }$$

Thus, the functional variation of the fluctuation correction is

$$\begin{equation} \delta\mathcal{H}^{^{\prime}} = \frac{\varepsilon}{2}\mathrm{tr}\left[\left(I+\frac{\varepsilon G_{0}U}{k_\mathrm{B}T}\right)^{-1} \bar{\delta} G_{0}U\right]. \end{equation} \tag{ 83 }$$

Therefore, we obtain

$$\begin{equation} \delta\mathcal{H}^{^{\prime}} = \frac{\varepsilon}{2}\mathrm{tr}\left[\left(I+\frac{\varepsilon G_{0}U}{k_\mathrm{B}T}\right)^{-1} \frac{\varepsilon}{k_\mathrm{B}T}G_{0}\Delta^{\prime} G_{0}U\right]. \end{equation} \tag{ 84 }$$

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

$$\begin{equation} \frac{\varepsilon}{k_\mathrm{B}T}\left(-\Delta+U\right)G = I. \end{equation} \tag{ 85 }$$

Then, using equations (80) and (85), we can get

$$\begin{equation} \left(I+\frac{\varepsilon G_{0}U}{k_\mathrm{B}T}\right)^{-1} = GG_{0}^{-1}. \end{equation} \tag{ 86 }$$

Thus, the functional variation is

$$\begin{equation} \delta\mathcal{H}^{^{\prime}} = \frac{\varepsilon^{2}}{2k_\mathrm{B}T}\mathrm{tr}\left(G\Delta^{\prime} G_0U\right). \end{equation} \tag{ 87 }$$

Using equation (86) and taking into account equation (80), we arrive at

$$\begin{equation} \frac{\varepsilon G_{0}U}{k_\mathrm{B}T} = G_{0}G^{-1}-I. \end{equation} \tag{ 88 }$$

Thus, by using the invariance of the trace with respect to cyclic permutations of operators in equation (87), we can obtain the following:

$$\begin{equation} \begin{aligned} \delta\mathcal{H}^{^{\prime}} = -\frac{\varepsilon}{2}\mathrm{tr}\left(\Delta^{\prime}\bar{G}\right), \end{aligned} \end{equation} \tag{ 89 }$$

where

$$\begin{equation} \bar{G} = G-G_0. \end{equation} \tag{ 90 }$$

Expression (78) can be rewritten as

$$\begin{equation} \begin{aligned} \Delta^{\prime} f = Df- \frac{1}{2}g^{mn}\delta g_{mn}\Delta f, \end{aligned} \end{equation} \tag{ 91 }$$

where we have introduced the following auxiliary differential operator:

$$\begin{equation} \begin{aligned} Df = \frac{1}{2\sqrt{g}}\partial_{i}\left(\sqrt{g}g^{mn}\delta g_{mn}g^{ij}\partial_{j}f\right) -\frac{1}{\sqrt{g}}\partial_{i}\left(\sqrt{g}\delta g_{mn}g^{im}g^{jn}\partial_{j}f\right). \end{aligned} \end{equation} \tag{ 92 }$$

Then, we have

$$\begin{equation} \begin{aligned} \delta\mathcal{H}^{^{\prime}} = -\frac{\varepsilon}{2}\mathrm{tr}\left(D\bar{G}\right)+\frac{\varepsilon}{4}\mathrm{tr}\left(g^{mn}\delta g_{mn}\Delta \bar{G}\right), \end{aligned} \end{equation} \tag{ 93 }$$

where

$$\begin{equation} \begin{aligned} \mathrm{tr}\left(DG\right) = \int \textrm{d}\textbf{r} \int \textrm{d}\textbf{r}^{\prime} \sqrt{g\left(\textbf{r}^{\prime}\right)}G\left(\textbf{r}^{\prime},\textbf{r}\right)\sqrt{g\left(\textbf{r}\right)}D_\textbf{r}\delta\left(\textbf{r}-\textbf{r}^{\prime}\right) \end{aligned} \end{equation} \tag{ 94 }$$

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

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

$$\begin{align} \mathrm{tr}\left(DG\right) = -\int \textrm{d}\textbf{r}^{\prime} \int \textrm{d}\textbf{r} \sqrt{g\left(\textbf{r}^{\prime}\right)}\partial_{i}G\left(\textbf{r}^{\prime},\textbf{r}\right)\sqrt{g\left(\textbf{r}\right)}\delta g_{mn}\left(\textbf{r}\right)\partial_{j}\delta\left(\textbf{r}-\textbf{r}^{\prime}\right)\gamma^{ijmn}\left(\textbf{r}\right), \end{align} \tag{ 95 }$$

where we have introduced the following tensor function:

$$\begin{equation} \gamma^{ijmn}\left(\textbf{r}\right) = \frac{1}{2}g^{mn}\left(\textbf{r}\right)g^{ij}\left(\textbf{r}\right) -g^{im}\left(\textbf{r}\right)g^{jn}\left(\textbf{r}\right). \end{equation} \tag{ 96 }$$

Thus, we have

$$\begin{align} \int \textrm{d}\textbf{r}^{\prime} \sqrt{g\left(\textbf{r}^{\prime}\right)}\partial_{i}G\left(\textbf{r}^{\prime},\textbf{r}\right)\partial_{j}\delta\left(\textbf{r}-\textbf{r}^{\prime}\right) = \lim_{\textbf{r}^{\prime}\rightarrow \textbf{r}}\partial_{j}\int \textrm{d}\textbf{r}^{^{\prime\prime}} \sqrt{g\left(\textbf{r}^{^{\prime\prime}}\right)}\partial_{i}^{\prime} G\left(\textbf{r}^{^{\prime\prime}},\textbf{r}^{\prime}\right)\delta\left(\textbf{r}-\textbf{r}{^{\prime\prime}}\right) \end{align} \tag{ 97 }$$

so that

$$\begin{equation} \begin{aligned} \int \textrm{d}\textbf{r}^{\prime} \sqrt{g\left(\textbf{r}^{\prime}\right)}\partial_{i}G\left(\textbf{r}^{\prime},\textbf{r}\right)\partial_{j}\delta\left(\textbf{r}-\textbf{r}^{\prime}\right) = \lim_{\textbf{r}^{\prime}\rightarrow \textbf{r}}\partial_{j}\partial_{i}^{\prime} G\left(\textbf{r},\textbf{r}^{\prime}\right). \end{aligned} \end{equation} \tag{ 98 }$$

Therefore, we obtain

$$\begin{equation} \mathrm{tr}\left(DG\right) = -\int \textrm{d}\textbf{r} \sqrt{g\left(\textbf{r}\right)}\delta g_{mn}\left(\textbf{r}\right)\gamma^{ijmn}\left(\textbf{r}\right)\lim_{\textbf{r}^{\prime}\rightarrow \textbf{r}}\partial_{j}\partial_{i}^{\prime} G\left(\textbf{r},\textbf{r}^{\prime}\right). \end{equation} \tag{ 99 }$$

The second term in equation (93) is

$$\begin{align} \mathrm{tr}\left(g^{mn}\delta g_{mn}\Delta G\right) = \int \textrm{d}\textbf{r}^{^{\prime} {mn}}\left(\textbf{r}^{\prime}\right)\delta g_{mn}\left(\textbf{r}^{\prime}\right)\sqrt{g\left(\textbf{r}^{\prime}\right)}\int \textrm{d}\textbf{r}\sqrt{g\left(\textbf{r}\right)}G\left(\textbf{r},\textbf{r}^{\prime}\right)\Delta_{\textbf{r}^{\prime}}\delta\left(\textbf{r}^{\prime}-\textbf{r}\right). \end{align} \tag{ 100 }$$

Using the expression

$$\begin{align} \int \textrm{d}\textbf{r}\sqrt{g\left(\textbf{r}\right)}G\left(\textbf{r},\textbf{r}^{\prime}\right)\Delta_{\textbf{r}^{\prime}}\delta\left(\textbf{r}^{\prime}-\textbf{r}\right) = \lim_{\textbf{r}^{^{\prime\prime}}\rightarrow \textbf{r}^{\prime}}\Delta_{\textbf{r}^{\prime}}\int \textrm{d}\textbf{r}\delta\left(\textbf{r}^{\prime}-\textbf{r}\right)\sqrt{g\left(\textbf{r}\right)}G\left(\textbf{r},\textbf{r}^{^{\prime\prime}}\right), \end{align} \tag{ 101 }$$

we can write

$$\begin{equation} \mathrm{tr}\left(g^{mn}\delta g_{mn}\Delta G\right) = \\ \int \textrm{d}\textbf{r}^{\prime} g^{mn}\left(\textbf{r}^{\prime}\right)\delta g_{mn}\left(\textbf{r}^{\prime}\right)\sqrt{g\left(\textbf{r}^{\prime}\right)}\lim_{\textbf{r}^{^{\prime\prime}}\rightarrow \textbf{r}^{\prime}}\left(\Delta_{\textbf{r}^{\prime}}G\left(\textbf{r}^{\prime},\textbf{r}^{^{\prime\prime}}\right)\right). \end{equation} \tag{ 102 }$$

Then, using equation (68), we can get

$$\begin{align} \sigma^{\left(1\right)}_{ij}\left(\textbf{r}\right) = \frac{\varepsilon}{2}\lim_{\textbf{r}^{\prime}\rightarrow \textbf{r}}\left(\delta_{ij}\partial_k\partial_k \bar{G}\left(\textbf{r},\textbf{r}^{\prime}\right)+\delta_{ij}\partial_k\partial_k^{\prime} \bar{G}\left(\textbf{r},\textbf{r}^{\prime}\right)-\partial_i\partial_j^{\prime} \bar{G}\left(\textbf{r},\textbf{r}^{\prime}\right)-\partial_j\partial_i^{\prime} \bar{G}\left(\textbf{r},\textbf{r}^{\prime}\right)\right). \end{align} \tag{ 103 }$$

The latter can be expressed in more compact form

$$\begin{equation} \begin{aligned} \sigma^{\left(1\right)}_{ij}\left(\textbf{r}\right) = \frac{\varepsilon}{2}\lim_{\textbf{r}^{\prime}\rightarrow \textbf{r}}\hat{D}_{ij}\bar{G}\left(\textbf{r},\textbf{r}^{\prime}\right), \end{aligned} \end{equation} \tag{ 104 }$$

where the following differential operator:

$$\begin{align} \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} \end{align} \tag{ 105 }$$

has been introduced.

The divergence of this tensor is

$$\begin{align} \partial_{i}\sigma^{(1)}_{ij}(\textbf{r}) = \frac{\varepsilon}{2}\lim_{\textbf{r}^{\prime} \rightarrow \textbf{r}} &\big[(\partial_{i}+\partial_{i}^{\prime})(\delta_{ij}\partial_k\partial_k \bar{G}(\textbf{r},\textbf{r}^{\prime})+\delta_{ij}\partial_k\partial_k^{\prime}\bar{G}(\textbf{r},\textbf{r}^{\prime})\nonumber\\ &\quad -\partial_i\partial_j^{\prime}\bar{G}(\textbf{r},\textbf{r}^{\prime})-\partial_j\partial_i^{\prime}\bar{G}(\textbf{r},\textbf{r}^{\prime}))\big]. \end{align} \tag{ 106 }$$

Therefore, after some algebra, we can get

$$\begin{align} \partial_{i}\sigma^{\left(1\right)}_{ik}\left(\textbf{r}\right) = \frac{\varepsilon}{2}\lim_{\textbf{r}^{\prime}\rightarrow \textbf{r}} \left[\partial_{k}\partial_{i}\partial_{i}\bar{G}\left(\textbf{r},\textbf{r}^{\prime}\right)-\partial_{i}^{\prime}\partial_{i}^{\prime}\partial_{k}\bar{G}\left(\textbf{r},\textbf{r}^{\prime}\right)\right], \end{align} \tag{ 107 }$$

and using equations (80) and (85), we rewrite equation (107) in the form

$$\begin{equation} \begin{aligned} \partial_{i}\sigma^{\left(1\right)}_{ik}\left(\textbf{r}\right) = \frac{\varepsilon}{2}G\left(\textbf{r},\textbf{r}\right)\partial_{k}U\left(\varphi\right), \end{aligned} \end{equation} \tag{ 108 }$$

where

$$\begin{equation} \partial_{k}U\left(\varphi\right) = \frac{\partial U\left(\varphi\right)}{\partial\varphi}\partial_{k}\varphi. \end{equation} \tag{ 109 }$$

Thus, using equations (65) and (108), we obtain

$$\begin{equation} \partial_{i}\sigma_{ik}\left(\textbf{r}\right) = \frac{\delta\Omega\left[\varphi\right]}{\delta\varphi\left(\textbf{r}\right)}\partial_{k}\varphi\left(\textbf{r}\right) \end{equation} \tag{ 110 }$$

where

$$\begin{equation} \frac{\delta\Omega\left[\varphi\right]}{\delta\varphi\left(\textbf{r}\right)} = -\frac{\partial P\left(\varphi\right)}{\partial\varphi}+\varepsilon\Delta\varphi+\frac{\varepsilon}{2}\frac{\partial U\left(\varphi\right)}{\partial\varphi\left(\textbf{r}\right)}G\left(\textbf{r},\textbf{r}\right). \end{equation} \tag{ 111 }$$

Thereby, the local mechanical equilibrium condition [3, 20, 38],

$$\begin{equation} \partial_{i}\sigma_{ik} = 0, \end{equation} \tag{ 112 }$$

is satisfied if ϕ is a solution of the EL equation (25). Thus, we have demonstrated the exactness of equation (110) with respect to the terms of any order in $1/\beta$ for the functional (55). However, as the functional (55) 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 (103) can be rewritten in the form

$$\begin{equation} \begin{aligned} \sigma^{\left(1\right)}_{ij}\left(\textbf{r}\right) = \sigma^\mathrm{(cor)}_{ij}\left(\textbf{r}\right) = \frac{\varepsilon}{2}\left(U\left(\psi\right)Q\left(\textbf{r}\right) +\mathcal{D}_{kk}\left(\textbf{r}\right)\right)\delta_{ij}-\varepsilon \mathcal{D}_{ij}\left(\textbf{r}\right), \end{aligned} \end{equation} \tag{ 113 }$$

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

$$\begin{equation} Q\left(\textbf{r}\right) = \bigg\langle \eta^2\left(\textbf{r}\right) \bigg\rangle = \lim_{\textbf{r}^{\prime}\rightarrow \textbf{r}}G\left(\textbf{r},\textbf{r}^{\prime}|\psi\right), \end{equation} \tag{ 114 }$$

$$\begin{equation} \mathcal{D}_{ij}\left(\textbf{r}\right) = \lim_{\textbf{r}^{\prime}\rightarrow \textbf{r}} \partial_i\partial_j^{\prime} G\left(\textbf{r},\textbf{r}^{\prime}|\psi\right) \end{equation} \tag{ 115 }$$

where the expectation values are calculated in the Gaussian approximation derived from the expansion near the mean-field configuration.

Thus, we arrive at the same expression for the correlation stress tensor (48), which has been obtained within the Noether theorem-based approach. Note that, as was done earlier, we removed the terms that include the bare Green function, $G_{0}(\textbf{r},\textbf{r}^{\prime})$, because they pertain to identically divergenceless tensors and do not contribute to mechanical forces.

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 Lifshitz [39, 40]. Nevertheless, our theory deals with purely classical electrostatic fluctuations, as opposed to the electromagnetic quantum fluctuations presented in the Lifshitz theory.

In this section, we would like to specify the correlation stress tensor obtained earlier for a practically important one-dimensional (1D) 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

$$\begin{equation} G\left(\textbf{r},\textbf{r}^{\prime}\right) = G\left({\boldsymbol{\rho}}-{\boldsymbol{\rho}}^{\prime}|z,z^{\prime}\right), \end{equation} \tag{ 116 }$$

where ρ 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

$$\begin{equation} G\left({\boldsymbol{\rho}}-{\boldsymbol{\rho}}^{\prime}|z,z^{\prime}\right) = \int \frac{\mathrm{d}^2\textbf{q}}{\left(2\pi\right)^2}\mathrm{e}^{-\mathrm{i}\textbf{q}\left({\boldsymbol{\rho}}-{\boldsymbol{\rho}}^{\prime}\right)}G\left(q|z,z^{\prime}\right), \end{equation} \tag{ 117 }$$

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

$$\begin{equation} \mathcal{D}_{xy}\left(z\right) = \lim_{z^{\prime}\rightarrow z}\int \frac{\mathrm{d}^2\textbf{q}}{\left(2\pi\right)^2}q_xq_yG\left(q|z,z^{\prime}\right), \end{equation} \tag{ 118 }$$

$$\begin{equation} \mathcal{D}_{xz}\left(z\right) = \lim_{z^{\prime}\rightarrow z}\int \frac{\mathrm{d}^2\textbf{q}}{\left(2\pi\right)^2}iq_x\partial_zG\left(q|z,z^{\prime}\right), \end{equation} \tag{ 119 }$$

$$\begin{equation} \mathcal{D}_{yz}\left(z\right) = \lim_{z^{\prime}\rightarrow z}\int \frac{\mathrm{d}^2\textbf{q}}{\left(2\pi\right)^2}iq_y\partial_zG\left(q|z,z^{\prime}\right). \end{equation} \tag{ 120 }$$

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

$$\begin{align} \mathcal{D}_{xy}\left(z\right) = \mathcal{D}_{xz}\left(z\right) = \mathcal{D}_{yz}\left(z\right) = 0. \end{align} \tag{ 121 }$$

Diagonal elements are

$$\begin{equation} \mathcal{D}_{xx}\left(z\right) = \int \frac{\mathrm{d}^2\textbf{q}}{\left(2\pi\right)^2}q_x^2Q\left(q,z\right),~\mathcal{D}_{yy}\left(z\right) = \int \frac{\mathrm{d}^2\textbf{q}}{\left(2\pi\right)^2}q_y^2Q\left(q,z\right), \end{equation} \tag{ 122 }$$

$$\begin{equation} \mathcal{D}_{zz}\left(z\right) = \int \frac{\mathrm{d}^2\textbf{q}}{\left(2\pi\right)^2}\mathcal{D}_{zz}\left(q,z\right), \end{equation} \tag{ 123 }$$

where

$$\begin{equation} \mathcal{D}_{zz}\left(q,z\right) = \lim_{z^{\prime}\rightarrow z}\partial_z\partial_{z^{\prime}}G\left(q|z,z^{\prime}\right), \end{equation} \tag{ 124 }$$

$$\begin{equation} Q\left(q,z\right) = \lim_{z^{\prime}\rightarrow z}G\left(q|z,z^{\prime}\right). \end{equation} \tag{ 125 }$$

Thus, the trace is

$$\begin{equation} \mathcal{D}_{kk}\left(z\right) = \int \frac{\mathrm{d}^2\textbf{q}}{\left(2\pi\right)^2}\left(q^2Q\left(q,z\right)+\mathcal{D}_{zz}\left(q,z\right)\right). \end{equation} \tag{ 126 }$$

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

$$\begin{equation} \varepsilon \beta \left(-\partial^2_z+q^2+\varkappa^2\left(z\right)\right)G\left(q|z,z^{\prime}\right) = \delta\left(z-z^{\prime}\right). \end{equation} \tag{ 127 }$$

The correlation stress tensor elements are defined by the following expressions:

$$\begin{equation} \begin{aligned} \sigma^\mathrm{(cor)}_{ij}\left(z\right) = \int \frac{\mathrm{d}^2\textbf{q}}{\left(2\pi\right)^2}\sigma^\mathrm{(cor)}_{ij}\left(q,z\right), \end{aligned} \end{equation} \tag{ 128 }$$

$$\begin{equation} \begin{aligned} \sigma^\mathrm{(cor)}_{xx}\left(q,z\right) = \frac{\varepsilon}{2}\left(Q\left(q,z\right)\varkappa^2\left(z\right) +\mathcal{D}_{zz}\left(q,z\right)\right), \end{aligned} \end{equation} \tag{ 129 }$$

$$\begin{equation} \begin{aligned} \sigma^\mathrm{(cor)}_{yy}\left(q,z\right) = \sigma^\mathrm{(cor)}_{xx}\left(q,z\right), \end{aligned} \end{equation} \tag{ 130 }$$

$$\begin{equation} \begin{aligned} \sigma^\mathrm{(cor)}_{zz}\left(q,z\right) = \frac{\varepsilon}{2}\left(Q\left(q,z\right)\left(q^2+\varkappa^2\left(z\right)\right) -\mathcal{D}_{zz}\left(q,z\right)\right).\end{aligned}\end{equation} \tag{ 131 }$$

Note that the formulas derived from formulas (128)–(131) in the case of Coulomb gas (point-like ions) for very wide pores yield expressions for the stresses in the bulk phase fluid, i.e. Debye–Hückel expression (see appendix B).

We developed a first-principles statistical field theory for analyzing macroscopic forces in spatially inhomogeneous Coulomb fluids. Additionally, we generalized Noether's first theorem to account for a fluctuating order parameter and derived the total stress tensor using this formulation. By employing this methodology, we obtained the stress tensor for the Coulomb fluid, which included both the previously derived mean-field stress tensor and fluctuation corrections at the one-loop correction level. These fluctuation corrections incorporate thermal fluctuations of the local electrostatic potential and field around the mean-field configuration. The correlation stress tensor, derived from these fluctuation corrections, reflects how the electrostatic correlation influences the local stresses in a nonuniform Coulomb fluid. Furthermore, by combining the previously formulated general covariant approach with the functional Legendre transform method, we reproduced the same results for the stress tensor.

The developed formalism is interesting for modeling not only inhomogeneous Coulomb fluids but also bulk fluids. In particular, it is intriguing to analyze how this formalism can predict the screening length for concentrated electrolyte solutions, considering the excluded volume of the ions, in comparison with other methods of statistical physics [41, 42]. This issue might be the subject of future research.

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

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.

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

$$\begin{equation} \int \mathrm{d}\textbf{r}\sqrt{\tilde{g}\left(\textbf{r}\right)}\tilde{\delta}\left(\textbf{r}-\textbf{r}^{\prime}\right)f\left(\textbf{r}\right) = f\left(\textbf{r}^{\prime}\right) \end{equation} \tag{ A1 }$$

where we consider metric variation as the type of function variation (not a variation under coordinate transformations)

$$\begin{equation} \tilde{g}_{ij}\left(\textbf{r}\right) = g_{ij}\left(\textbf{r}\right)+\delta g_{ij}\left(\textbf{r}\right), \end{equation} \tag{ A2 }$$

and then we can define the variation

$$\begin{equation} \delta^{\prime}\left(\textbf{r}-\textbf{r}^{\prime}\right) = \tilde{\delta}\left(\textbf{r}-\textbf{r}^{\prime}\right)-\delta\left(\textbf{r}-\textbf{r}^{\prime}\right). \end{equation} \tag{ A3 }$$

Then, using equation (60), we can show that the equality (A1) is fulfilled if the variation of delta function is determined by the following identity:

$$\begin{equation} \delta^{\prime}\left(\textbf{r}-\textbf{r}^{\prime}\right) = -\frac{1}{2}g^{ij}\left(\textbf{r}^{\prime}\right)\delta g_{ij}\left(\textbf{r}^{\prime}\right)\delta\left(\textbf{r}-\textbf{r}^{\prime}\right). \end{equation} \tag{ A4 }$$

Note that using known delta-function properties, we can write

$$\begin{equation} g^{ij}\left(\textbf{r}^{\prime}\right)\delta g_{ij}\left(\textbf{r}^{\prime}\right)\delta\left(\textbf{r}-\textbf{r}^{\prime}\right) = g^{ij}\left(\textbf{r}\right)\delta g_{ij}\left(\textbf{r}\right)\delta\left(\textbf{r}-\textbf{r}^{\prime}\right). \end{equation} \tag{ A5 }$$

It is evident that in this case

$$\begin{equation} \delta\left(\sqrt{g\left(\textbf{r}^{\prime}\right)}\delta\left(\textbf{r}-\textbf{r}^{\prime}\right)\right) = 0. \end{equation} \tag{ A6 }$$

Under transformations (A2), the operator G₀ is infinitesimally transformed as follows:

$$\begin{equation} \tilde{G_{0}} = G_{0}+\delta G_{0}, \end{equation} \tag{ A7 }$$

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

$$\begin{equation} -\frac{\varepsilon}{k_\mathrm{B}T}\left(\Delta_{\textbf{r}}+\Delta^\prime_{\textbf{r}}\right)\tilde{G_{0}}\left(\textbf{r},\textbf{r}^{\prime}\right) = \tilde{\delta}\left(\textbf{r}-\textbf{r}^{\prime}\right). \end{equation} \tag{ A8 }$$

Equation (39) can be written in the form

$$\begin{equation} -\frac{\varepsilon}{k_\mathrm{B}T}\Delta_{\textbf{r}}G_{0}\left(\textbf{r},\textbf{r}^{\prime}\right) = \delta\left(\textbf{r}-\textbf{r}^{\prime}\right). \end{equation} \tag{ A9 }$$

Then, subtracting equation (A9) from expression (A8) with account of equation (A4) and holding only first-order terms, we can get

$$\begin{equation} -\frac{\varepsilon}{k_\mathrm{B}T}\left(\Delta^\prime_{\textbf{r}}G_{0}\left(\textbf{r},\textbf{r}^{\prime}\right)+\Delta_{\textbf{r}}\delta G_{0}\left(\textbf{r},\textbf{r}^{\prime}\right)\right) = -\frac{1}{2}g^{ij}\left(\textbf{r}^{\prime}\right)\delta g_{ij}\left(\textbf{r}^{\prime}\right)\delta\left(\textbf{r}-\textbf{r}^{\prime}\right). \end{equation} \tag{ A10 }$$

Using equation (A9) again, we can rewrite this equality in the form

$$\begin{align} -\frac{\varepsilon}{k_\mathrm{B}T}\left(\Delta^\prime_{\textbf{r}} G_{0}\left(\textbf{r},\textbf{r}^{\prime}\right)+\Delta_{\textbf{r}}\delta G_{0}\left(\textbf{r},\textbf{r}^{\prime}\right)+\frac{1}{2}g^{ij}\left(\textbf{r}^{\prime}\right)\delta g_{ij}\left(\textbf{r}^{\prime}\right)\Delta_{\textbf{r}} G_{0}\left(\textbf{r},\textbf{r}^{\prime}\right)\right) = 0. \end{align} \tag{ A11 }$$

Using definition (75), equation (A11) can be reduced to equation (81).

Let us consider the following expressed diffeomorphic transformations:

$$\begin{equation} \begin{aligned} \tilde{x}^{k} = x^{k}+\xi^{k}\left(\textbf{r}\right),\\ \tilde{g}_{ij}\left(\tilde{\textbf{r}}\right) = g_{ij}\left(\textbf{r}\right)+\tilde{\delta} g_{ij}\left(\textbf{r}\right), \end{aligned} \end{equation} \tag{ A12 }$$

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

$$\begin{equation} \int d\tilde{\textbf{r}}_1\sqrt{\tilde{g}\left(\tilde{\textbf{r}}_1\right)}\tilde{\delta}\left(\tilde{\textbf{r}}_1-\tilde{\textbf{r}}_2\right)f\left(\tilde{\textbf{r}}_1\right) = f\left(\tilde{\textbf{r}}_2\right) \end{equation} \tag{ A13 }$$

we can show that the delta function is transformed under equation (A12) in accordance with the expression

$$\begin{equation} \tilde{\delta}\left(\tilde{\textbf{r}}_1-\tilde{\textbf{r}}_2\right) = \left(1-\partial_{k} \xi^{k}\left(\textbf{r}_1\right)- \frac{1}{2}g^{ij}\left(\textbf{r}_1\right)\tilde{\delta} g_{ij}\left(\textbf{r}_1\right)\right)\delta\left(\textbf{r}_1-\textbf{r}_2\right). \end{equation} \tag{ A14 }$$

If we put that

$$\begin{equation} \tilde{\delta}g^{ik}\left(\textbf{r}\right) = \xi^{i,k}\left(\textbf{r}\right)+\xi^{k,i}\left(\textbf{r}\right), \end{equation} \tag{ A15 }$$

where

$$\begin{equation} \xi^{i,k}\left(\textbf{r}\right) = g^{kj}\partial_{j} \xi^{i}\left(\textbf{r}\right), \end{equation} \tag{ A16 }$$

then in this special case

$$\begin{equation} \tilde{\delta}\left(\tilde{\textbf{r}}_1-\tilde{\textbf{r}}_2\right) = \delta\left(\textbf{r}_1-\textbf{r}_2\right). \end{equation} \tag{ A17 }$$

The delta function thus defined is invariant under a subgroup of transformations (A12) that correspond to equation (A15) (under standard diffeomorphisms). Our approach requires us to consider another subgroup of equation (A12) 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 (A4).

By analyzing the asymptotic behavior of the normal stresses for extremely large pore thicknesses, H, it is possible to determine whether the formulated formalism can accurately describe the stresses in the bulk Coulomb fluids. Specifically, the analysis would focus on the behavior of the normal stresses at $H\to \infty$. In this case, $\varkappa(z) = \kappa =\,$const (inverse screening length of bulk Coulomb fluid [43]) and

$$\begin{equation} G\left(q|z,z^{^{\prime}}\right)\simeq \frac{\exp\left(-\kappa_q|z-z^{^{\prime}}|\right)}{2\beta\varepsilon\kappa_q}, \end{equation} \tag{ B1 }$$

so that we have

$$\begin{equation} Q\left(q,z\right)\simeq\frac{1}{2\beta\varepsilon \kappa_q},~\mathcal{D}_{zz}\left(q,z\right)\simeq -\frac{\kappa_q}{2\beta\varepsilon}, \end{equation} \tag{ B2 }$$

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

$$\begin{equation} \sigma_{zz}^\mathrm{(cor)} = \frac{k_\mathrm{B}T}{2}\int\frac{\mathrm{d}^2\textbf{q}}{\left(2\pi\right)^2}\left(\kappa_{q}-q-\frac{\kappa^2}{2q}\right) = -\frac{k_\mathrm{B}T\kappa^3}{12\pi}, \end{equation} \tag{ B3 }$$

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 equation (66))

$$\begin{equation} \sigma_{zz}^{\left(0\right)} = -k_\mathrm{B}T\sum\limits_{\alpha}c_{\alpha}-\sum\limits_{\alpha}c_{\alpha}\mu_{\alpha}^{\left(1\right)}, \end{equation} \tag{ B4 }$$

where c_α 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

$$\begin{equation} P\left(\left\{\bar{\mu}_{\alpha}\right\}\right) = k_\mathrm{B}T\sum\limits_{\alpha}\Lambda_{\alpha}^{-3}\mathrm{e}^{\beta\bar{\mu}_{\alpha}}, \end{equation} \tag{ B5 }$$

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

$$\begin{align} \bar{c}_{\alpha}\left(\textbf{r}\right) = \frac{\delta \Omega}{\delta u_{\alpha}\left(\textbf{r}\right)} & = \Lambda_{\alpha}^{-3}\mathrm{e}^{\beta\bar{\mu}_{\alpha}}+\frac{\varepsilon}{2}\frac{\partial U\left(\varphi\right)}{\partial u_{\alpha}\left(\textbf{r}\right)}G\left(\textbf{r},\textbf{r}|\varphi\right)\nonumber\\ & =\Lambda_{\alpha}^{-3}\mathrm{e}^{\beta\bar{\mu}_{\alpha}}\left(1-\frac{q_{\alpha}^2}{2\left(k_\mathrm{B}T\right)^2}G\left(\textbf{r},\textbf{r}|\varphi\right)\right) \end{align} \tag{ B6 }$$

that yields

$$\begin{equation} \bar{\mu}_{\alpha} = k_\mathrm{B}T\ln\left(\Lambda_{\alpha}^3\bar{c}_{\alpha}\right)+\frac{q_{\alpha}^2}{2k_\mathrm{B}T}G\left(\textbf{r},\textbf{r}|\varphi\right). \end{equation} \tag{ B7 }$$

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

$$\begin{equation} \mu_{\alpha} = k_\mathrm{B}T\ln\left(\Lambda_{\alpha}^3{c}_{\alpha}\right)+\frac{q_{\alpha}^2}{2k_\mathrm{B}T}\left(G\left(0\right)-G_0\left(0\right)\right). \end{equation} \tag{ B8 }$$

We subtracted the infinite constant $q_{\alpha}^2G_0(0)/{2k_\mathrm{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, Green function is translation invariant, i.e. $G(\textbf{r},\textbf{r}^{^{\prime}}|0) = G(\textbf{r}-\textbf{r}^{^{\prime}})$. The second term in equation (B8) is the fluctuation correction to the ideal gas bulk chemical potential, which can be rewritten as follows:

$$\begin{equation} \mu_{\alpha}^{\left(1\right)} = \frac{q_{\alpha}^2}{2k_\mathrm{B}T}\int\frac{\mathrm{d}^2\textbf{q}}{\left(2\pi\right)^2}\left(G\left(q|z,z\right)-G_{0}\left(q|z,z\right)\right), \end{equation} \tag{ B9 }$$

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

$$\begin{equation} \sigma_{zz} = -P_\mathrm{o} = -k_\mathrm{B}T\sum\limits_{\alpha}c_{\alpha}+\frac{k_\mathrm{B}T\kappa^3}{24\pi}, \end{equation} \tag{ B10 }$$

where $P_\mathrm{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

$$\begin{equation} \sigma^\mathrm{(cor)}_{xx} = \sigma^\mathrm{(cor)}_{yy} = \frac{k_\mathrm{B}T}{4}\int\frac{\mathrm{d}^2\textbf{q}}{\left(2\pi\right)^2}\left(q-\frac{q^2}{\sqrt{q^2+\kappa^2}}-\frac{\kappa^2}{2q}\right) = -\frac{k_\mathrm{B}T\kappa^3}{12\pi}. \end{equation} \tag{ B11 }$$

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_\mathrm{o}\delta_{ij}$.