# Field theory of structured liquid dielectrics

Ralf Blossey  University of Lille, Unité de Glycobiologie Structurale et
Fonctionnelle (UGSF) CNRS UMR8576, 59000 Lille, France  Rudolf Podgornik
<EMAIL_ADDRESS>School of Physical Sciences and Kavli Institute for
Theoretical Sciences, University of Chinese Academy of Sciences, Beijing
100049, China and CAS Key Laboratory of Soft Matter Physics, Institute of
Physics, Chinese Academy of Sciences, Beijing 100190, China
and Wenzhou Institute of the University of Chinese Academy of Sciences,
Wenzhou, Zhejiang 325000, China  Also affiliated with: Department of
Theoretical Physics, Jožef Stefan Institute, SI-1000 Ljubljana, Slovenia and
Department of Physics, Faculty of Mathematics and Physics, University of
Ljubljana, SI-1000 Ljubljana, Slovenia

###### Abstract

We develop of a field-theoretic approach for the treatment of both the non-
local and the non-linear response of structured liquid dielectrics. Our
systems of interest are composed of dipolar solvent molecules and simple salt
cations and anions. We describe them by two independent order parameters, the
polarization field for the solvent and the charge density field for the ions,
including and treating the non-electrostatic part of the interactions
explicitly and consistently. We show how to derive functionals for the
polarization and the electrostatic field of increasingly finer scales and
solve the resulting mean-field saddle-point equations in the linear regime. We
derive criteria for the character of their solutions that depend on the
structural lengths and the polarity of the solvent. Our approach provides a
systematic way to derive generalized polarization theories.

field theory, soft matter electrostatics, polarization

##  I Introduction

Phenomenological functionals of polarization have featured prominently in the
development of the electrostatic theory of structured dielectrics starting
from the seminal contributions of Marcus in the 1950’s on electron transfer
theory based on local dielectric response theory [1, 2]. Aqueous dielectrics,
de rigueur in soft- and bio-matter systems [3], pose another challenge as the
dielectric response of water is strongly non-local and cannot be exhaustively
described by a local dielectric approximation [4]. This has been the driving
force behind the phenomenological non-local dielectric function methodology
developed extensively by Kornyshev and collaborators [5, 6, 4], as well as the
Landau-Ginzburg-type free energy functionals describing the non-local
dielectric response in spatially confined dielectrics [7, 8, 9, 10, 11], the
two approaches being of course equivalent in the bulk but differing in systems
with boundaries.

While it would be of course desirable to have more microscopic theories,
aqueous dielectrics have consistently eluded descriptions based on more
detailed explicit molecular interactions. As a successful example of a
microscopic non-local dielectric theory one can consider the statistical
mechanical formulation of the non-local dielectric properties of ice by
Onsager and Dupuis [12, 13], later further elaborated and generalized by Gruen
and Marčelja [14, 15]. Of all the methodologies applied to non-local
dielectrics this is the first to be found in the literature that implements
not only the non-local but also the non-linear dielectric response of a
system. It leads to two coupled equations for the electrostatic field and the
dielectric polarization field, and in general cannot be reduced to either the
non-local dielectric function methodology or the local Landau-Ginzburg
polarization functional. The coupling between the polarization and
electrostatic fields is provided by the structural interaction mediated by the
Bjerrum structural defects in the ideal hydrogen bonded lattice of ice. In its
absence the basic equations reduce to the standard dielectric response theory.
Being based on the crystalline lattice and hydrogen-bonding defects there is
no simple way to use Onsager-Dupuis theory straightforwardly for a liquid
aqueous dielectric.

An important area of application for the non-local dielectric response of
structured dielectrics was and remains the theory of hydration forces [16, 17]
where the dielectric response of ice [18] as well as the non-local
electrostatics of liquid water [19, 6] were applied to water confined between
substrate surfaces with the propensity for local ordering of vicinal water
molecules. This type of water structural interaction modeling has been applied
to solvation of phospholipid membranes [20], DNA molecules [21], linear
polysaccharides [22] or to the solvation of proteins [23, 24]. More recent
approaches extend the ideas of non-local dielectric response of confined water
[25] to order parameter description beyond the polarization [26], which are
able to describe additional features of confined water structuring [27, 28].
The order parameter description can be arguably viewed as a generalization of
the Landau-Ginzburg-type free energy functionals beyond the polarization
formulation [29], describing additional degrees of freedom pertinent to the
complicated hydrogen-bond phenomenology [30]. The hydration force theory also
sheds light on the fact that part of the water polarization is purely non-
electrostatic in nature, being a separate field, independent of the imposed
electrostatic field. This is in fact one of the most important conclusion of
the Onsager-Dupuis theory of the dielectric response of ice [13, 14].

Motivated by the Onsager-Dupuis approach we develop a systematic field theory
for the non-local as well as non-linear dielectric response of structured
liquid dielectrics based upon the electrostatic field and dielectric
polarization as separate order parameters. From the outset, our approach
decouples electrostatic and non-electrostatic structural contributions to
these order parameters, allowing for the consistent inclusion of the effects
that could not be implemented in previous approaches. On the linear level our
theory reduces to two equations for the electrostatic and polarization
potentials, which are coupled by the non-electrostatic interaction potentials
and assume a form completely consistent with the linearized Onsager-Dupuis
theory of ice. We explicitly discuss three examples of linear models differing
by the degree of structural detailed included, and show how they relate to
previously studied simplified models.

##  II Order parameters and microscopic interactions

The systems we consider are composed of three components: rigid rod-like
solvent with a dipolar moment (N), positive simple salt ions (+) and negative
simple salt ions (-). We first introduce two separate order parameters: charge
density and polarization (dipolar charge density). For the solvent (N)
component the order parameter is the polarization field defined as

$\hat{\cal\bf P}({\bf x})\equiv\sum_{(N)}p\,{\bf n}~{}\delta\left({\bf x}-{\bf
x}_{n}\right)$ (1)

where $p$ is the dipolar moment of the molecule, so that the bound charge
density field of the solvent material is

$\hat{\rho}_{(N)}({\bf x})\equiv p\sum_{N}{\bf
n}\cdot\mbox{\boldmath$\nabla$}~{}\delta\left({\bf x}-{\bf
x}_{n}\right)=\mbox{\boldmath$\nabla$}\\!\cdot\hat{\cal\bf P}({\bf x}),$ (2)

with dipoles considered to be point-like. For the simple salt component we can
define the total charge density $\hat{\rho}_{(\pm)}({\bf x})$ as the sum of
cation and anion density fields

$\hat{\rho}_{(\pm)}({\bf x})\equiv e\sum_{N_{+}}\delta({\bf x}-{\bf
x}_{n})-e\sum_{N_{-}}\delta({\bf x}-{\bf x}_{n}),$ (3)

where $e$ is the elementary charge of the salt ions, $N$ is the number of
dipoles while $N_{\pm}$ are the numbers of positive and negative ions,
respectively. The total charge density field is then

$\hat{\rho}({\bf x})=\hat{\rho}_{(\pm)}({\bf
x})+\mbox{\boldmath$\nabla$}\\!\cdot\hat{\cal\bf P}({\bf x}).$ (4)

We will now first write the interactions in terms of these microscopic order
parameters and then in the form of macroscopic collective coordinates, in
exactly the same way as this is done for standard Coulomb fluids and/or for
nematic Coulomb fluids [31]. The interactions in the system have the following
components: the non-electrostatic interactions in the harmonic approximation

${{\cal H}_{SR}}={\textstyle\frac{1}{2}}\int\\!\\!\int_{V}d{\bf x}d{\bf
x}^{\prime}~{}\hat{\cal\bf P}_{i}({\bf x})\tilde{u}_{ij}({\bf x}-{\bf
x}^{\prime})\hat{\cal\bf P}_{j}({\bf x}^{\prime}),$ (5)

where we assumed that the non-electrostatic tensorial part
$\tilde{u}_{ij}({\bf x}-{\bf x}^{\prime})$ is a short-range, non-electrostatic
potential; the electrostatic interactions are given by the standard Coulomb
form

${{\cal H}_{C}}={\textstyle\frac{1}{2}}\int\\!\\!\int_{V}d{\bf x}d{\bf
x}^{\prime}~{}{\hat{\rho}}({\bf x})u({\bf x}\\!-\\!{\bf
x}^{\prime}){\hat{\rho}}({\bf x}^{\prime}),$ (6)

with the proviso that we take the Coulomb potential with the non-
configurational dielectric constant, or the high frequency dielectric constant
$\varepsilon_{\infty}$ as in the Onsager-Dupuis model, i.e.

$u({\bf x}-{\bf
x}^{\prime})=({1}/{4\pi{\textstyle\varepsilon_{\infty}\varepsilon_{0}})~{}|{\bf
x}-{\bf x}^{\prime}|}^{-1}\,,$ (7)

where $\varepsilon_{\infty}$ is associated with all the relaxation mechanisms
apart from the dipolar one; further, we allow for non-electrostatic
interactions due to the hydration shell of the ions, which generalizes the
role played by the ion-bound Bjerrum defects in the Onsager-Dupuis theory of
ice [18]. This effect corresponds to the coupling between the ion density and
$\hat{\rho}_{(\pm)}({\bf x})$ and $\hat{\cal\bf P}_{j}({\bf x})$. To the
lowest order this coupling can be written via the hydration energy of the form
(for details see Appendix B)

${{\cal H}_{HY}}=\int\\!\\!\int_{V}d{\bf x}d{\bf x}^{\prime}~{}\hat{\rho}({\bf
x})\tilde{u}({\bf x}-{\bf
x}^{\prime})\mbox{\boldmath$\nabla$}\\!\cdot\hat{\cal\bf P}({\bf
x}^{\prime}),$ (8)

where the potential $\tilde{u}({\bf x}-{\bf x}^{\prime})$ is again a short-
range, non-electrostatic potential for the hydration shell interactions. Here
we assume for simplicity that the polarization in the hydration shell of
anions and cations is - apart from its direction - the same.

The total interaction energy functional is thus composed of the standard
electrostatic Coulomb interaction between all the charges in the system, but
in a medium with high frequency dielectric constant, and the non-electrostatic
short-range polarization interactions and ion polarization interactions, and
is written as

${\cal H}={\cal H}_{C}[{\hat{\rho}}({\bf x})]+{{\cal H}_{SR}}[{\hat{\bf
P}({\bf x})]+{{\cal H}_{HY}}[{\hat{\rho}}({\bf x}),{\hat{\bf P}}({\bf x})],}$
(9)

so that the interaction Hamiltonian is ${\cal H}\equiv{\cal
H}[{\hat{\rho}}({\bf x}),{\hat{\bf P}}({\bf x})]$. We note that the
interaction Hamiltonian need not be quadratic and our approach can easily
incorporate higher-order interactions.

##  III Collective description and field theory

From here we derive the partition function in terms of two order parameters
and two auxiliary fields (for details, see Appendix A) as

$\displaystyle{\Xi}$ $\displaystyle\equiv$ $\displaystyle\int{\cal D}[{\cal\bf
P}_{i}({\bf x})]{\cal D}[{\cal E}_{i}({\bf x})]~{}\int{\cal D}[\rho({\bf
x})]{\cal D}[\phi({\bf x})]~{}\times$ (10) $\displaystyle e^{-\beta{\cal
H}[{\cal\bf P}_{i}({\bf x}),\rho({\bf x});{\cal E}_{i}({\bf x}),\phi({\bf
x})]}.$

This expression differs from the standard Coulomb fluid field-theoretical
representation of the partition function [32, 33] by the presence of two
additional independent fields ${\cal\bf P}_{i}({\bf x})$ and ${\cal
E}_{i}({\bf x})$. In what follows we will not analyze the general features of
the partition function, but will limit ourselves immediately to its saddle-
point evaluation, which is taken with respect to the auxiliary fields ${\cal
E}_{i}({\bf x})\longrightarrow i{\cal E}^{*}_{i}({\bf x})$ and $\phi({\bf
x})\longrightarrow i\phi^{*}({\bf x})$. The partition function is then
approximated by $-\ln{{\Xi}}\longrightarrow\beta{\cal H}[{\cal\bf P}_{i}({\bf
x}),\rho({\bf x});i{\cal E}^{*}_{i}({\bf x}),i\phi^{*}({\bf x})]$ where the
effective-field Hamiltonian is obtained from the partial trace over coordinate
and orientation degrees of freedom for the ions and dipoles in the partition
function interacting via the interaction energy functional Eq. (9). It can be
obtained in the form

$\beta{\cal H}[{\cal\bf P}_{i}({\bf x}),\rho({\bf x});i{\cal E}^{*}_{i}({\bf
x}),i\phi^{*}({\bf x})]={\cal H}[{\cal\bf P}_{i}({\bf x}),\rho({\bf
x})]-\int_{V}d{\bf x}~{}{\cal\bf P}_{i}({\bf x}){\cal E}^{*}_{i}({\bf
x})-\int_{V}d{\bf x}~{}\rho({\bf x})\phi^{*}({\bf x})-\int_{V}d^{3}{\bf
x}~{}v^{*}\left({\cal E}^{*}_{i}({\bf x}),\phi^{*}({\bf x})\right),$ (11)

where the first term is the interaction Hamiltonian, Eq. (9).

The expression for the effective-field interaction potential $v^{*}({\cal
E}_{i}({\bf x}),\phi^{*}({\bf x}))=v(i{\cal E}^{*}_{i}({\bf x}),i\phi^{*}({\bf
x}))$ depends crucially on the model used for the dipolar fluid. We
specifically obtain its form for two models: the case of a dipolar model (DM)
is obtained without any constraints on the summation [34], while in the case
of the dipolar Langevin model (DLM) one assumes an underlying lattice of
spacing $a$ (roughly equal to the molecular size) [35], imposing the condition
that each site of the lattice is occupied by only one of the species
(incompressibility condition).

For the DM-model the effective-field interaction potential can be written
explicitly (see Appendix A) as

$\displaystyle v_{DM}^{*}({\cal E}^{*}_{i}({\bf x}),\phi^{*}({\bf
x}))=\lambda~{}\Big{(}\frac{\sinh{u}}{u}+2\tilde{\lambda}_{s}\cosh{w}\Big{)}$
(12)

where $u=\beta p~{}|{\cal E}^{*}({\bf
x})-\mbox{\boldmath$\nabla$}\phi^{*}({\bf x})|$ and $w=\beta e\phi^{*}({\bf
x})$. The DM model implies no constraints on the local density of either ions
or dipoles, and thus exhibits a form of an ideal van’t Hoff gas of ions and
dipoles [36, 34]. For the DLM model we have
$v^{*}_{DLM}=\frac{1}{a^{3}}\ln\Big{(}v^{*}_{DM}/\lambda\Big{)}\,$.

An important new feature of both models is the dependence of the effective
field interaction potential on the field difference ${\cal E}^{*}({\bf
x})-\mbox{\boldmath$\nabla$}\phi^{*}({\bf x})$. The auxiliary field ${\cal
E}^{*}({\bf x})$ here plays the role of a coupling field that quantifies the
degree of non-electrostatic coupling between polarization and electrostatic
fields and is a distinct feature of our approach.

##  IV Mean-field (saddle-point) equations

The saddle-point equations can be obtained for all the variables, i.e. the
order parameters as well as the auxiliary fields. For the two order parameters
we obtain

$\displaystyle{\cal E}^{*}_{i}({\bf x})$ $\displaystyle=$
$\displaystyle-\int_{V}d{\bf x}^{\prime}\tilde{u}_{ij}({\bf x}-{\bf
x}^{\prime}){\cal\bf P}_{j}({\bf x}^{\prime})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+\mbox{\boldmath$\nabla$}_{i}\int_{V}d{\bf
x}^{\prime}~{}\tilde{u}({\bf x}-{\bf x}^{\prime}){\rho}({\bf x}^{\prime}),$
$\displaystyle\phi^{*}({\bf x})$ $\displaystyle=$ $\displaystyle\int_{V}d{\bf
x}^{\prime}u({\bf x}-{\bf x}^{\prime}){\rho}({\bf x}^{\prime})$ (13)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+\int_{V}d{\bf x}^{\prime}\tilde{u}({\bf
x}-{\bf x}^{\prime})\mbox{\boldmath$\nabla$}\cdot{{\cal\bf P}}({\bf
x}^{\prime}).$

Of the two equations the first one contains only non-electrostatic
contributions, while the second one contains combined electrostatic - non-
electrostatic contributions.

After observing that the effective field interaction potential $v^{*}[{\cal
E}^{*}_{i}({\bf x}),\phi^{*}({\bf x})]=v^{*}[u,w]$, the two coupled saddle-
point equations for the auxiliary fields can be written as

$\displaystyle{\cal\bf P}({\bf x})$ $\displaystyle=$ $\displaystyle-(\beta
p)^{2}\left(\frac{1}{u}\frac{\partial v^{*}}{\partial u}\right)\left({\cal
E}^{*}({\bf x})-\mbox{\boldmath$\nabla$}\phi^{*}({\bf x})\right)$ (14)
$\displaystyle{\rho}({\bf x})-\mbox{\boldmath$\nabla$}\cdot{\cal\bf P}({\bf
x})=-(\beta e)\frac{\partial v^{*}}{\partial w}.$

Eqs. (14) together with Eqs. (IV) and (B) constitute the mean-field
description of the structured dielectric.

In order to fully formulate our field equations we need to define the
structural potentials. For this we expand the non-electrostatic potential in
terms of its gradients and consider explicitly the local, second-order and
fourth-order forms

$\displaystyle\tilde{u}_{ij}({\bf x}-{\bf
x}^{\prime})=\tilde{u}_{P}(0)\Big{(}\delta_{ij}\delta({\bf x}-{\bf
x}^{\prime})+\xi^{2}~{}\mbox{\boldmath$\nabla$}^{\prime}_{j}\mbox{\boldmath$\nabla$}_{i}\delta({\bf
x}-{\bf
x}^{\prime})+\zeta^{4}\mbox{\boldmath$\nabla$}^{\prime}_{k}\mbox{\boldmath$\nabla$}^{\prime}_{j}\mbox{\boldmath$\nabla$}_{k}\mbox{\boldmath$\nabla$}_{i}\delta({\bf
x}-{\bf x}^{\prime})\Big{)}+\dots$ (15)

where $\tilde{u}_{P}(0),\xi$ and $\zeta$ are the material constants in the
constitutive relations. The length $\xi$ corresponds to the solvent particle
correlation length and $\zeta$ to the solvent molecular size.

Eq. (14) for ${\bf P}({\bf x})$ can be expressed for the DM and DLM-models as

${\bf P}({\bf x})=-\lambda_{0}(\beta p)^{2}Q(u)\left({\cal E}^{*}({\bf
x})-\mbox{\boldmath$\nabla$}\phi^{*}({\bf x})\right)\,,$ (16)

where $\lambda_{0}\equiv\lambda$ for the dipolar model and $\lambda_{0}\equiv
1/a^{3}$ for the dipolar Langevin model. Writing

${\cal E}^{*}({\bf x})=u_{P}(0){\bf P}({\bf x})-{\cal
F}(\mbox{\boldmath$\nabla$}\cdot{\bf P}({\bf x}))$ (17)

we can express the polarization field as

${\bf P}({\bf x})=\frac{\lambda_{0}(\beta
p)^{2}Q(u)}{1+\lambda_{0}u_{P}(0)(\beta
p)^{2}Q(u)}(\mbox{\boldmath$\nabla$}\phi^{*}({\bf x})+{\cal
F}(\mbox{\boldmath$\nabla$}\cdot{\bf P}))\,.$ (18)

The function $Q(u)$ for the DM-model can be written in the form

$Q(u)={\cal G}(u)\equiv\frac{1}{u}\frac{d}{du}\left(\frac{\sinh(u)}{u}\right)$
(19)

while for the DLM-model it reads as

$Q(u)={\cal H}(u)\equiv\frac{u}{\sinh(u)}{\cal G}(u)\,.$ (20)

The highly nonlinear dependence of the field equations on the polarization
field and the electrostatic field via the function $Q(u)$ is illustrated in
Figure 1 for the DM-model. It shows the function

${\cal G}^{*}(u)=\frac{\lambda(\beta p)^{2}G(u)}{1+\lambda u_{P}(0)(\beta
p)^{2}G(u)}$ (21)

together with the lowest-order nonlinear approximation from its Taylor-
expansion around its minimum.

Figure 1: ${\cal G}^{*}(u)$ in the DM model (blue line) and the lowest-order
nonlinear approximation to ${\cal G}^{*}(u)$ in an expansion around the
minimum (brown line); linear approximation to the central region of the
function (red line).

${\cal G}^{*}(u)$ saturates at both small and large arguments and linearly
interpolates between the two limiting values, see Figure 1.

##  V Linear model equations

In view of the complexity of the fully nonlinear equations, we restrict our
calculations in the following to the linear case, identical for the DM- and
DLM-models, corresponding to the lowest-order term in Eq. (21). Truncating the
non-electrostatic structural interaction potentials Eq. (15) after the second-
order and fourth-order derivative terms, respectively, we refer to the
resulting expressions as Models 1 and 2. As the energy Eq. (8) is already
first order in the derivatives, the simplest possible local from for the ion
density-polarization coupling is

$\tilde{u}({\bf x}-{\bf x}^{\prime})=\alpha~{}\delta({\bf x}-{\bf
x}^{\prime})+\dots$ (22)

where $\alpha$ can be interpreted as the strength of the hydration shell of
the ions. Considering an ion with effective radius $a$, charge $e$ and
polarization of the hydration shell $P_{n}$ with a hydration (free) energy
${\cal F}_{HY}$, we are led to identify $\alpha={\textstyle\frac{1}{3}}{\cal
F}_{HY}a/(P_{n}e)$. Including the ion density-polarization non-electrostatic
coupling given to the lowest order, we term this expression Model 3. We now
formulate the equations of these models explicitly.

Model 1. Stopping at the second-order derivative term in Eq. (15), we have
from Eqs. (IV) for the structural part the relation

${\cal E}^{*}({\bf x})=u_{P}(0)\Big{(}{\bf P}({\bf
x})-\xi^{2}\mbox{\boldmath$\nabla$}\left(\mbox{\boldmath$\nabla$}\cdot{\bf
P}({\bf x})\right)\Big{)}\,.$ (23)

Employing the linearized saddle-point Eq. (14) with

$\displaystyle v^{*}[{\cal E}^{*}({\bf x}),\phi^{*}({\bf
x}),\nabla\phi^{*}({\bf x})]={\textstyle(\beta e)^{2}}\phi^{*}({\bf
x})^{2}+{\textstyle\frac{\lambda(\beta p)^{2}}{6}}u({\bf x})^{2}$

we can eliminate the field ${\cal E}^{*}({\bf x})$ and find

$u_{P}(0)\xi^{2}\mbox{\boldmath$\nabla$}\left(\mbox{\boldmath$\nabla$}\cdot{\bf
P}({\bf x})\right)-\Big{(}u_{P}(0)+\frac{3}{\lambda p^{2}}\Big{)}{\bf P}({\bf
x})=-\mbox{\boldmath$\nabla$}\phi^{*}({\bf x})\,$ (24)
$\displaystyle\mbox{\boldmath$\nabla$}\left(\varepsilon_{\infty}\varepsilon_{0}\mbox{\boldmath$\nabla$}\phi^{*}({\bf
x}){+{\cal\bf P}({\bf x})}\right)={2(\beta e)^{2}\lambda_{s}\phi^{*}({\bf
x})}\,.$ (25)

The coupled equations for $\phi^{*}({\bf x})$ and ${\bf P}({\bf x})$
correspond to the equations derived in Ref. [37] with an opposite sign of the
polarization field.

As in Refs. [14, 37] we introduce the polarization potential
$\phi^{\dagger}({\bf x})$ which can be defined via

${\bf P}({\bf
x})=(\varepsilon-\varepsilon_{\infty})\varepsilon_{0}\mbox{\boldmath$\nabla$}\phi^{\dagger}({\bf
x})\,,$ (26)

making use of the identification

$\varepsilon\varepsilon_{0}\equiv\varepsilon_{\infty}\varepsilon_{0}+\frac{{\textstyle\frac{1}{3}}\lambda
p^{2}}{1+u_{P}(0){\textstyle\frac{1}{3}}\lambda p^{2}}$ (27)

that relates the structural coupling $u_{P}(0)$ to the dielectric constants
and the strength of the water dipole. Integrating the equation for ${\bf P}$,
using the definition Eq. (27) and defining

$\widehat{\xi}^{2}\equiv(\varepsilon-\varepsilon_{\infty})\varepsilon_{0}u_{P}(0)\,\xi^{2},$
(28)

the mean-field equations for our Model 1 read as

$\displaystyle\mbox{\boldmath$\nabla$}^{2}\phi^{\dagger}({\bf x})$
$\displaystyle=$ $\displaystyle{\widehat{\xi}^{-2}}\left(\phi^{\dagger}({\bf
x})-\phi^{*}({\bf x})\right),$ (29)
$\displaystyle\mbox{\boldmath$\nabla$}^{2}\phi^{*}({\bf x})$ $\displaystyle=$
$\displaystyle~{}{\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}\kappa_{D}^{2}\phi^{*}({\bf
x})$
$\displaystyle+\,\,{\widehat{\xi}^{-2}}\left({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1\right)\left(\phi^{*}({\bf
x})-\phi^{\dagger}({\bf x})\right)\,,$

in which the inverse square of Debye length is defined by
$\kappa_{D}^{2}\equiv{2(\beta
e)^{2}\lambda_{s}}/{\varepsilon\varepsilon_{0}}$. Except for the sign of
$\phi^{\dagger}$ these equations are the same as the Onsager-Dupuis equations
[12].

Model 2. Keeping the terms up to the fourth order in derivatives in Eq. (15)
we find, introducing a second redefined length $\widehat{\zeta}_{0}$

$\widehat{\zeta_{0}}^{4}\equiv(\varepsilon-\varepsilon_{\infty})\varepsilon_{0}u_{P}(0)\,\zeta^{4}$
(30)

the set of equations

$\displaystyle\widehat{\zeta_{0}}^{4}\mbox{\boldmath$\nabla$}^{2}\mbox{\boldmath$\nabla$}^{2}\phi^{\dagger}({\bf
x})-{\widehat{\xi}^{2}}\mbox{\boldmath$\nabla$}^{2}\phi^{\dagger}({\bf x})$
$\displaystyle=$ $\displaystyle\phi^{*}({\bf x})-\phi^{\dagger}({\bf x}),$
$\displaystyle\mbox{\boldmath$\nabla$}^{2}\phi^{*}({\bf
x})+\left({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1\right)\mbox{\boldmath$\nabla$}^{2}\phi^{\dagger}({\bf
x})$ $\displaystyle=$
$\displaystyle~{}{\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}\kappa_{D}^{2}\phi^{*}({\bf
x}),$ (31)

which depend on three characteristic lengths, $1/\kappa_{D},\xi$ and
$\zeta_{0}$.

Model 3. The ion density-polarization coupling is incorporated via the lowest
order expansion Eq. (22) and gives rise to the following equations for the
electrostatic and polarization potentials (for details on the derivation, see
Appendix B)

$\displaystyle\widehat{\zeta}^{4}\mbox{\boldmath$\nabla$}^{4}\phi^{\dagger}({\bf
x})-\widehat{\xi}^{2}\mbox{\boldmath$\nabla$}^{2}\phi^{\dagger}({\bf
x})+\widehat{\alpha}^{2}\mbox{\boldmath$\nabla$}^{2}\phi^{*}({\bf
x})=\phi^{*}({\bf x})-\phi^{\dagger}({\bf x}),$
$\displaystyle\left(\textstyle{\frac{\varepsilon}{\varepsilon_{\infty}}}-1\right)[1-\widehat{\alpha}^{2}\mbox{\boldmath$\nabla$}^{2}]\mbox{\boldmath$\nabla$}^{2}\phi^{\dagger}+\mbox{\boldmath$\nabla$}^{2}\phi^{*}=\textstyle{\frac{\varepsilon}{\varepsilon_{\infty}}}\kappa_{D}^{2}\phi^{*}$
(32)

with the additional length $\widehat{\zeta}$

$\widehat{\zeta}^{4}\equiv\left(u_{P}(0)\zeta^{4}-\alpha^{2}\varepsilon_{\infty}\varepsilon_{0}\right)(\varepsilon-\varepsilon_{\infty})\varepsilon_{0}=\widehat{\zeta}^{4}_{0}-\widehat{\alpha}^{4}\left({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1\right)\\\
$ (33)

where the definition
$\widehat{\alpha}^{2}\equiv\alpha\varepsilon_{0}\varepsilon_{\infty}$ has been
introduced. This model obviously depends on all four characteristic lengths
that we consider in this model catalogue.

##  VI Decoupling and solving the equations

Since we will consider one-dimensional geometries like a slit geometry $-L\leq
z\leq L$ or a single-plate geometry $0\leq z<\infty$ we assume ${\bf
x}\equiv(0,0,z)$ and the mean-field equations reduce to ordinary differential
equations. The first step in solving these equations is then the decoupling of
the electrostatic and polarization potentials, $\phi^{*}$ and
$\phi^{\dagger}$, respectively.

###  VI.1 Model 1

In the case of Model 1 the equations can be written in matrix form
$\frac{d^{2}}{dz^{2}}\Phi_{1}(z)={\cal M}_{1}\Phi_{1}(z)$ with the
introduction of the composite field
$\Phi_{1}(z)\equiv(\phi^{*}(z),\phi^{\dagger}(z))$. The coupling matrix ${\cal
M}_{1}$ is given by

$\displaystyle{\cal
M}_{1}=\left(\begin{array}[]{cc}\frac{\varepsilon}{\varepsilon_{\infty}}\kappa_{D}^{2}+\left(\frac{\varepsilon}{\varepsilon_{\infty}}-1\right)\widehat{\xi}^{-2}&-\left(\frac{\varepsilon}{\varepsilon_{\infty}}-1\right)\widehat{\xi}^{-2}\\\
-\widehat{\xi}^{-2}&\widehat{\xi}^{-2}\end{array}\right)\,.$ (36)

Diagonalizing the matrix and rescaling the eigenvalues with
$\lambda\longrightarrow\lambda/\widehat{\xi}^{2}$ leads to a quadratic
eigenvalue equation given by

$\lambda^{2}-\lambda~{}{\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}\left(1+(\kappa_{D}\widehat{\xi})^{2}\right)+{\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}(\kappa_{D}\widehat{\xi})^{2}=0,$
(37)

where $\kappa^{\pm}=\sqrt{\lambda_{\pm}}$ are then the two inverse decay
lengths corresponding to the two eigenvalues, which are both real and
positive. They correspond to the decay lengths of the electrostatic and the
polarization potentials, respectively, and coincide with the result of the
linearized Onsager-Dupuis theory [15, 37]. The decay lengths are shown in
Figure 2 as functions of $\kappa_{D}\hat{\xi}$ in accord with earlier results,
see [38].

Figure 2: The eigenvalues of the matrix ${\cal M}_{1}$, with $\lambda_{1}\leq
1$, as a function of dimensionless parameter $\kappa_{D}\hat{\xi}$. The
horizontal line corresponds to $\lambda=1$. Insert: the ratio of $\lambda_{1}$
and $\lambda_{2}$.

The key feature of our description, being the decoupling of the degrees of
freedom of the electrostatic field and the polarization field, translates into
the requirement of separate boundary conditions for both. In order to solve
the linear model equations, we therefore need to specify the boundary
conditions for $-\mbox{\boldmath$\nabla$}\phi^{*}\cdot{\bf n}$ and ${\bf
P}\cdot{\bf n}$ at the bounding surface(s) with normal vector ${\bf n}$.
Alternatively we can assume that the bounding surfaces carry ‘polarization
charges’ $\sigma_{P}$

$\displaystyle\mbox{\boldmath$\nabla$}\cdot{\bf P}({\bf x})=\rho_{P}({\bf
x})=\sigma_{P}\delta(z-z_{0}),$ (38)

as well as the standard electric charges $\sigma$

$\displaystyle\varepsilon_{\infty}\varepsilon_{0}\mbox{\boldmath$\nabla$}\cdot{\bf
E}({\bf x})=\rho({\bf x})=\sigma\delta(z-z_{0}),$ (39)

where $z$ is the coordinate perpendicular to the boundary. Figure 3 shows a
solution of Model 1 for the case of prescribed electric field and polarization
at the boundary of a slit of width $2L$, i.e., $z_{0}=\pm L$. The solutions
can be compared with the results derived in [38].

Figure 3: Exemplary solutions of Model 1 for $\phi^{*}(z)$ (blue line) and
$\phi^{\dagger}(z)$ (yellow line). The boundary conditions are:
$(\phi^{*})^{\prime}(-L)=-E_{0}$, $(\phi^{*})^{\prime}(L)=E_{0}$,
$(\phi^{\dagger})^{\prime}(-L)=(\phi^{\dagger})^{\prime}(L)=0$.

###  VI.2 Model 2

For Model 2 (and later Model 3), we can proceed in a similar manner. In this
case we first have to reduce the order of the equations from four to two by
introducing an auxiliary potential $\psi$ which in the one-dimensional
geometry leads to

$\frac{d^{2}\phi^{\dagger}(z)}{dz^{2}}=\frac{1}{\widehat{\xi}^{2}}\psi(z)\,.$
(40)

In terms of the three-potential composite field
$\Phi_{2}(z)\equiv(\phi^{*}(z),\phi^{\dagger}(z),\psi(z))$ we now have the
expression $\frac{d^{2}}{dz^{2}}\Phi_{2}(z)={\cal M}_{2}\Phi_{2}(z)$, with the
coupling matrix

$\displaystyle{\cal
M}_{2}=\left(\begin{array}[]{ccc}\frac{\varepsilon}{\varepsilon_{\infty}}\kappa_{D}^{2}&0&-\left(\frac{\varepsilon}{\varepsilon_{\infty}}-1\right)\widehat{\xi}^{-2}\\\
0&0&\widehat{\xi}^{-2}\\\
\frac{\widehat{\xi}^{2}}{\widehat{\zeta_{0}^{4}}}&-\frac{\widehat{\xi}^{2}}{\widehat{\zeta_{0}}^{4}}&\frac{\widehat{\xi}^{2}}{\widehat{\zeta_{0}}^{4}}\end{array}\right)\,.$
(44)

The corresponding eigenvalue equations are thus of third order, depending on
four parameters. One can scale out one of them, remaining with the ratios.
E.g., with the substitution
$\lambda\longrightarrow\lambda/\widehat{\zeta_{0}}^{2}$ the characteristic
polynomial is given by an expression containing only the ratios of the model
parameters

$\displaystyle\lambda^{3}$ $\displaystyle-$
$\displaystyle\lambda^{2}\left[\left({\textstyle\frac{\widehat{\xi}}{\widehat{\zeta_{0}}}}\right)^{2}+{\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}(\kappa_{D}\widehat{\zeta_{0}})^{2}\right]+$
(45)
$\displaystyle+~{}\lambda{\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}\left(1+(\kappa_{D}\widehat{\xi})^{2}\right)-{\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}(\kappa_{D}\widehat{\zeta_{0}})^{2}=0\,.$

In order to see the character of the roots we need to look at the discriminant
$\Delta$ of the cubic polynomial. Since both $(\kappa_{D}\widehat{\xi})^{2}\ll
1$ and $(\kappa_{D}\widehat{\zeta_{0}})^{2}\ll 1$, the discriminant is
approximately given by the expression

$\displaystyle\Delta=\left({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}\right)^{2}\left(\left({\textstyle\frac{\widehat{\xi}}{\widehat{\zeta_{0}}}}\right)^{4}-4{\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}\right)$
(46)
$\displaystyle-4\left({\textstyle\frac{\widehat{\xi}}{\widehat{\zeta_{0}}}}\right)^{2}\left(\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}\right)\left(\left({\textstyle\frac{\widehat{\xi}}{\widehat{\zeta_{0}}}}\right)^{4}-4.5{\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}\right)(\kappa_{D}\widehat{\zeta_{0}})^{2}$
$\displaystyle-27\left({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}\right)^{2}(\kappa_{D}\widehat{\zeta_{0}})^{4}\,.$

In the limit $\kappa_{D}\widehat{\zeta_{0}}\rightarrow 0$ one sees that the
sign of the discriminant depends on the difference
$\left({\textstyle\frac{\widehat{\xi}}{\widehat{\zeta_{0}}}}\right)^{4}-4{\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}$
i.e. in a highly nonlinear way on the length-ratio $\xi/\zeta_{0}$ and the
polarity of the solvent.

Figure 4 shows the discriminant, Eq. (49), as a function of $\xi/\zeta_{0}$
and $\varepsilon/\varepsilon_{\infty}$. The function is cut by the zero level
such that for positive $\Delta$ there are three real roots, while for negative
$\Delta$ there is one real and two complex conjugate roots. The discriminant
of Model 2 therefore signals a change in the character of the solutions from
exponential decay to damped oscillatory behaviour.

Figure 4: The discriminant $\Delta$ for Model 2 as a function of
$\xi/\zeta_{0}$ and $\varepsilon/\varepsilon_{\infty}$ in relevant parameter
ranges; $\kappa_{D}\widehat{\zeta_{0}}=10^{-4}$. For $\Delta>0$ there are
three real roots of the eigenvalue equation, while for $\Delta<0$, one real
and two complex conjugate roots arise.

Figure 5 displays the numerical solution of the equations of Model 2 for the
case the plate geometry. The boundary conditions for the fields $\phi^{*}$ and
$\phi^{\dagger}$ have been chosen in an antagonistic fashion, which show that
both fields behave as separate entities near the wall. For distances $z\gg 0$,
both fields merge into each other and fulfill the definition of the
polarization field of Eq. (26) in the bulk.

Figure 5: Solutions to the equations of Model 2. Blue: $\phi^{*}(z)$, Brown:
$\phi^{\dagger}(z)$, Green: auxiliary field $\psi(z)$. For $z\gg 0$, the
solutions of $\phi^{*}$ and $\phi^{\dagger}$ merge, establishing the validity
of the standard relation between the polarization and electric fields in the
bulk.

Neglecting ionic screening altogether, i.e. $\kappa_{D}=0$, the equations for
Model 2 can be simplified further and reduced to a single fourth-order
equation for the auxiliary field $\psi(z)$ in the form

$\widehat{\zeta_{0}}^{4}\psi^{iv}(z)-{\widehat{\xi}^{2}}\psi^{\prime\prime}(z)+({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}})\psi(z)=0\,.$
(47)

Apart from the definition of the constants, the basic equation then reduce to
the form introduced in Ref. [11]. We note that fourth-order equations can be
obtained also by considering quadrupole coupling terms pertinent to the
quadrupolarizable solvent models [39, 40].

Making the substitution $\kappa\longrightarrow\kappa/\hat{\zeta}_{0}$, so that
the decay length is in units of $\hat{\zeta}$, the characteristic polynomial
for the inverse decay length, Eq.(47), can be obtained in the simple form

$\displaystyle\kappa^{4}-\kappa^{2}~{}\left({\textstyle\frac{\widehat{\xi}}{\widehat{\zeta_{0}}}}\right)^{2}+{\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}=0.$
(48)

The solution for the eigenvalues are then of two types: real and complex,
depending on the sign of the discriminant

$\Delta=\left({\textstyle\frac{\widehat{\xi}}{\widehat{\zeta_{0}}}}\right)^{4}-4{\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}.$
(49)

which amounts to the identical dependence found in the discriminant $\Delta$
of the full system of equations. In the case of
$\widehat{\xi}^{2}>\sqrt{{4\varepsilon}/{\varepsilon_{\infty}}}~{}\widehat{\zeta_{0}}^{2}$,
the eigenvalues are real and obtained as

$\displaystyle\kappa_{1,2,3,4}=\pm\frac{1}{\sqrt{2}}\left(\left({\textstyle\frac{\widehat{\xi}}{\widehat{\zeta_{0}}}}\right)^{2}\pm{\delta}^{1/2}\right)^{1/2},$
(50)

while for
$\widehat{\xi}^{2}\leq\sqrt{{4\varepsilon}/{\varepsilon_{\infty}}}~{}\widehat{\zeta_{0}}^{2}$,
the eigenvalues are complex conjugate pairs and given by

$\displaystyle\kappa_{1,2,3,4}$ $\displaystyle=$
$\displaystyle\pm\left(\kappa_{d}\pm i\kappa_{o}\right)$ (51)

where

$\displaystyle\kappa_{d,o}={\textstyle\frac{1}{2}\left({\left({\left({\textstyle\frac{\widehat{\xi}}{\widehat{\zeta_{0}}}}\right)^{4}+(-\delta)}\right)^{1/2}\pm\left({\textstyle\frac{\widehat{\xi}}{\widehat{\zeta_{0}}}}\right)^{2}}\right)^{1/2}}.$
(52)

This last case is also the one considered in Ref. [11]. Clearly, even with the
fourth-order derivative polarization functional the corresponding
characteristic lengths can be real, so that periodic solutions are not
universal but depend on the sign of the discriminant Eq. (49). Since
$\zeta_{0}$ is comparable to the size of the molecule and $\xi$ is the
macroscopic correlation length, the periodic solution are expected for
sufficiently polar media $4{\varepsilon}/{\varepsilon_{\infty}}\gg 1$.

###  VI.3 Model 3

A similar, albeit approximate solution can be found for the case of Model 3
when one neglects the solvent molecular length $\zeta$ compared to the ion
density-polarization coupling $\alpha$. We then have the coupling matrix for
the composite potential
$\Phi_{3}(z)\equiv(\phi^{*}(z),\phi^{\dagger}(z),\psi(z))$ by the expression

$\displaystyle{\cal
M}_{3}=\left(\begin{array}[]{ccc}\widehat{\alpha}^{-2}&-\widehat{\alpha}^{-2}&\widehat{\alpha}^{-2}\\\
0&0&\widehat{\xi}^{-2}\\\
\frac{\left(\frac{\widehat{\xi}}{\widehat{\alpha}}\right)^{2}-\frac{\varepsilon}{\varepsilon_{\infty}}(\kappa_{D}\xi)^{2}}{({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1)\widehat{\alpha}^{2}}&-\frac{\left(\frac{\widehat{\xi}}{\widehat{\alpha}}\right)^{2}}{({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1)\widehat{\alpha}^{2}}&\frac{\left(\frac{\widehat{\xi}}{\widehat{\alpha}}\right)^{2}+({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1)}{({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1)\widehat{\alpha}^{2}}\end{array}\right)\,.$
(56) (57)

The corresponding characteristic polynomial is again of third order and
contains four model parameters. Choosing one of them, e.g. $\hat{\alpha}$, and
rescaling the roots by $\lambda\longrightarrow\lambda/\hat{\alpha}^{2}$ the
characteristic polynomial comes out as

$\displaystyle\lambda^{3}-$
$\displaystyle\lambda^{2}\left[\frac{\left(\frac{\widehat{\xi}}{\widehat{\alpha}}\right)^{2}+2\left(\frac{\varepsilon}{\varepsilon_{\infty}}-1\right)}{({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1)}\right]+$
(58)
$\displaystyle+~{}\lambda\frac{\frac{\varepsilon}{\varepsilon_{\infty}}\left(1+(\kappa_{D}\widehat{\xi})^{2}\right)}{({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1)}-\frac{\frac{\varepsilon}{\varepsilon_{\infty}}(\kappa_{D}\widehat{\alpha})^{2}}{({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1)}=0\,.$

Comparing this expression with Eq. (45) we see that apart from the denominator
$({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1)$ there is an almost
complete correspondence between the molecular length $\zeta$ and the secondary
hydration structural length $\alpha$. The discussion of the discriminant
$\Delta$ of the cubic polynomial for Model 3 therefore shows the same
characteristic behaviour as for Model 2. The criterion for the change of the
eigenvalue spectrum from real to real and complex conjugates is for this model
replaced by the expression

${\textstyle\left(\frac{\widehat{\xi}}{\widehat{\alpha}}\right)^{2}}<2\left[1-\textstyle{\frac{\varepsilon}{\varepsilon_{\infty}}}\left(1-{\textstyle(1-\frac{\varepsilon_{\infty}}{\varepsilon}})^{1/2}\right)\right],$
(59)

which in the limit $\varepsilon\gg\varepsilon_{\infty}$ reduces to
$\xi^{2}/\widehat{\alpha}^{2}<1$. Since by definition $\widehat{\alpha}^{2}$
can be negative, in this case the eigenvalues are always real and complex
conjugates; in the case $\widehat{\alpha}^{2}>0$, if the correlation length is
larger than the ion-density polarization coupling, all eigenvalues are real.
The plot of the discriminant for the case $\widehat{\alpha}^{2}>0$ is shown in
Figure 6.

Figure 6: The discriminant $\Delta$ for Model 3 as a function of
$\widehat{\xi}/\widehat{\alpha}$ and $\varepsilon/\varepsilon_{\infty}$ in
relevant parameter ranges; $\kappa_{D}\widehat{\zeta}=10^{-4}$. For $\Delta>0$
there are three real roots of the eigenvalue equation, while for $\Delta<0$,
one real and two complex conjugate roots arise.

Profiting further from our insights from Model 2, we also consider the case of
a vanishing Debye-screening length in Model 3. We then find the expression

$\widehat{\zeta_{0}}^{4}(\phi^{\dagger})^{iv}-\left(2\widehat{\alpha}^{2}({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1)+\widehat{\xi}^{2}\right)(\phi^{\dagger})^{{}^{\prime\prime}}+{\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}\phi^{\dagger}=0\,.$
(60)

Clearly the hydration coupling corresponds to secondary hydration effects with
the square of the effective structural correlation length given by
$2\widehat{\alpha}^{2}({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1)+\widehat{\xi}^{2}$.
Making the substitution
$\kappa\longrightarrow\kappa/\left(({\textstyle\frac{\varepsilon_{\infty}}{\varepsilon}})^{1/4}\widehat{\zeta_{0}}\right)$,
the characteristic polynomial of Eq.(58) can be obtained in the simple form

$\displaystyle\kappa^{4}-\kappa^{2}{\left({\textstyle\frac{\varepsilon_{\infty}}{\varepsilon}}\right)^{1/2}\left(2\left({\textstyle\frac{\widehat{\alpha}}{\zeta_{0}}}\right)^{2}({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1)+\left({\textstyle\frac{\widehat{\xi}}{\widehat{\zeta_{0}}}}\right)^{2}\right)}+1=0.$

Equivalently, one can choose a different rescaling
$\kappa\longrightarrow\kappa/\left(\left({\textstyle\frac{\varepsilon_{\infty}}{\varepsilon}}\right)\left(2\widehat{\alpha}^{2}({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1)+\widehat{\xi}^{2}\right)\right)^{1/2}$
which then leads to the characteristic equation

$\left({\textstyle\frac{\varepsilon_{\infty}}{\varepsilon}}\right)^{3}\widehat{\zeta_{0}}^{4}\left(2\widehat{\alpha}^{2}({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1)+\widehat{\xi}^{2}\right)^{2}\kappa^{4}-\kappa^{2}+1=0\,.$
(62)

The case
$\left({\textstyle\frac{\varepsilon_{\infty}}{\varepsilon}}\right)^{3}\widehat{\zeta_{0}}^{4}\left(2\widehat{\alpha}^{2}({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1)+\widehat{\xi}^{2}\right)^{2}>1/2$
has been discussed in [41] in the context of ionic liquids but is formally
equivalent to the above case. The exact profile for the case of a planar
charged surface $0\leq z<\infty$ is given in Appendix D. The corresponding
profiles for $\phi^{*}(z)$ and $\phi^{\dagger}(z)$ are shown in Figure 7.

Figure 7: Solutions of the equations of the reduced Model 3, Eq. (60);
$\phi^{*}(z)$ (blue line); $\phi^{\dagger}(z)$ (brown line); see Appendix D
for the analytic expressions for $\phi^{*}(z)$ and $\phi^{\dagger}(z)$.

##  VII Discussion and Conclusions

In this paper we have derived a field-theoretic description of structured
liquid dielectrics, taking the Onsager-Dupuis theory of ice as a motivation.
Starting from the introduction of two separate order parameters for the
polarization field of the solvent and the charge density of the solvated ions,
we provide a general formalism to introduce a family of models of different
degree of coarse-graining. Considering a harmonic approximation for the free
energy of the polarization field and a dipolar model for the solvent, we
derive novel nonlinear equations that emerge from the present theory.

As a first application we discuss three models defined in the linear limit of
the theory, which lead to coupled differential equation systems for the
electrostatic and polarization potentials. The generic behaviour of the
solutions is determined from a decoupling of these equations. The resulting
discriminants of the eigenvalue equations define the solution behaviour as
functions of model structural lengths and solvent polarity. Since the model
equations allow for independent boundary conditions for the electrostatic and
polarization potentials they have a high flexibility in the discussion of
physical systems, which we illustrate for plate and slab geometries. In the
case of the ion charge density-solvent polarization coupling our approach
allowed us to describe the case of secondary hydration effects where the
structural solvent correlation length has solvent as well as ion
contributions.

We believe that, beyond the introduction of our approach and first results in
the linear regime, our theory is capable of addressing a wide class of
continuum models for polarization phenomena in confined solvents, in
particular phenomena pertinent to the nano-confined aqueous systems [25, 26].
Further extension of our work to allow for nonlinear polarization functionals
in the bulk and the incorporation of surface polarization functionals are also
immediately feasible.

###### Acknowledgements.

R.P. acknowledges the support of the University of Chinese Academy of Sciences
and the funding from the NSFC under grant No. 12034019.

##  VIII Appendices

##  Appendix A Collective description and field theory

To express the partition function in the form of a field theory, we take the
same approach as with anisotropic Coulomb fluids [31], the outlines of which
we briefly reiterate below.

We write the total energy with two collective variables, the polarization and
the charge density, Eqs. (1) and (4), by introducing the decomposition of
unity

$\displaystyle 1\equiv\int{\cal D}[{\cal\bf P}_{i}({\bf x})]~{}\Pi_{{\bf
x}}\delta({\cal\bf P}_{i}({\bf x})-\hat{\cal\bf P}_{i}({\bf x}))\times$
$\displaystyle\int{\cal D}[\rho({\bf x})]~{}\Pi_{{\bf x}}\delta(\rho({\bf
x})-\hat{\rho}({\bf x}))$ (63)

together with the functional integral representation of the delta functions in
terms of their respective auxiliary fields ${\cal E}_{i}(x)$, $\phi(x)$:

$\displaystyle\Pi_{{\bf x}}\delta({\cal\bf P}_{i}({\bf x})-\hat{\cal\bf
P}_{i}({\bf x}))\longrightarrow\int{\cal D}[{\cal E}_{i}({\bf
x})]~{}e^{-i{\cal E}_{i}({\bf x})({\cal\bf P}_{i}({\bf x})-\hat{\cal\bf
P}_{i}({\bf x}))}$ $\displaystyle\Pi_{{\bf x}}\delta(\rho({\bf
x})-\hat{\rho}({\bf x}))\longrightarrow\int{\cal D}[\phi({\bf
x})]~{}e^{-i\phi({\bf x})(\rho({\bf x})-\hat{\rho}({\bf x}))}.$ (64)

Notably, we introduce the microscopic polarization as a separate order
parameter, independently of the ion charge density. This will eventually lead
to a two-order parameter representation of the partition function.

Using these relations in the definition of the grand canonical partition
function

$\displaystyle{\Xi}$ $\displaystyle\equiv$
$\displaystyle\sum_{N}\frac{\lambda^{N}}{N!}\sum_{N_{+}}\sum_{N_{-}}\frac{\lambda_{(+)}^{N_{+}}\lambda_{(-)}^{N_{-}}}{N_{+}!N_{-}!}$
(65) $\displaystyle\int{\cal D}[{\bf x}_{N}]{\cal D}[{\bf n}_{N}]\int{\cal
D}[{\bf x}_{N_{\pm}}]~{}e^{-\beta{\cal H}[{\bf x}_{N},{\bf x}_{N_{\pm}}]},$

with $\lambda$ and $\lambda_{(\pm)}$ as the absolute activities of the dipolar
and ionic components, the introduction of collective variables ${\cal\bf
P}_{i}({\bf x})$ and $\rho({\bf x})$ leads to the following field-
representation of the partition function

$\displaystyle{\Xi}$ $\displaystyle=$ $\displaystyle\int{\cal D}[{\cal\bf
P}_{i}({\bf x})]{\cal D}[{\cal E}_{i}({\bf x})]\int{\cal D}[\rho({\bf
x})]{\cal D}[\phi({\bf x})]~{}\times$ (66) $\displaystyle e^{-\beta\tilde{\cal
H}[{\cal\bf P}_{i}({\bf x}),~{}\rho({\bf x});{\cal E}_{i}({\bf x}),\phi({\bf
x})]-V[{\cal E}_{i}({\bf x}),\phi({\bf x})]}.$

The field action above decouples into two terms, of which the first one
$\beta\tilde{\cal H}[{\cal\bf P}_{i}({\bf x}),~{}\rho({\bf x});{\cal
E}_{i}({\bf x}),\phi({\bf x})]$ depends on the collective order parameters
${\cal\bf P}_{i}({\bf x})$ and $\rho({\bf x})$ as well as the auxiliary fields

$\displaystyle\beta\tilde{\cal H}[{\cal\bf P}_{i}({\bf x}),~{}\rho({\bf
x});{\cal E}_{i}({\bf x}),\phi({\bf x})]=$
$\displaystyle~{}~{}~{}~{}~{}{\textstyle\frac{1}{2}}\int\\!\\!\int_{V}d{\bf
x}d{\bf x}^{\prime}~{}{\cal\bf P}_{i}({\bf x})\tilde{u}_{ij}({\bf x}-{\bf
x}^{\prime}){\cal\bf P}_{j}({\bf x}^{\prime})+$
$\displaystyle~{}~{}~{}~{}~{}{\textstyle\frac{1}{2}}\int\\!\\!\int_{V}d{\bf
x}d{\bf x}^{\prime}~{}{\rho}({\bf x})u({\bf x}-{\bf x}^{\prime}){\rho}({\bf
x}^{\prime})+$ $\displaystyle~{}~{}~{}~{}~{}i\int_{V}d{\bf x}~{}{\cal\bf
P}_{i}({\bf x}){\cal E}_{i}({\bf x})+i\int_{V}d{\bf x}~{}\rho({\bf
x})\phi({\bf x}),$ (67)

and is a sum of the interaction Hamiltonian for non-electrostatic polarization
interactions and Coulomb interactions between mobile and bound charges.

The second component, the effective field interaction potential, corresponding
to the integration over the configurational variables ${\bf x}_{N}$, ${\bf
n}_{N}$ and ${\bf x}_{N_{\pm}}$ of the solvent dipoles and electrolyte ions,
depends only on the auxiliary fields ${\cal E}_{i}({\bf x})$ and $\phi({\bf
x})$ as

$\displaystyle V[{\cal E}_{i}({\bf x}),\phi({\bf x})]$ $\displaystyle\equiv$
$\displaystyle\log{\Big{(}\sum_{N}\
\sum_{N^{+}}{\sum_{N^{-}}}\frac{\lambda^{N}}{N!}\frac{\lambda_{(+)}^{N^{+}}\lambda_{(-)}^{N^{-}}}{N^{+}!N^{-}!}\int{\cal
D}[{\bf x}_{N}]{\cal D}[{\bf n}_{N}]~{}\int{\cal D}[{\bf
x}_{N_{\pm}}]~{}e^{-\beta\tilde{\cal H}^{*}[{\cal E}_{i}({\bf
x}_{i}),\phi({\bf x}_{i})]}\Big{)}}$ $\displaystyle=$
$\displaystyle\lambda\int_{V}d^{3}{\bf x}\int_{\Omega}d{\bf
n}~{}e^{-ip~{}\\!{\bf n}_{i}~{}\\!\big{(}{\cal E}_{i}({\bf
x})-\mbox{\boldmath$\nabla$}_{i}\phi({\bf
x})\big{)}}+2\lambda_{s}\int_{V}d^{3}{\bf x}\cos{\beta e\phi({\bf x})}$
$\displaystyle=$ $\displaystyle\int_{V}d^{3}{\bf
x}~{}\Big{(}\lambda\int_{\Omega}\sin{\theta}d\theta~{}e^{-ip\cos{\theta}~{}|{\cal
E}_{i}({\bf x})-\mbox{\boldmath$\nabla$}_{i}\phi({\bf
x})|}+2\lambda_{s}\cos{\beta e\phi({\bf x})}\Big{)}$ $\displaystyle=$
$\displaystyle\lambda\int_{V}d^{3}{\bf x}~{}\Big{(}\frac{\sin{p|{\cal
E}_{i}({\bf x})-\mbox{\boldmath$\nabla$}_{i}\phi({\bf x})|}}{p|{\cal
E}_{i}({\bf x})-\mbox{\boldmath$\nabla$}_{i}\phi({\bf
x})|}+2\tilde{\lambda}_{s}\cos{\beta e\phi({\bf x})}\Big{)}=\int_{V}d^{3}{\bf
r}~{}v\left({\cal E}_{i}({\bf x}),\phi({\bf x})\right),$ (68)

where the single-particle configurational Hamiltonian is

$\displaystyle\beta\tilde{\cal H}^{*}[{\cal E}_{i}({\bf x}_{i}),\phi({\bf
x}_{i})]$ $\displaystyle=$ $\displaystyle-i~{}p\sum_{N}{\bf
n}_{i}~{}\\!\big{(}{\cal E}_{i}({\bf
x}_{i})-\mbox{\boldmath$\nabla$}_{i}\phi({\bf x}_{i})\big{)}-$ (69)
$\displaystyle-i\sum_{N_{+}}~{}\phi({\bf x}_{i})+i\sum_{N_{-}}~{}\phi({\bf
x}_{i}).~{}$

We have rescaled all the interaction energies in terms of the thermal energy
and assumed that the salt ion activities are the same,
$\lambda_{(+)}=\lambda_{(-)}=\lambda_{s}$ [36, 34]. Note that above we also
redefined $\lambda_{s}/\lambda\longrightarrow\tilde{\lambda}_{s}$.

Eq. (68) describes an ideal van’t Hoff gas of ions and dipoles, corresponding
to the dipolar model (DM). This expression differs from the standard Coulomb
fluid field-theoretical representation of the partition function [32, 33] by
the presence of two additional independent fields ${\cal\bf P}_{i}({\bf x})$
and ${\cal E}_{i}({\bf x})$. The saddle-point evaluation of the partition
function is taken with respect to the auxiliary fields ${\cal E}_{i}({\bf
x})\longrightarrow i{\cal E}^{*}_{i}({\bf x})$ and $\phi({\bf
x})\longrightarrow i\phi^{*}({\bf x})$. The partition function is then
approximated by

$\displaystyle-\ln{{\Xi}}\longrightarrow\beta{\cal H}[{\cal\bf P}_{i}({\bf
x}),\rho({\bf x});i{\cal E}^{*}_{i}({\bf x}),i\phi^{*}({\bf x})],$ (70)

where the effective-field Hamiltonian is obtained from the partial trace over
coordinate and orientation degrees of freedom for the ions and dipoles in the
partition function interacting via the interaction energy functional Eq. (9).

The expression for the effective field interaction potential $V^{*}[{\cal
E}_{i}({\bf x}),\phi^{*}({\bf x})]=V[i{\cal E}^{*}_{i}({\bf x}),i\phi^{*}({\bf
x})]$ depends crucially on the model used for the dipolar fluid, see main
text. The saddle-point form of the DM effective field interaction potential
$V^{*}[i{\cal E}^{*}_{i}({\bf x}),i\phi^{*}({\bf x})]$ can be finally written
explicitly in the form of Eq. (68) but with $\sin{\beta p|{\cal E}^{*}({\bf
x})-\mbox{\boldmath$\nabla$}\phi^{*}({\bf x})|}\longrightarrow\sinh{\beta
p|{\cal E}^{*}({\bf x})-\mbox{\boldmath$\nabla$}\phi^{*}({\bf x})|}$ and
$\cos{\beta e\phi^{*}({\bf x})}\longrightarrow\cosh{\beta e\phi^{*}({\bf
x})}$.

##  Appendix B Hydration shell coupling and the derivation of Model 3

Since it couples the ion density and polarization, the microscopic hydration
energy of Model 3 is actually of the form

$\displaystyle{{\cal H}_{HY}}$ $\displaystyle=$
$\displaystyle{\textstyle\frac{1}{2}}\int\\!\\!\int_{V}d{\bf x}d{\bf
x}^{\prime}~{}\hat{\rho}_{(\pm)}({\bf x})\tilde{u}({\bf x}\\!-\\!{\bf
x}^{\prime})\mbox{\boldmath$\nabla$}\\!\cdot\\!\hat{\cal\bf P}({\bf
x}^{\prime}),$ (71)

where the potential $\tilde{u}({\bf x}-{\bf x}^{\prime})$ is a short-range,
non-electrostatic potential describing the hydration shell interactions. We
assumed that the hydration polarization for anions and cations is - apart from
the direction - the same. This could be of course elaborated further. Notice,
however, that the hydration energy can be recast as

$\displaystyle{{\cal H}_{HY}}$ $\displaystyle=$
$\displaystyle{\textstyle\frac{1}{2}}\\!\\!\int\\!\\!\int_{V}\\!\\!d{\bf
x}d{\bf x}^{\prime}\\!\\!\left({\hat{\rho}}({\bf
x})-\mbox{\boldmath$\nabla$}\\!\cdot\\!\hat{\cal\bf P}({\bf
x})\right)\tilde{u}({\bf x}\\!-\\!{\bf
x}^{\prime})\mbox{\boldmath$\nabla$}\\!\cdot\\!\hat{\cal\bf P}({\bf
x}^{\prime}),$ (72)

so that the last term would just rescale the potential $\tilde{u}_{ij}({\bf
x}-{\bf x}^{\prime})$ in Eq. 5

$\displaystyle\tilde{u}_{ij}({\bf x}-{\bf
x}^{\prime})\longrightarrow\tilde{u}_{ij}({\bf x}-{\bf
x}^{\prime})+\tilde{u}({\bf x}-{\bf
x}^{\prime})\mbox{\boldmath$\nabla$}_{i}\mbox{\boldmath$\nabla$}^{\prime}_{j},$
(73)

and if we take the non-electrostatic potential to the second order in
derivatives this simply means a rescaling of the corresponding constants,
which is irrelevant. The hydration energy can thus be simply taken as

$\int\\!\\!\int_{V}d{\bf x}d{\bf x}^{\prime}~{}\hat{\rho}({\bf
x})\tilde{u}({\bf x}-{\bf
x}^{\prime})\mbox{\boldmath$\nabla$}\cdot{\hat{\cal\bf P}}({\bf x}^{\prime}),$
(74)

which is what we use in the main text, Eq. (8). If the hydration shell of
anions and cations differ, the proper Ansatz would then be

$\displaystyle\int\\!\\!\int_{V}d{\bf x}d{\bf
x}^{\prime}~{}\hat{\rho}_{(+)}({\bf x})\tilde{u}_{+}({\bf x}-{\bf
x}^{\prime})\mbox{\boldmath$\nabla$}\cdot{\cal\bf\hat{P}}({\bf x}^{\prime})+$
$\displaystyle+\int\\!\\!\int_{V}d{\bf x}d{\bf
x}^{\prime}~{}\hat{\rho}_{(-)}({\bf x})\tilde{u}_{-}({\bf x}-{\bf
x}^{\prime})\mbox{\boldmath$\nabla$}\cdot{\cal\bf\hat{P}}({\bf x}^{\prime}),$
(75)

where the $\tilde{u}$’s are now the non-electrostatic hydration shell
potentials, different for different types of ions.

In the case of the mean-field approximation, the Model 3 saddle-point
equations are reduced to

$\displaystyle\phi^{*}({\bf x})=\int_{V}d{\bf x}^{\prime}~{}u({\bf x}-{\bf
x}^{\prime}){\rho}({\bf x}^{\prime})+$ (76) $\displaystyle\int_{V}d{\bf
x}^{\prime}~{}\tilde{u}({\bf x}-{\bf
x}^{\prime})\mbox{\boldmath$\nabla$}\cdot{{\cal\bf P}}({\bf x}^{\prime})$

and

$\displaystyle{\cal E}^{*}({\bf x})=-\alpha\mbox{\boldmath$\nabla$}{\rho}({\bf
x})+u_{P}(0){\bf P}({\bf x})-{\cal F}(\mbox{\boldmath$\nabla$}\cdot{\bf
P}({\bf x}))$ (77)

where ${\cal F}$ is considered up to fourth order. With

$\displaystyle\tilde{u}({\bf x}-{\bf x}^{\prime})=\alpha~{}\delta({\bf x}-{\bf
x}^{\prime})+\dots$ (78)

the expression for $\phi^{*}({\bf x})$ simplifies to

$\displaystyle\phi^{*}({\bf x})=\int_{V}d{\bf x}^{\prime}~{}u({\bf x}-{\bf
x}^{\prime}){\rho}({\bf
x}^{\prime})+\alpha\mbox{\boldmath$\nabla$}\cdot{{\cal\bf P}}({\bf x})\,.$
(79)

Applying the Laplacian on the last equation yields

$\displaystyle{\rho}({\bf
x})=\alpha\varepsilon_{\infty}\varepsilon_{0}\nabla^{2}\mbox{\boldmath$\nabla$}\cdot{{\cal\bf
P}}({\bf x})-\varepsilon_{\infty}\varepsilon_{0}\nabla^{2}\phi^{*}({\bf x}).$
(80)

Combining everything together into the final set of saddle-point equations for
the ion density polarization coupling case we remain with

${\cal E}^{*}({\bf x})=u_{P}(0)-{\cal F}(\mbox{\boldmath$\nabla$}\cdot{\bf
P}({\bf x}))$ (81)
$\displaystyle-\alpha^{2}\varepsilon_{\infty}\varepsilon_{0}\nabla^{2}\mbox{\boldmath$\nabla$}\left(\mbox{\boldmath$\nabla$}\cdot{\bf
P}({\bf
x})\right)+\alpha\varepsilon_{\infty}\varepsilon_{0}\nabla^{2}(\mbox{\boldmath$\nabla$}\phi^{*}({\bf
x}))\,.$

The $\alpha$-coupling therefore modifies the fourth-order coefficient of the
expansion of the non-electrostatic potential.

##  Appendix C Model parametrization

For numerical calculations of the solutions shown in the paper we employ
parameters chosen according to the ranges indicated in Table I.

Symbol | Meaning | Value
---|---|---
$1/\kappa_{D}$ | Debye length |  $\sim$ 10 nm
$\xi$ | bulk water correlation length | 1-10 nm
$\zeta$ | water molecular size | 0.3 nm
$\varepsilon$ | bulk dielectric constant | 5-80
$\varepsilon_{\infty}$ | dielectric constant first dispersion | 4-10
$\alpha$ | ion density-polarization coupling |  $\pm$ 0.5 - 5 nm
Table 1: Model parameters

The coupling parameter of the structural interactions,
$u_{P}(0)\varepsilon_{0}$, can be obtained from Eq. (26) in the main text,
slightly rewritten in order to make the dimensional dependencies explicit. It
then reads as

${\textstyle\left(\frac{\varepsilon}{\varepsilon_{0}}-\frac{\varepsilon_{\infty}}{\varepsilon_{0}}\right)}=\frac{{\textstyle\frac{\lambda
p^{2}}{3\varepsilon_{0}}}}{1+u_{P}(0)\varepsilon_{0}{\textstyle\frac{\lambda
p^{2}}{3\varepsilon_{0}}}},$ (82)

where $\lambda$ is the fugacity and $p$ the dipole moment in the dipolar
model. The form of this relation would remain unchanged if one considers the
dipolar Langevin model, but the values of the constants would have a different
meaning. We consider these differences as irrelevant at this stage, since we
are mainly interested here in qualitative features of the theory. By
dimensional analysis, one has

$[\lambda]=J^{-1}m^{-3}\,,\,\,\,[p]=Cm\,,\,\,\,[\varepsilon_{0}]=\frac{C^{2}}{Jm}\,.$
(83)

In order to connect the value for $u_{P}(0)\varepsilon_{0}$ with a measurable
macroscopic property, we invoke an additional relation

$\displaystyle{1+\frac{\lambda
p^{2}}{3\varepsilon_{0}}=\frac{\varepsilon}{\varepsilon_{0}}}$ (84)

which has been used [36]. For the standard values for water, $\varepsilon=80$,
$\varepsilon_{\infty}=4$, we have $u_{P}(0)\varepsilon_{0}=5\cdot 10^{-4}$.
The hatted parameters can then be computed from the bare lengths by
multiplication with a factor $\approx 0.038$.

##  Appendix D The analytic solution for the reduced Model 3 in the plate
geometry

The analytic solution for the reduced Model 3 in the slab geometry can be
obtained by making recourse to the SI of [41]. Defining

$\delta_{c}^{2}\equiv\left({\textstyle\frac{\varepsilon_{\infty}}{\varepsilon}}\right)^{3}\widehat{\zeta_{0}}^{4}\left(2\widehat{\alpha}^{2}({\textstyle\frac{\varepsilon}{\varepsilon_{\infty}}}-1)+\widehat{\xi}^{2}\right)^{2}$
(85)

one has, for $\delta_{c}>1/2$, the exact solution in the slab geometry given
by

$\phi^{\dagger}(z)=P_{c}e^{-\kappa_{1}z}[\cos(\kappa_{2}z)+A\sin(\kappa_{2}z)]$
(86)

with

$\kappa_{1}=\frac{\sqrt{2\delta_{c}+1}}{2\delta_{c}}\,,\,\,\kappa_{2}=\frac{\sqrt{2\delta_{c}-1}}{2\delta_{c}}$
(87)

and

$A=-\frac{\sqrt{2\delta_{c}+1}(\delta_{c}-1)}{\sqrt{2\delta_{c}-1}(\delta_{c}+1)}\,.$
(88)

From this solution, $\phi^{*}(z)$ follows via the expression

$\phi^{*}(z)=-\left(\frac{\varepsilon}{\varepsilon_{\infty}}-1\right)\left(1-\widehat{\alpha}^{2}\gamma\frac{d^{2}}{dz^{2}}\right)\phi^{\dagger}(z)\,$
(89)

with

$\gamma\equiv{\textstyle\left(2\widehat{\alpha}^{2}\left(\frac{\varepsilon}{\varepsilon_{\infty}}-1\right)+\widehat{\zeta}^{2}\right)^{-1}}\,.$
(90)

## References

  * [1] R. A. Marcus, On the theory of oxidation and reduction reactions involving electron transfer. i, The Journal of Chemical Physics 24, 966 (1956).
  * [2] R. A. Marcus, Electrostatic free energy and other properties of states having nonequilibrium polarization. i, The Journal of Chemical Physics 24, 979 (1956).
  * [3] C. Holm, P. Kekicheff, and R. Podgornik, Electrostatic effects in soft matter and biophysics, 2nd ed., NATO science series, Series II, Mathematics, physics, and chemistry, Vol. 46 (Kluwer Academic Publishers, Dordrecht, Boston, 2001).
  * [4] P. A. Bopp, A. A. Kornyshev, and G. Sutmann, Static nonlocal dielectric function of liquid water, Phys. Rev. Lett. 76, 1280 (1996).
  * [5] A. A. Kornyshev, A. I. Rubinshtein, and M. A. Vorotyntsev, Model nonlocal electrostatics. i, 11, 3307 (1978).
  * [6] A. A. Kornyshev, On the non-local electrostatic theory of hydration force, Journal of Electroanalytical Chemistry and Interfacial Electrochemistry 204, 79 (1986).
  * [7] A. C. Maggs and R. Everaers, Simulating nanoscale dielectric response, Phys. Rev. Lett. 96, 230603 (2006).
  * [8] E. L. Mertz, Conditions for insensitivity of the microscopic-scale dielectric response to structural details of dipolar liquids, Phys. Rev. E 76, 062503 (2007).
  * [9] H. Berthoumieux and F. Paillusson, Dielectric response in the vicinity of an ion: A nonlocal and nonlinear model of the dielectric properties of water, The Journal of Chemical Physics 150, 094507 (2019).
  * [10] H. Berthoumieux, G. Monet, and R. Blossey, Dipolar Poisson models in a dual view, The Journal of Chemical Physics 155, 024112 (2021).
  * [11] G. Monet, F. Bresme, A. Kornyshev, and H. Berthoumieux, Nonlocal dielectric response of water in nanoconfinement, Phys. Rev. Lett. 126, 216001 (2021).
  * [12] L. Onsager and M. Dupuis, The electrical properties of ice, in Rend. Sc. Int. Fis. ‘Enrico Fermi”, Corso X Varenna 1959 (N. Zanichelli, Bologna, 1960) pp. 294-315.
  * [13] L. Onsager and M. Dupuis, The electrical properties of ice, in Electrolytes, edited by B. Pesce (Pergamon, Oxford, 1962) pp. 27-46.
  * [14] D. W. R. Gruen and S. Marčelja, Spatially varying polarization in ice, J. Chem. Soc. Faraday Trans. 2 79, 211 (1983).
  * [15] D. W. R. Gruen and S. Marčelja, Spatially varying polarization in water. A model for the electric double layer and the hydration force, J. Chem. Soc. Faraday Trans. 2 79, 225 (1983).
  * [16] S. Leikin, V. A. Parsegian, D. C. Rau, and R. P. Rand, Hydration forces, Annual Review of Physical Chemistry 44, 369 (1993).
  * [17] V. Parsegian and T. Zemb, Hydration forces: Observations, explanations, expectations, questions, Current Opinion in Colloid $\&$ Interface Science 16, 618 (2011).
  * [18] D. Gruen, S. Marčelja, and B. Pailthorpe, Theory of polarization profiles and the “hydration force”, Chemical Physics Letters 82, 315 (1981).
  * [19] M. Belaya, M. Feigel’man, and V. Levadny, Hydration forces as a result of non-local water polarizability, Chemical Physics Letters 126, 361 (1986).
  * [20] S. Marčelja, and N.Radić, Repulsion of interfaces due to boundary water, Chemical Physics Letters 42, 129 (1976).
  * [21] D. C. Rau, B. Lee, and V. A. Parsegian, Measurement of the repulsive force between polyelectrolyte molecules in ionic solution: hydration forces between parallel DNA double helices, 81, 2621 (1984).
  * [22] D. C. Rau and V. A. Parsegian, Direct measurement of forces between linear polysaccharides xanthan and schizophyllan, Science 249, 1278 (1990).
  * [23] A. Hildebrandt, R. Blossey, S. Rjasanow, O. Kohlbacher and H.-P. Lenhof, Novel formulation of nonlocal electrostatics, Phys. Rev. Lett. 93, 108104 (2004).
  * [24] A. Hildebrandt, R. Blossey, S. Rjasanow, O. Kohlbacher, and H.-P. Lenhof, Electrostatic potentials of proteins in water: a structured continuum approach, Bioinformatics 23, e99 (2007).
  * [25] L. Fumagalli, A. Esfandiar, R. Fabregas, S. Hu, P. Ares, A. Janardanan, Q. Yang, B. Radha, T. Taniguchi, K. Watanabe, G. Gomila, K. S. Novoselov, and A. K. Geim, Anomalously low dielectric constant of confined water, Science 360, 1339 (2018).
  * [26] E. Papadopoulou, J. Zavadlav, R. Podgornik, M. Praprotnik, and P. Koumoutsakos, Tuning the dielectric response of water in nanoconfinement through surface wettability, ACS Nano 0, 0 (2021).
  * [27] M.Kanduc,A.Schlaich,E.Schneck,andR.R.Netz, Hydration repulsion between membranes and polar surfaces: simulation approaches versus continuum theories, Adv Colloid Interface Sci 208, 142 (2014).
  * [28] A. Schlaich, E. W. Knapp, and R. R. Netz, Water dielectric effects in planar confinement, Phys. Rev. Lett. 117, 1 (2016).
  * [29] V. S. Mitlin and M. M. Sharma, A local gradient theory for structural forces in thin fluid films, Journal of Colloid and Interface Science 157, 447 (1993).
  * [30] S. Varghese, S. K. Kannam, J. S. Hansen, and S. P. Sathian, Effect of hydrogen bonds on the dielectric properties of interfacial water, Langmuir 35, 8159 (2019).
  * [31] R. Podgornik, Theory of inhomogeneous rod-like Coulomb fluids, Symmetry 13 (2021).
  * [32] S. F. Edwards and A. Lenard, Exact statistical mechanics of a one-dimensional system with Coulomb forces. ii. the method of functional integration, Journal of Mathematical Physics 3, 778 (1962).
  * [33] R. Podgornik, Electrostatic correlation forces between surfaces with surface specific ionic interactions, The Journal of Chemical Physics 91, 5840 (1989).
  * [34] A. Levy, D. Andelman, and H. Orland, Dipolar Poisson-Boltzmann approach to ionic solutions: A mean field and loop expansion analysis, The Journal of Chemical Physics 139, 164909 (2013).
  * [35] R. M. Adar, T. Markovich, A. Levy, H. Orland, and D. Andelman, Dielectric constant of ionic solutions: Combined effects of correlations and excluded volume, The Journal of Chemical Physics 149, 054504 (2018).
  * [36] A. Abrashkin, D. Andelman, and H. Orland, Dipolar Poisson-Boltzmann equation: Ions and dipoles close to charge interfaces, Phys. Rev. Lett. 99, 077801 (2007).
  * [37] F. Paillusson and R. Blossey, Slits, plates, and Poisson-Boltzmann theory in a local formulation of nonlocal electrostatics, Phys. Rev. E 82, 052501 (2010).
  * [38] E. Ruckenstein and M. Manciu, The coupling between the hydration and double layer interactions, Langmuir 18, 7584 (2002).
  * [39] R. I. Slavchov, Quadrupole terms in the Maxwell equations: Debye-Hu?ckel theory in quadrupolarizable solvent and self-salting-out of electrolytes, The Journal of Chemical Physics 140, 164510 (2014).
  * [40] R. I. Slavchov and T. I. Ivanov, Quadrupole terms in the Maxwell equations: Born energy, partial molar volume, and entropy of ions, The Journal of Chemical Physics 140, 074503 (2014).
  * [41] M. Z. Bazant, B. D. Storey, and A. A. Kornyshev, Double layer in ionic liquids: Overscreening versus crowding, Phys. Rev. Lett. 106, 046102 (2011).


# Computation with Large Advice

Hiroki Morizumi
Shimane University, Japan
<EMAIL_ADDRESS>
###### Abstract

In this paper, we consider a new direction of computation, which we call
computation with large advice. We mainly consider constant space computation
with large advice in Turing machines, and prove the following facts.

  * •

The class of decision problems solvable by a constant space Turing machine
with polynomial-size advice includes nonuniform-NC1.

  * •

The class of decision problems solvable by a constant space Turing machine
with quasipolynomial-size advice equals nonuniform-polyL.

The facts mean constant space computation with large advice has an unexpected
computational power. On the other hand, we attempt to propose a concept of
“algorithms with large advice”. In the proposal, advice is precomputed data
for a problem and a fixed instance size, and we expect efficient algorithms by
large or huge advice.

##  1 Introduction

In this paper, we prove the following theorems.

######  Theorem 1.

The class of decision problems solvable by a constant space Turing machine
with polynomial-size advice includes nonuniform-NC1.

######  Theorem 2.

The class of decision problems solvable by a constant space Turing machine
with quasipolynomial-size advice equals nonuniform-polyL.

The advice of Turing machines is an extra input which depends only on the size
of the main input. Two theorems imply unexpected power of large advice in
computation.

###  1.1 The Contents of This Paper

The main results of this paper are Theorem 1 and Theorem 2. In Section 3, we
describe our ideas and contribution.

Theorem 1 and Theorem 2 raise questions and interests for the power of large
advice in computation. In this paper, we also mention the power of large
advice for algorithms. The two theorems do not directly imply usefulness of
algorithms with large advice, and we mention a conceptual proposal toward
future studies of algorithms. See Section 6 for the details.

##  2 Definitions

###  2.1 Branching Programs

A branching program is a directed acyclic graph. The nodes of out-degree 2 are
called inner nodes and labeled by a variable. The nodes of out-degree 0 are
called sinks and labeled by 0 or 1. For each inner node, one of the outgoing
edges is labeled by 0 and the other one is labeled by 1. There is a single
specific node called the start node. An assignment to the variables determines
a computation path from the start node to a sink node. The value of the sink
node is the output of the branching program. If the nodes are arranged into a
sequence of levels with edges going only from one level to the next, then the
width is the size of the largest level.

###  2.2 Boolean Circuits

Circuits are formally defined as directed acyclic graphs. The nodes of in-
degree 0 are called inputs, and each one of them is labeled by a variable or
by a constant 0 or 1. The other nodes are called gates, and each one of them
is labeled by a Boolean function. The fan-in of a node is the in-degree of the
node, and the fan-out of a node is the out-degree of the node. In this paper,
the gates are AND gates of fan-in two, OR gates of fan-in two, and NOT gates.
There is a single specific node called output. The size of a circuit is the
number of gates in the circuit.

###  2.3 Complexity Classes

NC1 is the class of decision problems solvable by a uniform family of Boolean
circuits with polynomial size, depth $O(\log n)$, and fan-in 2. polyL is the
class of decision problems solvable by a polylogarithmic space Turing machine.

###  2.4 The Sorting Problem

The sorting problem is defined as follows in this paper.

Input: $n$ numbers in the set $\\{0,1,\ldots,k\\}$.

Output: A permutation $\langle a_{1},a_{2},\ldots,a_{n}\rangle$ of the input
such that $a_{1}\leq a_{2}\leq\cdots\leq a_{n}$.

Note that $n$ numbers are restricted in the set $\\{0,1,\ldots,k\\}$.

##  3 The Key Ideas and Our Contribution

Our results of this paper are closely related to Barrington’s theorem which
has been well-known in computational complexity theory. In the 1980’s,
Barrington proved the following theorem.

######  Theorem 3 ([1]).

The class of decision problems solvable by a nonuniform family of polynomial-
size $5$-width branching programs equals nonuniform-NC1.

The theorem is a surprising result, since $5$-width branching programs have
only at most five states in each stage of computation. Thus, it has been known
that computation with constant working space can have an unexpected
computational power. In this paper, we represent the computation in
Barrington’s theorem by Turing machines with large advice, which is our main
contribution and gives a new insight to computation.

##  4 Proof of Theorem 1

We firstly prove a relation between constant space Turing machines with
polynomial-size advice and polynomial-size $5$-width branching programs.

######  Lemma 4.

Constant space Turing machines with polynomial-size advice simulate
polynomial-size $5$-width branching programs.

###### Proof.

We prepare an encoding of the branching program on the tape of advice. The
execution of the branching program is simulated with sequential access to the
encoding. The encoding for variables of branching programs needs to be
careful. The encoding for variable $x_{i}$ consists of one starting symbol and
$n$ symbols, the $i$’th of which is marked. ∎

We prove Theorem 1.

###### Proof of Theorem 1.

By Theorem 3 and Lemma 4, the theorem holds. ∎

##  5 Proof of Theorem 2

The proof outline of Theorem 2 is similar to the one of Theorem 1. To prove
Theorem 2, we prove the quasipolynomial-size variant of Barrington’s theorem.

###  5.1 The Quasipolynomial-Size Variant of Barrington’s Theorem

The following theorem is the quasipolynomial-size variant of Barrington’s
theorem.

######  Theorem 5.

The class of decision problems solvable by a nonuniform family of
quasipolynomial-size $5$-width branching programs equals nonuniform-polyL.

To prove Barrington’s theorem, Barrington has proved the following lemma. The
lemma is useful also for the proof of Theorem 5.

######  Lemma 6 (Theorem 1 of [1] (with some modifications)).

Any circuit with depth $d$ and fan-in 2 can be converted to a $5$-width
branching program with length at most $4^{d}$.

We also need the following lemma.

######  Lemma 7.

Any quasipolynomial-size branching program can be converted to a circuit of
polylogarithmic depth.

###### Proof.

By a standard conversion to a parallel execution, the lemma holds. ∎

We prove Theorem 5.

###### Proof of Theorem 5.

By Lemma 6 and Lemma 7, any quasipolynomial-size branching program can be
converted to a quasipolynomial-size $5$-width branching program. Since the
class of decision problems solvable by a nonuniform family of quasipolynomial-
size branching programs equals nonuniform-polyL, the theorem holds. ∎

###  5.2 Proof of Theorem 2

This subsection is almost similar to Section 4.

######  Lemma 8.

Constant space Turing machines with quasipolynomial-size advice simulate
quasipolynomial-size $5$-width branching programs.

###### Proof.

We can prove the lemma by the similar way of Lemma 4. ∎

We prove Theorem 2.

###### Proof of Theorem 2.

Nonuniform-polyL obviously includes the class of decision problems solvable by
a constant space Turing machine with quasipolynomial-size advice. By Theorem 5
and Lemma 8, the theorem holds. ∎

##  6 Algorithms with Large Advice

In this section, we attempt to propose a concept of “algorithms with large
advice”. In our proposal, advice is precomputed data for a problem and a fixed
instance size, and we expect efficient algorithms by large or huge advice. (We
call the data advice, since the data correspond to advice of Turing machines.)
In the following example, the merge sort algorithm for the sorting problem
accelerates from $O(n\log n)$ time to $O(n\log\log n)$ time by huge advice.

######  Theorem 9.

A modified merge sort algorithm runs in $O(n\log\log n)$ time, if $n$ is fixed
and $\tilde{O}(k^{\frac{n}{\log n}})$ advice (i.e., precomputed data of
$\tilde{O}(k^{\frac{n}{\log n}})$ space) is given.

###### Proof.

We execute the $\log\log n$ usual recursive steps of the merge sort, and sort
$\frac{n}{\log n}$ numbers by precomputed data. Precomputed data has all
patterns of sequences of $\frac{n}{\log n}$ numbers and the sorted sequences.
∎

The sorting problem is an example for the framework of algorithms with large
advice. We are interested in algorithms for computationally hard problems such
as NP-hard problems, e.g., the SAT problem.

Some remarks. In one of situations which we suppose, advice (i.e., precomputed
data) is stored in some places of the internet as public data. The current
developed internet offers good environment for algorithms with large advice.
In this paper, we described a minimal framework for algorithms with large
advice. For actual execution, we need an efficient algorithm to compute advice
(i.e., precomputed data), although the algorithm runs once and the obtained
advice can be stored.

## References

  * [1] David A. Mix Barrington.  Bounded-width polynomial-size branching programs recognize exactly those languages in nc1.  J. Comput. Syst. Sci., 38(1):150–164, 1989.