# Contribution of dipolar bridging to phospholipid membrane interactions: a
mean-field analysis

Sahin Buyukdagli1111email<EMAIL_ADDRESS>and Rudolf
Podgornik2,3,4222email<EMAIL_ADDRESS>
Also affiliated: Department of Physics, Faculty of Mathematics and Physics,
University of Ljubljana, SI-1000 Ljubljana, Slovenia. 1Department of Physics,
Bilkent University, Ankara 06800, Turkey
2School of Physical Sciences and Kavli Institute for Theoretical Sciences,
University of Chinese Academy of Sciences, Beijing 100049, China
3CAS Key Laboratory of Soft Matter Physics, Institute of Physics, Chinese
Academy of Sciences (CAS), Beijing 100190, China
4 Wenzhou Institute of the University of Chinese Academy of Sciences, Wenzhou,
Zhe-jiang 325000, China.

###### Abstract

We develop a model of interacting zwitterionic membranes with rotating surface
dipoles immersed in a monovalent salt, and implement it in a field theoretic
formalism. In the mean-field regime of monovalent salt, the electrostatic
forces between the membranes are characterized by a non-uniform trend: at
large membrane separations, the interfacial dipoles on the opposing sides
behave as like-charge cations and give rise to repulsive membrane
interactions; at short membrane separations, the anionic field induced by the
dipolar phosphate groups sets the behavior in the intermembrane region. The
attraction of the cationic nitrogens in the dipolar lipid headgroups leads to
the adhesion of the membrane surfaces via dipolar bridging. The underlying
competition between the opposing field components of the individual dipolar
charges leads to the non-uniform salt ion affinity of the zwitterionic
membrane with respect to the separation distance; large inter-membrane
separations imply anionic excess while small, nanometer size separations,
favor cationic excess. This complex ionic selectivity of zwitterionic
membranes may have relevant repercussions on nanofiltration and nanofluidic
transport techniques.

###### pacs:

41.20.Cv,82.45.Gj,87.16.Dg

##  I Introduction

Macro-molecular forces driven by electrostatics play a primary role in the
stability of living matter as well as the regulation of the underlying
biological processes Isr ; PodgornikRev ; biomatter . From the aggregation of
like-charged macromolecules Molina2014 ; exp2 ; Cruz  to the emergence of
anionic streaming currents through negatively charged pores Heyden2006  and
anionic polymer mobility along the applied electric field Qiu2015 ,
electrostatic interactions have been also at the origin of numerous seemingly
counterintuitive phenomena. Among these unconventional electrostatic effects,
the repulsion between overall neutral dipolar membranes has been identified
first in Parsegian’s experiments with egg lecithin bilayers as hydration
interaction Parsegian1979 ; Rand1981 . Subsequently, these experiments have
been extended to understand the effect of monovalent salt McIntosh1990 ;
Parsegian1991 ; Kekicheff2014 , divalent buffers Petrache2011 , and
multivalent ions Averbakh2000 ; Fink2020  on zwitterionic membrane
interactions.

Early theoretical analyses proposed the interactions mediated by perturbed
structural order of water molecules Marcelja  leading to repulsive hydration
interactions Leikin . Later, the underlying mechanism behind these repulsive
hydration forces was surmized to be the image dipoles associated with the
interfacial zwitterionic charges Jonsson1983 ; Kjellander1984 . The structural
solvation interaction was further detailed in an effective solvent model based
on a non-local dielectric response formalism that explained the dipolar
membrane repulsion in terms of the coupling between the solvent and surface
dipoles Belaya1 ; Belaya2 . A thermodynamic model of hydration interactions
was formulated for complex lipid bilayers with dipolar headgroups based on
several spatially separated charges and internal conformational degrees of
freedom in Ref. Bohinc2017 . Finally, recent solvent explicit simulations
indicated that although the direct electrostatic interaction of the dipolar
charges on the opposing membrane walls is attractive, this attraction is
reversed by water-mediated hydration forces contributing the major part to the
zwitterionic membrane repulsion sim1 ; sim2 ; sim3 . These studies also showed
that various additional effects such as the entropy of the interacting lipid
headgroups associated with their annealed degrees of freedom play a non-
negligible role in the inter-membrane coupling.

Due to the complexity of the interaction picture emerging from the simulation
results, numerical studies have to be completed by simpler theoretical
formalisms including the essential details of the physical system. Motivated
by this point, we develop a new model of zwitterionic membranes and implement
it in a field theoretic formalism to study the electrostatic mean-field (MF)
interaction regime in monovalent salt bathing solutions, where the model
provides us with analytically transparent predictions.

We introduce a dipolar membrane model in Sec. II, composed of an electrolyte
confined by two membrane walls coated by rotating surface dipoles. By
including this orientational degree of freedom, the present formalism
generalizes our earlier work on fixed dipolar membrane interactions to
realistic zwitterionic charges DipFix . In the model, the electrolyte ions in
chemical equilibrium with the external bulk reservoir are treated in the
grand-canonical ensemble, while the surface dipoles of fixed number are
incorporated within the canonical ensemble. We then derive the saddle-point
equation within the field theoretic formalism and establish the corresponding
form of the contact value theorem, that provides us with the MF-level
disjoining pressure. Based upon this formalism we characterize in Sec. III the
interaction of zwitterionic membranes in MF conditions of monovalent salt and
reveal a membrane complexation mechanism driven by the competition between
large separation repulsive membrane interactions and short separation
zwitterionic charge attractions mediated by dipolar bridging interactions. Our
main results are summarized and potential extensions of the zwitterionic
membrane model are discussed in Conclusions.

##  II Theory and Model

###  II.1 Field theoretic model

Figure 1: (Color online) Side view of the lipid membrane model including the
slit pore between two membranes at separation $d$. The surface dipoles (D) are
separated by the lateral distance $l$, and located at the distance $p$ from
the membrane surface of total area $S$. The membrane and electrolyte
permittivities are $\varepsilon_{\rm m}$ and $\varepsilon_{\rm w}$,
respectively.

In this part, we introduce our dipolar membrane model depicted in Fig. 1. The
two membrane surfaces are located in the $x-y$ plane and are separated by the
distance $d$ along the $z$ axis. The surfaces connected to an ion reservoir
confine an electrolyte mixture at temperature $T=300$ K. The electrolyte is
composed of $s$ ionic species. The ions of the species $n$ with total number
$N_{n}$ have valency $q_{n}$, fugacity $\lambda_{n}$, and reservoir
concentration $\rho_{{\rm b}n}$. The system is characterized by the dielectric
permittivity profile

$\varepsilon(\mathbf{r})=\varepsilon_{\rm
m}\left[\theta(-z)+\theta(z-d)\right]+\varepsilon_{\rm
w}\theta(z)\theta(d-z),$ (1)

where $\varepsilon_{\rm m}$ and $\varepsilon_{\rm w}=80$ stand for the
dielectric permittivity of the membrane and solvent, respectively.

Each membrane wall of total area $S$ carries surface dipoles (D). The charge
structure of each dipole is described by an anion ($-Q_{-}$), fixed at the
separation distance $p$ from the wall, and a cation ($Q_{+}$) at a distance
$a$ from the anion, free to rotate around the anion at a fixed $a$. The wall
distribution of the dipoles is illustrated in Fig. 2(a). The dipoles separated
from each other by the lateral distance $l$ are located at the discrete mesh
points of square symmetry, numbered by the indices $i$ and $j$ along the $x$
and $y$ axis, respectively. Each row or column on the membrane surfaces
contains a total of $N$ dipoles, i.e. $1\leq i\leq N$ and $1\leq j\leq N$.
Thus, the position vectors $\mathbf{r}_{ij}$ and $\mathbf{r}^{\prime}_{ij}$ of
the dipoles respectively located on the left and right surfaces and
characterized by the solid angles $\mathbf{\Omega}_{ij}$ and
$\mathbf{\Omega}^{\prime}_{ij}$ read

$\displaystyle\mathbf{r}_{ij}$ $\displaystyle=$
$\displaystyle(i-1)l\hat{u}_{x}+(j-1)l\hat{u}_{y}+p\hat{u}_{z},$ (2)
$\displaystyle\mathbf{r}^{\prime}_{ij}$ $\displaystyle=$
$\displaystyle(i-1)l\hat{u}_{x}+(j-1)l\hat{u}_{y}+(d-p)\hat{u}_{z}.$ (3)

The total charge density operator is

$\displaystyle\hat{\rho}_{\rm c}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{s}q_{n}\sum_{m=1}^{N_{n}}\delta(\mathbf{r}-\mathbf{R}_{nm})$
$\displaystyle+\sum_{i=1}^{N}\sum_{j=1}^{N}\left[\hat{F}\left(\mathbf{r};\mathbf{r}_{ij},\mathbf{\Omega}_{ij}\right)+\hat{F}\left(\mathbf{r};\mathbf{r}^{\prime}_{ij},\mathbf{\Omega}^{\prime}_{ij}\right)\right],$

where $\mathbf{R}_{nm}$ stands for the mobile ion coordinate, and
$\hat{F}\left(\mathbf{r};\mathbf{r}_{ij},\mathbf{\Omega}_{ij}\right)$ is the
charge structure factor of the dipole $(i,j)$. The canonical partition
function of the system reads

$Z_{\rm
c}=\prod_{n=1}^{s}\prod_{m=1}^{N_{n}}\int\mathrm{d}^{3}\mathbf{R}_{nm}\prod_{i=1}^{N}\prod_{j=1}^{N}\int{\color[rgb]{0,0,0}\frac{\mathrm{d}^{2}\mathbf{\Omega}_{ij}}{4\pi}\frac{\mathrm{d}^{2}\mathbf{\Omega}^{\prime}_{ij}}{4\pi}}e^{-\beta
E_{\rm c}-\beta E_{\rm n}},$ (5)

where the electrostatic interaction energy is given by

$\beta E_{\rm
c}=\frac{1}{2}\int\mathrm{d^{3}}\mathbf{r}\mathrm{d^{3}}\mathbf{r}^{\prime}\rho_{\rm
c}(\mathbf{r})v_{\rm c}(\mathbf{r},\mathbf{r}^{\prime})\rho_{\rm
c}(\mathbf{r}^{\prime}),$ (6)

In Eq. (6), the Coulomb Green’s function $v_{\rm
c}(\mathbf{r},\mathbf{r}^{\prime})$ is defined in terms of its inverse,
$v_{\rm c}^{-1}(\mathbf{r},\mathbf{r}^{\prime})=-(k_{\rm
B}T/e^{2})\nabla\cdot\varepsilon(\mathbf{r})\nabla\delta(\mathbf{r}-\mathbf{r}^{\prime})$.
Moreover, In Eq. (5), the total steric energy of the ions and dipoles is
defined as

$\displaystyle\beta E_{\rm n}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{s}\sum_{m=1}^{N_{n}}V_{\rm n}(\mathbf{R}_{nm})$
$\displaystyle+\sum_{i=1}^{N}\sum_{j=1}^{N}\left[V_{\rm
l}(\mathbf{r}_{ij},\mathbf{\Omega}_{ij})+V_{\rm
r}(\mathbf{r}^{\prime}_{ij},\mathbf{\Omega}^{\prime}_{ij})\right],$

where we introduced the ionic onsite potential $V_{\rm n}(\mathbf{R}_{nm})$,
and the dipolar onsite potentials $V_{\rm
l}(\mathbf{r}_{ij},\mathbf{\Omega}_{ij})$ and $V_{\rm
r}(\mathbf{r}^{\prime}_{ij},\mathbf{\Omega}^{\prime}_{ij})$, with the indices
${\rm l}$ and ${\rm r}$ referring to the dipoles located on the left and right
walls, respectively. These one-body potentials will be used to generate the
ionic charge densities and the conditional probabilities of the surface
dipoles.

At this point, performing an Hubbard-Stratonovich transformation, we introduce
the fluctuating electrostatic potential $\psi(\mathbf{r})$ that allows us to
linearize Eqs. (5)-(6) in the charge density (II.1). Then, via the replacement
$\sum_{i=1}^{N}\sum_{j=1}^{N}\to\sigma_{\rm
s}\int\mathrm{d}^{2}\mathbf{r}_{\parallel}$, we pass to the continuously
distributed dipole limit, with the position vector
$\mathbf{r}_{\parallel}=x\hat{u}_{x}+y\hat{u}_{y}$ within the membrane
surface, and the surface dipole density $\sigma_{\rm s}=l^{-2}$. Finally,
taking into account the ion exchange driven by the chemical equilibrium
between the slit region and the reservoir, the grand-canonical partition
function of the system takes the functional integral form

$Z_{\rm
G}=\prod_{n=1}^{s}\sum_{N_{n}=1}^{\infty}\frac{\lambda_{n}^{N_{n}}}{N_{n}!}Z_{\rm
c}=\int\mathcal{D}{\color[rgb]{0,0,0}\psi}\;e^{-\beta
H[{\color[rgb]{0,0,0}\psi}]},$ (8)

with the Hamiltonian functional

$\displaystyle\beta H[{\color[rgb]{0,0,0}\psi}]$ $\displaystyle=$
$\displaystyle\frac{k_{\rm
B}T}{2e^{2}}\int\mathrm{d}^{3}\mathbf{r}\;\varepsilon(\mathbf{r})\left[\nabla{\color[rgb]{0,0,0}\psi}(\mathbf{r})\right]^{2}$
$\displaystyle-\sum_{n=1}^{s}\lambda_{n}\int\mathrm{d}^{3}\mathbf{r}e^{iq_{n}{\color[rgb]{0,0,0}\psi}(\mathbf{r})-V_{n}(\mathbf{r})}$
$\displaystyle-\sigma_{\rm
s}\int\mathrm{d}^{3}\mathbf{r}\left[\delta(z-p)+\delta(z-d+p)\right]$
$\displaystyle\hskip
19.91692pt\times\ln\int\frac{\mathrm{d}^{2}\mathbf{\Omega}}{4\pi}e^{i\int\mathrm{d}^{3}\mathbf{r}^{\prime}\hat{F}\left(\mathbf{r}^{\prime};\mathbf{r},\mathbf{\Omega}\right){\color[rgb]{0,0,0}\psi}(\mathbf{r}^{\prime})-V_{\rm
d}(\mathbf{r},\mathbf{\Omega})}.$

In Eq. (II.1), the dipolar potential is defined as $V_{\rm
d}(\mathbf{r},\mathbf{\Omega})=V_{\rm l}(\mathbf{r},\mathbf{\Omega})$ if $z=p$
and $V_{\rm d}(\mathbf{r},\mathbf{\Omega})=V_{\rm
r}(\mathbf{r},\mathbf{\Omega})$ if $z=d-p$.

Figure 2: (Color online) (a) Lateral view of the lipid membrane surfaces lying
in the $x-y$ plane, and located at $z=0$ and $z=d$. The surface dipoles (D)
separated by the distance $l$ are located at the discrete mesh points numbered
by the indices $i$ and $j$ along the $x$ and $y$ axis, respectively. (b) The
inner charge structure of the surface dipoles of length $a$, with an annealed
coordinate of the cation that is allowed to rotate about the anionic charge
fixed at $z=p$, with $a\leq p$ and $d\geq 2p+a$.

###  II.2 Derivation of the saddle-point equations

Fig. 2(b) displays the inner charge structure of the rotating membrane dipoles
with total length $a$ and terminal charges $\pm Q_{\pm}$. The cationic charge
$Q_{+}$ stands for the valency of the nitrogen (${\rm N}^{+}$) located at the
trimethyl ammonium moiety of the lipid headgroup, and the fixed anionic charge
$-Q_{-}$ is the valency of the phosphate group. The charge structure factor of
the surface dipole rotating about the anionic charge located at $\mathbf{r}$
reads

$\hat{F}(\mathbf{r}^{\prime};\mathbf{r},\mathbf{\Omega})=-Q_{-}\delta(\mathbf{r}^{\prime}-\mathbf{r})+Q_{+}\delta(\mathbf{r}^{\prime}-\mathbf{r}-{\boldsymbol{a}}).$
(10)

In order to prevent the rotating cations from touching the membrane surface
and the tail of the dipoles on the opposing membrane surface, we will restrict
ourselves to the tail length and pore size regimes of $p\geq a$ and $d\geq
2p+a$. Moreover, with the aim to prohibit the contact of the neighbouring
dipoles located on the same wall, we will limit our study to the dipolar
distance regime of $l>2a$. For the characteristic dipole size $a=2$ Å
considered in our work, this assumption sets for the surface dipole density
the reasonable constraint of $\sigma_{\rm s}\lesssim 6$ ${\rm nm}^{-2}$.

Within the MF approximation, the ion density follows from the thermodynamic
relation $\rho_{n}(\mathbf{r})=-\delta\ln Z_{\rm G}/\delta
V_{n}(\mathbf{r})\approx\delta H/\delta V_{n}(\mathbf{r})$ as
$\rho_{n}(\mathbf{r})=\lambda_{n}e^{iq_{n}{\color[rgb]{0,0,0}\psi}(\mathbf{r})-V_{n}(\mathbf{r})}$.
In the bulk reservoir where
${\color[rgb]{0,0,0}\psi}(\mathbf{r})=V_{n}(\mathbf{r})=0$, this relates the
fugacity to the bulk concentration as $\lambda_{n}=\rho_{\rm nb}$. Passing to
the real dimensionless potential $\phi(\mathbf{r})$ via the transformation
$\psi(\mathbf{r})\to i\phi(\mathbf{r})$, and accounting for the planar
symmetry $\phi(\mathbf{r})=\phi(z)$, the average ion density becomes

$\rho_{n}(z)=\rho_{n{\rm b}}e^{-q_{n}\phi(z)-V_{n}(z)}.$ (11)

In order to simplify the notation, from now on, the ionic and dipolar steric
potentials will be omitted. Plugging now the dipolar structure factor (10)
into Eq. (II.1), and evaluating the saddle-point condition $\delta
H[\phi]/\delta\phi(z)=0$, after some algebra, one obtains the zwitterionic
Poisson-Boltzmann (ZPB) Eq. in the form

$\displaystyle\frac{k_{\rm
B}T}{e^{2}}\partial_{z}{\color[rgb]{0,0,0}\left[\varepsilon(z)\partial_{z}\phi(z)\right]}+\sum_{n=1}^{s}q_{n}\rho_{n}(z)+Q_{+}\left[\rho^{+}_{\rm
l}(z)+\rho^{+}_{\rm r}(z)\right]$ $\displaystyle=\sigma_{\rm
s}Q_{-}\left[\delta(z-p)+\delta(z-d+p)\right].$ (12)

In Eq. (II.2), we introduced the dipolar cation densities

$\rho^{+}_{\rm l,r}(z)=\frac{\sigma_{\rm s}}{2aI_{\rm
l,r}}e^{-Q_{+}\phi(z)}\theta_{\rm l,r}(z)$ (13)

derived in Appendix A. The auxiliary functions

$\displaystyle\theta_{\rm l}(z)$ $\displaystyle=$
$\displaystyle\theta(z-p+a)\theta(a-z+p),$ (14) $\displaystyle\theta_{\rm
r}(z)$ $\displaystyle=$ $\displaystyle\theta(z-d+p+a)\theta(d-p+a-z),$ (15)

including the Heaviside step function $\theta(z)$ account for the finite
extension of the interfacial layer scanned by the rotating dipoles. Then, the
integrals incorporating the constraint on the fixed lateral dipole density
read

$I_{\rm
l}=\int_{-a}^{a}\frac{\mathrm{d}a_{z}}{2a}e^{-Q_{+}\phi(p+a_{z})};\;I_{\rm
r}=\int_{-a}^{a}\frac{\mathrm{d}a_{z}}{2a}e^{-Q_{+}\phi(d-p+a_{z})},$ (16)

with the dipole projection $a_{z}=a\cos\theta$. Using the mirror symmetry with
respect to the mid-pore, one can show that these integrals are equal, $I_{\rm
l}=I_{\rm r}$. Therefore, when the specification of the membrane location is
not necessary, these integrals will be denoted by the symbol $I$.

In Appendix A, we show that the ZPB Eq. (II.2) can be expressed as a
generalized Poisson identity. The iterative solution of the ZPB Eq. can be
carried out by setting first the integral (16) to $I=1$, and updating the
value of $I$ at each iterative step until numerical convergence is achieved.
This solution scheme should be applied with the electrostatic boundary
conditions

$\phi^{\prime}(z_{\rm c}^{+})-\phi^{\prime}(z_{\rm c}^{-})=2/\mu_{-};\hskip
14.22636pt\phi^{\prime}(z_{\rm s})=0.$ (17)

In Eq. (17), $z_{\rm c}=\\{p,d-p\\}$ is the position of the fixed dipolar
anions. Moreover, the anionic and cationic Gouy-Chapman (GC) lengths are
defined as Gouy ; Chapman

$\mu_{\pm}=\left(2\pi Q_{\pm}\ell_{\rm B}\sigma_{\rm s}\right)^{-1}.$ (18)

Finally, the boundaries of the finite field region read $z_{\rm
s}=\\{p-a,d-p+a\\}$ for the salt-free liquid, and $z_{\rm s}=\\{0,d\\}$ in the
presence of added salt.

Figure 3: (Color online) Interaction pressure between the dipolar membrane
walls as a function of the membrane separation $d$ for varying (a) dipolar
cation valency, (b) salt concentration, and (c) dipolar anion valency. Solid
lines are from Eq. (II.3) with the numerical solution of the ZPB Eq. (II.2),
the open circles in (a) display the Debye-Hückle-level pressure (III.2), and
the red circles in the inset of (b) correspond to Eq. (39) for pure solvent.
The model parameters are $p=5$ Å, $a=3$ Å, and $\sigma_{\rm s}=0.1$
${\color[rgb]{0,0,0}\rm nm^{-2}}$.

###  II.3 Contact value theorem

Here, we calculate the inter-membrane disjoining pressure and derive the
contact value theorem. In the remainder of this work, we consider a symmetric
monovalent salt with ionic valencies $q_{\pm}=\pm 1$ and bulk concentrations
$\rho_{\pm{\rm b}}=\rho_{\rm b}$. The MF grand potential of the system
corresponds to the Hamiltonian (II.1) evaluated with the solution of the
saddle-point Eq. (II.2),

$\displaystyle\frac{\beta\Omega_{\rm mf}}{S}$ $\displaystyle=$
$\displaystyle-\int_{0}^{d}\frac{\mathrm{d}z}{8\pi\ell_{\rm
B}}\left[\phi^{\prime}(z)\right]^{2}-\int_{0}^{d}\mathrm{d}z\left[\rho_{+}(z)+\rho_{-}(z)\right]$
(19) $\displaystyle-\sigma_{\rm
s}Q_{-}\left[\phi(p)+\phi(d-p)\right]-\sigma_{\rm s}\ln\left[I_{\rm l}I_{\rm
r}\right].$

The MF pressure is obtained from the grand potential (19) via the identity
$P_{\rm in}=-\delta\Omega_{\rm mf}/\delta(Sd)$. Subtracting from the latter
the bulk pressure $\beta P_{\rm out}=2\rho_{\rm b}$, after some algebra, the
net pressure $P=P_{\rm in}-P_{\rm out}$ follows as

$\displaystyle\beta P$ $\displaystyle=$ $\displaystyle\frac{\sigma_{\rm
s}Q_{-}}{2}\left\\{\phi^{\prime}\left[(d-p)^{+}\right]+\phi^{\prime}\left[(d-p)^{-}\right]\right\\}$
$\displaystyle+\rho_{+}(d)+\rho_{-}(d)-2\rho_{\rm b}$
$\displaystyle+\rho^{+}_{\rm r}(d-p+a)-\rho^{+}_{\rm r}(d-p-a).$

Eq. (II.3) generalizes the contact value theorem, originally introduced for
charged membranes Hen1 ; Hen2 ; Isr , to membranes decorated with a dipolar
charge. The first line on the r.h.s. of Eq. (II.3) corresponds to the
electrostatic pressure, while the second line is the standard osmotic pressure
associated with the contact ion density at the membrane surface. Finally, the
third line is the osmotic pressure contribution stemming from the surface
dipoles. Interestingly, one sees that the dipolar pressure component involves
the density contrast of the surface dipoles between the orientational
configurations $\theta=\pi$ and $\theta=0$. The discrepancy between this form
and the standard ion contact density stems from the fact that the dipolar
cations, unable to touch the membrane surface, transmit to the latter the
mechanical stress originating from their rotational kinetic energy via their
tail.

In the standard PB framework the contact theorem and the first integral of the
PB equation are connected. We now analyze their relation in the present model.
In Appendix B, we derive the first integral of the ZPB Eq. (II.2) and obtain
the following more general expression

$\displaystyle\beta P$ $\displaystyle=$
$\displaystyle-\frac{\phi^{\prime}(z)^{2}}{8\pi\ell_{\rm
B}}+\sum_{n=1}^{s}\left[\rho_{n}(z)-\rho_{n{\rm b}}\right]+\rho^{+}_{\rm
r}(z)+\rho^{+}_{\rm l}(z)$ (21) $\displaystyle+\rho^{+}_{\rm
r}(d-p+a)\theta(z-d+p-a)$ $\displaystyle-\rho^{+}_{\rm
r}(d-p-a)\theta(z-d+p+a)$ $\displaystyle+\rho^{+}_{\rm
l}(p-a)\theta(p-a-z)-\rho^{+}_{\rm l}(p+a)\theta(p+a-z)$
$\displaystyle+\frac{\sigma_{\rm
s}Q_{-}}{2}\left\\{\phi^{\prime}\left[(d-p)^{+}\right]+\phi^{\prime}\left[(d-p)^{-}\right]\right\\}$
$\displaystyle\hskip 35.56593pt\times\theta(z-d+p)$
$\displaystyle-\frac{\sigma_{\rm
s}Q_{-}}{2}\left[\phi^{\prime}(p^{+})+\phi^{\prime}(p^{-})\right]\theta(p-z),$

for the disjoining pressure. In agreement with the mechanical equilibrium
constraint, Eq. (21) corresponding to the first integral, is of course
independent of the coordinate $z$.

While on first sight the first integral differs from the contact value
theorem, it is straightforward to verify that at the point $z=d$, Eq. (21)
yields consistently exactly the contact value pressure (II.3), derived
previously from the variation of the thermodynamic grand potential (19). Then,
at the mid pore, $z=d/2$, one obtains from Eq. (21) the following alternative
form of the pressure

$\displaystyle\beta P$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{s}\left[\rho_{n}(d/2)-\rho_{n{\rm b}}\right]$
$\displaystyle+\left[\rho_{{\rm l}+}(d/2)-\rho_{{\rm l}+}(p+a)\right.$
$\displaystyle\hskip 8.53581pt\left.+\rho_{{\rm r}+}(d/2)-\rho_{{\rm
r}+}(d-p-a)\right]\theta(p+a-d/2).$

One sees that at large separation distances $d>2(p+a)$ where the dipolar terms
in Eq. (II.3) vanish, the membrane interactions are purely due to the osmotic
pressure difference of the salt ions in the mid-pore and the bulk reservoir.
Then, in the close distance regime $d\leq 2(p+a)$ where the dipolar layers on
the opposing membranes overlap, the interaction pressure is augmented by the
dipolar density contrast between the mid-pore and the perpendicular dipole
configuration $\theta=0$. Below, we probe the consequences of this dipolar
interaction mechanism on the disjoining pressure of the zwitterionic
membranes.

##  III Results

Fig. 3(a) illustrates the effect of freely rotating surface dipoles on
membrane interactions. The solid curves show the separation dependence of the
disjoining pressure (II.3) computed from the numerical solution of the ZPB Eq.
(II.2) for different values of the dipolar cation valency $Q_{+}$. One sees
that in the standard case of purely monopolar charged membranes (black curve
for $Q_{+}=0$), the approach of the membranes leads to the monotonous rise of
the uniformly repulsive pressure ($d\downarrow P\uparrow$) according to the DH
law Isr

$\beta P\approx 2\pi\ell_{\rm B}\sigma_{\rm
s}^{2}Q_{-}^{2}\frac{\cosh^{2}\left(\kappa_{\rm
s}p\right)}{\sinh^{2}\left(\kappa_{\rm s}d/2\right)}$ (23)

(see the open circles). Then, as the dipolar cation valency rises towards the
zwitterionic dipole limit $Q_{+}=Q_{-}$, the short distance branch of the
pressure gradually drops and finally turns from repulsive to strongly
attractive. Hence, the interaction of neutral surface dipoles gives rise to
non-monotonic inter-membrane forces, characterized by a large-separation
repulsive regime and a short distance attractive branch, separated by a
pressure barrier. This membrane adhesion effect occurring on the level of MF
electrostatic conditions is the key prediction of our work.

Fig. 3(b) illustrates the influence of added salt on the interaction pressure
of a purely zwitterionic membrane with $Q_{\pm}=1$. One notes that the
increment of salt moderately weakens the short-range attractive interactions
but strongly amplifies the repulsive pressure branch. Indeed, one sees that in
the salt-free case (red curve at $\rho_{\rm b}=0$), the repulsive branch dies
out and the intermembrane interactions become purely attractive. This suggests
that in zwitterionic membranes, the inter-membrane repulsions occurring at
large separations are mediated by salt ions. In order to gain a better
physical insight into the features revealed so far, we consider below the
analytically tractable limits of the ZPB Eq. (II.2).

###  III.1 Explaining the membrane adhesion mechanism

In order to explain the mechanism driving the membrane attraction regime, we
consider here the salt-free case ($\rho_{n{\rm b}}=0$), where the ZPB Eq.
(II.2) possesses an exact analytical solution. The derivation of this solution
is detailed in Appendix D. As the interfacial dipole layers start to overlap
at $d\leq 2(p+a)$, the solution of Eq. (II.2) has to be considered in two
distinct regimes.

####  III.1.1 Decoupled dipole regime $d\geq 2(p+a)$

In the separation regime $d\geq 2(p+a)$ where the rotating dipoles on
different membrane surfaces do not overlap, the solution of Eq. (II.2) has the
piecewise form

$\displaystyle\phi(z)$ $\displaystyle=$
$\displaystyle\phi_{0}\theta(z)\theta(p-a-z)+\phi_{1}(z)\theta(z-p+a)\theta(p-z)$
$\displaystyle+\phi_{2}(z)\theta(z-p)\theta(p+a-z)$
$\displaystyle+\phi_{0}\theta(z-p-a)\theta(d/2-z)$
$\displaystyle+\phi(d-z)\theta(z-d/2)\theta(d-z),$

In Eq. (III.1.1), $\phi_{0}$ is a constant reference potential, the non-
uniform potential components

$\displaystyle\phi_{1}(z)=\phi_{0}+\frac{1}{Q_{+}}\ln\left\\{\cos^{2}\left[\frac{\kappa_{0}}{\sqrt{2}}(z-p+a)\right]\right\\},$
(25)
$\displaystyle\phi_{2}(z)=\phi_{0}+\frac{1}{Q_{+}}\ln\left\\{\cos^{2}\left[\frac{\kappa_{0}}{\sqrt{2}}(z-p-a)\right]\right\\}$
(26)

are similar in form to the monopolar membrane potential in counterion-only
liquids Isr , and the screening parameter $\kappa_{0}$ satisfies the
transcendental equation

$\frac{\kappa_{0}a}{\sqrt{2}}\tan\left(\frac{\kappa_{0}a}{\sqrt{2}}\right)=\frac{Q_{+}a}{2\mu_{-}}.$
(27)

Moreover, the fourth line of Eq. (III.1.1) indicates the mirror symmetry of
the potential solution with respect to the point $z=d/2$, i.e.
$\phi(d-z)=\phi(z)$. We finally note that due to the absence of dipolar
contact at $d>2(p+a)$, the electrostatic potential (III.1.1) does not depend
on the membrane separation distance $d$.

Figure 4: (Color online) (a) Electrostatic potential (III.1.1) and (b)
orientational probability (29) in the large distance regime $d\geq 2(p+a)$ of
the salt-free system ($\rho_{\rm b}=0$), where the surface dipoles are
electrostatically decoupled. (c)-(d) The plots in (a)-(b) in the short
separation regime $d<2(p+a)$, where the dipolar layers interact. The black and
blue arrows indicate the electric field lines $\mathbf{E}=-\nabla\phi(z)$ and
the resulting force $\mathbf{F}=Q_{+}\mathbf{E}$ on the dipolar cations
generated by the dipolar anions $-Q_{-}$ at $z=p$ and $z=d-p$, respectively.
The model parameters are $p=5$ Å, $a=3$ Å, $Q_{\pm}=1$, and $\sigma_{\rm
s}=1.0$ ${\rm nm}^{-2}$.

Fig. 4(a) displays the potential (III.1.1) corresponding to the dipolar cation
energy (see Eq.(13)) close to the wall at $z=0$. One sees that the electric
field lines $\mathbf{E}=-\nabla\phi(z)$ pointing towards their anionic charge
sources are oriented perpendicular to the membrane (black arrows). The
resulting force $\mathbf{F}=Q_{+}\mathbf{E}$ drives the dipolar cations
towards the potential minimum at $z=p$. The corresponding dipolar orientation
can be identified via the orientational probability density

$n(a_{\rm z})=\frac{1}{I}e^{-Q_{+}\phi(p+a_{z})}$ (28)

derived in Appendix C. Plugging Eq. (III.1.1) into Eq. (28), the probability
function follows as

$n(a_{\rm z})=\frac{\kappa a/\sqrt{2}}{\tan\left(\kappa
a/\sqrt{2}\right)}\sec^{2}\left[\frac{\kappa}{\sqrt{2}}\left(|a_{\rm
z}|-a\right)\right].$ (29)

Fig. 4(b) shows that the probability (29) is maximized at the angle
$\theta=\pi/2$. This indicates that in the separation regime $d>2(p+a)$, the
electrostatic force $\mathbf{F}$ allignes the dipoles parallel with the
membrane surface (see the inset).

In order to calculate the membrane interaction force corresponding to this
configuration, we substitute the potential solution (III.1.1) into Eq. (II.3).
This yields $P=0$, which confirms that the red pressure curve in Fig. 3(b)
vanishes at $d\geq 2(p+a)$. The cancellation of the interaction force at large
separation is due to the vanishing of the electrostatic field
$E_{z}(z)=-\partial_{z}\phi_{0}=0$ in the inter-membrane region $p+a<z<d-p-a$
separating the dipolar layers. This leads to the electrostatic decoupling of
the membranes, and the invariance of their individual electrostatic energy
under the variation of their separation.

####  III.1.2 Coupled dipole regime $2(p+a)\geq d\geq 2p+a$

For pore sizes located in the region $2(p+a)\geq d\geq 2p+a$ where the dipolar
regions overlap, the piecewise solution of Eq. (II.2) reads

$\displaystyle\phi(z)$ $\displaystyle=$
$\displaystyle\phi_{0}\theta(z)\theta(p-a-z)+\phi_{1}(z)\theta(z-p+a)\theta(p-z)$
(30) $\displaystyle+\phi_{2}(z)\theta(z-p)\theta(d-p-a-z)$
$\displaystyle+\phi_{3}(z)\theta(z-d+p+a)\theta(d/2-z)$
$\displaystyle+\phi(d-z)\theta(z-d/2)\theta(d-z).$

In Eq. (30), the potential components are

$\displaystyle\phi_{1}(z)=\phi_{0}+\frac{1}{Q_{+}}\ln\left\\{\cos^{2}\left[\frac{\kappa_{1}}{\sqrt{2}}(z-p+a)\right]\right\\},$
(31)
$\displaystyle\phi_{2}(z)=\phi_{0}+\frac{1}{Q_{+}}\ln\left\\{\frac{\kappa_{1}^{2}}{\kappa_{2}^{2}}\cos^{2}\left[\frac{\kappa_{2}}{\sqrt{2}}(z-z_{0})\right]\right\\},$
(32)
$\displaystyle\phi_{3}(z)=\phi_{0}+\frac{1}{Q_{+}}\ln\left\\{\frac{\kappa_{1}^{2}}{\kappa_{3}^{2}}\cos^{2}\left[\kappa_{3}\left(z-\frac{d}{2}\right)\right]\right\\},$
(33)

and the parameters $\kappa_{1}$, $\kappa_{2}$, $\kappa_{3}$, and $z_{0}$
satisfy the coupled transcendental Eqs.

$\displaystyle\cos\left[\frac{\kappa_{2}}{\sqrt{2}}(p-z_{0})\right]=\frac{\kappa_{2}}{\kappa_{1}}\cos\left(\frac{\kappa_{1}a}{\sqrt{2}}\right),$
(34)
$\displaystyle\tan\left[\frac{\kappa_{2}}{\sqrt{2}}(p-z_{0})\right]=\frac{\kappa_{1}}{\kappa_{2}}\tan\left(\frac{\kappa_{1}a}{\sqrt{2}}\right)-\frac{\sqrt{2}Q_{+}}{\kappa_{2}\mu_{-}},$
$\displaystyle\cos\left[\frac{\kappa_{2}}{\sqrt{2}}(d-p-a-
z_{0})\right]=\frac{\kappa_{2}}{\kappa_{3}}\cos\left[\kappa_{3}\left(\frac{d}{2}-p-a\right)\right],$
$\displaystyle\sin\left[\frac{\kappa_{2}}{\sqrt{2}}(d-p-a-
z_{0})\right]=\sqrt{2}\sin\left[\kappa_{3}\left(\frac{d}{2}-p-a\right)\right].$

In Appendix D.2.2, we show that by plugging the potential profile (30) into
Eq. (II.3), and using Eqs. (34)-(III.1.2), one obtains the interaction
pressure in the form

$\beta P=\frac{\kappa_{2}^{2}-\kappa_{3}^{2}}{2\pi\ell_{\rm B}Q_{+}^{2}}.$
(38)

Therein, by solving Eqs. (34)-(III.1.2) via a Taylor expansion in terms of the
dipole density $\sigma_{\rm s}$, we find that the pressure (38) takes the
approximate closed-form

$\beta P\approx-2\pi\ell_{\rm B}Q_{+}^{2}\sigma_{\rm
s}^{2}\frac{\left(p+a-d/2\right)^{2}}{a^{2}}$ (39)

Eq. (39) is displayed in the inset of Fig. 3(b) by red dots. This limiting law
indicates that as the membrane surfaces approach each other, the attractive
force is amplified algebraically with the width of the dipolar overlap region.
One also notes that as the dipolar docking occurs at $d=2p+a$, this attractive
pressure becomes independent of the dipolar size $a$, i.e. $\beta
P\approx-\pi\ell_{\rm B}Q_{+}^{2}\sigma_{\rm s}^{2}/2$.

In order to explain the origin of the attractive force in Eq. (39), we display
in Fig. 4(c) the electrostatic potential profile (30) together with the dipole
configuration close to the membrane surface at $z=0$. One sees that as the
dipolar layers overlap at $d<2(p+a)$, the attraction of the dipolar cations by
the electric field of the anions on the adjacent membrane wall comes into play
(blue arrows). This gives rise to a second potential minimum at $z=p+a$ (black
curve) that deepens with the decrease of the separation distance down to
$d=2p+a$ (red curve). To determine the consequence of this energy minimum on
the dipolar ordering, we substitute Eq. (30) into the orientational
probability density (28) to obtain

$\displaystyle n(a_{\rm z})$ $\displaystyle=$
$\displaystyle\frac{a\mu_{-}}{Q_{+}}\left\\{\kappa_{1}^{2}\sec^{2}\left[\kappa_{1}(a_{z}+a)/\sqrt{2}\right]\theta(-a_{z})\right.$
$\displaystyle\hskip
28.45274pt+\kappa_{2}^{2}\sec^{2}\left[\kappa_{2}(p+a_{z}-z_{0})/\sqrt{2}\right]$
$\displaystyle\hskip 36.98866pt\times\theta(a_{z})\theta(d-2p-a-a_{z})$
$\displaystyle\hskip
28.45274pt+\kappa_{3}^{2}\sec^{2}\left[\kappa_{3}(p+a_{z}-d/2)/\sqrt{2}\right]$
$\displaystyle\hskip 36.98866pt\left.\times\theta(a_{z}-d+2p+a)\right\\}.$

Eq. (III.1.2) plotted in Fig. 4(d) indicates that the dipolar cation
attraction towards the adjacent membrane wall translates into a second
probability peak in the perpendicular dipole orientation $\theta=0$ (black
curve). As the wall distance decreases further to $d=2p+a$ (red curve), this
probability peak reaches the height of the second maximum at the parallel
orientation $\theta=\pi/2$.

Hence, as the membrane surfaces approach each other, the amplification of the
dipolar attraction between the opposing walls is accompanied with the
reorientation of the dipoles from parallel to perpendicular configuration (see
the inset of Fig. 4(d)), akin to the formation of dipolar bridges across the
inter-membrane space. This dipolar bridging mechanism is precisely at the
origin of the attractive membrane interactions occurring at short separation
distances. It is noteworthy that such a short-ranged attractive energy
component has been identified in various MD simulations of dipolar membranes
sim1 ; sim2 . Next, we investigate the effect of added salt on this membrane
adhesion mechanism.

###  III.2 Understanding the salt-mediated repulsion

Having elucidated the mechanism governing the short-range attractive branch of
the pressure curves in Fig. 3(a), we investigate here the longe-range
repulsive interaction regime mediated by salt. To this aim, we focus on the DH
regime where the ZPB Eq. (II.2) can be linearized and solved exactly. The
derivation of this DH-level solution is presented in Appendix E. Injecting the
DH potential profile into Eq. (II.3), and Taylor-expanding the result at the
quadratic order in the dipole density $\sigma_{\rm s}$, the interaction
pressure follows for $d\geq 2p+a$ as

$\displaystyle\frac{\beta P}{2\pi\ell_{\rm B}\sigma_{\rm s}^{2}}$
$\displaystyle\approx$ $\displaystyle Q_{\rm
eff}^{2}\frac{\cosh^{2}\left(\kappa_{\rm
s}p\right)}{\sinh^{2}\left(\kappa_{\rm s}d/2\right)}$ $\displaystyle-
Q_{+}^{2}\frac{\sinh^{2}\left[\kappa_{\rm
s}\left(p+a-d/2\right)\right]}{\left(\kappa_{\rm
s}a\right)^{2}}\theta\left(p+a-d/2\right),$

where we introduced the salt screening parameter $\kappa_{\rm s}$ and the
effective dipole charge $Q_{\rm eff}$ defined as,

$\kappa_{\rm s}=\sqrt{8\pi\ell_{\rm B}\rho_{\rm b}};\hskip 28.45274ptQ_{\rm
eff}=Q_{+}\frac{\sinh\left(\kappa_{\rm s}a\right)}{\kappa_{\rm s}a}-Q_{-}.$
(42)

Fig. 3(a) shows that Eq. (III.2) can accurately approximate the numerical
pressure curves (open circles). We now note that Eq. (III.2) is composed of a
repulsive component (first term), and an attractive component (second term)
emerging in the short distance regime $d<2(p+a)$ studied in Sec. III.1.2.
Indeed, one sees that in salt-free liquids ($\kappa_{\rm s}=0$), the repulsive
component vanishes for $Q_{-}=Q_{+}$, and the attractive term tends to Eq.
(39).

In order to identify the origin of the repulsive term in Eq. (III.2), we note
that its magnitude is set by the effective charge $Q_{\rm eff}$ in Eq. (42).
Taylor-expanding the latter, and setting $Q_{-}=Q_{+}$, one obtains the
positive dipole charge $Q_{\rm eff}\approx Q_{+}\left(\kappa_{\rm
s}a\right)^{2}/6$ scaling linearly with the salt concentration. Thus, at large
distances, the interaction of zwitterionic membranes separated by a salt
solution is equivalent to the coupling of interfacially cationic walls not3 .

The behavior of the salt-dressed membrane dipoles as effective surface cations
stems from the screening of their elementary charges. Namely, in a salt-free
solution, and outside the region scanned by the rotating dipole, the field
created by the dipolar cation is exactly canceled by the dipolar anion of the
same magnitude. As a result, at $d>2(p+a)$, the intermembrane field and
interaction pressure vanish. However, in the presence of added salt, the field
of the screened dipolar anion can only partially cancel the field induced by
the rotating cation. This generates a net cationic inter-membrane field giving
rise to long-ranged repulsive membrane interactions.

Due to the competition between the attractive and repulsive interaction
components, the net pressure (III.2) vanishes at the characteristic separation
distance

$\displaystyle d^{*}$ $\displaystyle=$ $\displaystyle p+a+\kappa_{\rm
s}^{-1}{\rm arccosh}\left\\{\cosh\left[\kappa_{\rm s}(p+a)\right]\right.$
$\displaystyle\hskip 88.2037pt\left.-2\kappa_{\rm s}a\cosh(\kappa_{\rm
s}p)\frac{|Q_{\rm eff}|}{Q_{+}}\right\\}.$

Setting $Q_{-}=Q_{+}$ and Taylor-expanding Eq. (III.2) in terms of the
screening parameter, one obtains for $\kappa_{\rm s}a\lesssim 1$

$d^{*}\approx 2(p+a)-\frac{\kappa_{\rm s}a^{3}}{3(p+a)}.$ (44)

In agreement with Fig. 3(b), Eq. (44) indicates that added salt drops the
equilibrium separation ($\rho_{\rm b}\uparrow d^{*}\downarrow$) below the
close contact distance $d=2(p+a)$.

Figure 5: (Color online) (a) Cation and (b) anion density profiles from Eq.
(11) at various pore sizes in the separation regime $d\geq 2(p+a)$. (c)-(d)
Density plots in (a)-(b) in the dipolar contact regime $d\leq 2(p+a)$. The
parameters are $p=5$ Å, $a=3$ Å, $Q_{\pm}=1$, $\rho_{\rm b}=0.1$ M, and
$\sigma_{\rm s}=1.0$ ${\rm e}/{\rm nm}^{2}$.

###  III.3 Salt affinity of zwitterionic membranes and protonation effects on
membrane interactions

Figs. 5(a)-(b) illustrate the salt affinity of the zwitterionic membrane in
the decoupled dipole regime $d>2(p+a)$. One sees that due to the cooperative
effect of the anionic and cationic fields, the salt cations (anions) are
adsorbed (excluded) by the dipolar layers but excluded (adsorbed) by the
intermembrane region. Moreover, as the separation distance is decreased, the
amplification of the cationic intermembrane field enhances the cation
depletion ($d\downarrow\rho_{+}(z)\downarrow$) and the anion excess
($d\downarrow\rho_{-}(z)\uparrow$ ) in the mid-pore area.

Figs. 5(c)-(d) show that upon the further decrease of the separation distance
into the dipolar contact regime, the salt affinity of the intermembrane region
is entirely reversed. For $d\leq 2(p+a)$, the field of the dipolar cations is
dominated by the field of the dipolar anions on the opposing wall. The
resulting net anionic field responsible for the attractive dipolar bridging
mechanism in Fig. 4(d) switches the cation (anion) affinity of the mid-pore
area from depletion (excess) to excess (exclusion) state. This intricate
correlation between the slit size and the salt affinity of the dipolar
membrane may have significant impacts on nano-filtration and ion transport
techniques.

Previous electrophoretic transport experiments with zwitterionic lipid
vesicles revealed that the regulation of the solution pH via buffer addition
reduces the anionic surface charge strength of dipolar membranes Petrache2011
. The consequence of this protonation effect on membrane interactions is
qualitatively characterized in Fig. 3(c). One sees that upon the reduction of
$Q_{-}$, the attractive interaction regime is suppressed, and the net pressure
becomes more repulsive. The enhancement of the membrane repulsion by the
reduction of $Q_{-}$ stems from the increase of the effective dipolar charge
in Eq. (42), and the amplification of the repulsive force component in Eq.
(III.2), i.e. $Q_{-}\downarrow Q_{\rm eff}\uparrow P\uparrow$. Taylor-
expanding now Eq. (III.2) for $Q_{\rm eff}/Q_{+}\lesssim 1$, one obtains

$d^{*}\approx 2(p+a)-2a\frac{|Q_{\rm eff}|}{Q_{+}}\frac{\cosh(\kappa_{\rm
s}p)}{\sinh\left[\kappa_{\rm s}(p+a)\right]}.$ (45)

In agreement with Eq. (45), Fig. 3(c) indicates that the equilibrium distance
is expected to be reduced by protonation, i.e. $Q_{-}\downarrow Q_{\rm
eff}\uparrow d^{*}\downarrow$.

##  IV Conclusions

We developed a field theoretic model of electrolytes confined between dipolar
walls. Within the framework of this model, we investigated the electrostatic
interactions of zwitterionic membranes. We found that the intermembrane forces
are characterized by a non-uniform behavior governed by repulsive interactions
at large separation distances $d>2(p+a)$ and attractive interactions in the
short distance regime $d<2(p+a)$.

The repulsive interactions are mediated by the salt screening of the surface
dipoles that causes the opposing zwitterionic walls to behave as like-charged
cationic substrates. The attractive pressure branch is in turn triggered by a
dipolar bridging mechanism. Namely, as the membrane walls approach each other,
the dipolar cations attracted by the anions of the adjacent wall reorient
themselves towards the latter to minimize their distance. This leads to the
adhesion of the zwitterionic membrane walls via the electrostatic binding of
their surface dipoles.

We also observed that due to the competition between the individual field
components induced by the dipolar charges, the salt affinity of the
zwitterionic membrane, exhibits a complex behavior. In the large separation
regime where the field of the dipolar cations governs the inter-membrane
region, the latter zone is characterized by the salt anion excess and the salt
cation depletion. However, the drop of the membrane separation distance into
the dipolar overlap regime reverses the ion selectivity of the membrane; in
the inter-membrane region, the field of the dipolar anions taking over the
cationic field attracts the salt cations and excludes the salt anions. This
non-uniform trend of the membrane salt affinity may have significant
repercussions on water nanofiltration and nanofluidic transport techniques.

Our prediction of zwitterionic membrane attraction at short distances is
partly supported by previous MD simulations where the direct interaction of
the surface dipoles on the adjacent membrane walls has been indeed shown to be
attractive sim1 ; sim2 . However, the same works also revealed that this
short-range electrostatic attraction is overpowered by repulsive hydration
interactions mediated by restructuring of the solvent molecules. In order to
quantify the competition between the direct electrostatic interactions and
solvent-driven effects, we will generalize our model via the inclusion of
explicit solvent molecules together with beyond-MF charge correlations. The
inclusion of electrostatic fluctuations will also enable us to account for the
additional attractive effect of thermal van der Waals interactions. These
extensions will be considered in an upcoming article.

Acknowledgements RP acknowledges funding from the Key project #12034019 of the
National Natural Science Foundation of China.

Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Mapping Eq. (II.2) to the Poisson Eq.

Here, we show that Eq. (II.2) can be expressed as a generalized Poisson Eq.
Defining the dipolar on-site potentials as $V_{\rm
l,r}(\mathbf{r},\mathbf{\Omega})=V^{-}_{\rm l,r}(\mathbf{r})+V^{+}_{\rm
l,r}(\mathbf{r}+{\boldsymbol{a}})$, the MF-level number densities of the
dipolar anions $\rho^{-}_{\rm l,r}(\mathbf{r})$ and cations $\rho^{+}_{\rm
l,r}(\mathbf{r})$ can be obtained from the thermodynamic identity

$\rho^{\pm}_{\rm l,r}(\mathbf{r})=-\frac{\delta\ln Z_{\rm G}}{\delta
V^{\pm}_{\rm l,r}(\mathbf{r})}\approx\frac{\delta(\beta H)}{\delta
V^{\pm}_{\rm l,r}(\mathbf{r})}$ (46)

as

$\displaystyle\rho^{-}_{\rm l}(z)$ $\displaystyle=$ $\displaystyle\sigma_{\rm
s}\delta(z-p);\hskip 14.22636pt\rho^{-}_{\rm r}(z)=\sigma_{\rm
s}\delta(d-z+p),$ (47) $\displaystyle\rho^{+}_{\rm l}(z)$ $\displaystyle=$
$\displaystyle\frac{\sigma_{\rm s}}{2aI}e^{-Q_{+}\phi(z)}\theta_{\rm l}(z),$
(48) $\displaystyle\rho^{+}_{\rm r}(z)$ $\displaystyle=$
$\displaystyle\frac{\sigma_{\rm s}}{2aI}e^{-Q_{+}\phi(z)}\theta_{\rm r}(z).$
(49)

By identifying now the number densities (47)-(49) in the ZPB Eq. (II.2), one
can recast the latter as a Poisson Eq.,

$\frac{k_{\rm
B}T}{e^{2}}\partial_{z}\varepsilon(z)\partial_{z}\phi(z)+\sum_{n=1}^{s}q_{n}\rho_{n}(z)+\sum_{\alpha={\rm
l,r}}\sum_{\gamma=\pm}Q_{\gamma}\rho^{\gamma}_{\alpha}(z)=0.$ (50)

##  Appendix B Calculating the interaction pressure from the first integral of
the ZPB Eq. (II.2)

In this appendix, we calculate the first integral of the ZPB Eq. (II.2) and
derive a general form for the zwitterionic membrane pressure. To this aim, we
first multiply Eq. (II.2) with the local electric field $\phi^{\prime}(z)$.
This yields

$\displaystyle\phi^{\prime}(z)\phi^{\prime\prime}(z)+4\pi\ell_{\rm
B}\sum_{n=1}^{s}\rho_{{\rm b}n}q_{n}\phi^{\prime}(z)e^{-q_{n}\phi(z)}$
$\displaystyle-4\pi\ell_{\rm B}\sigma_{\rm
s}Q_{-}\phi^{\prime}(z)\left[\delta(z-p)+\delta(z-d+p)\right]$ (51)
$\displaystyle+\frac{2\pi\ell_{\rm B}}{aI}\sigma_{\rm
s}Q_{+}\phi^{\prime}(z)e^{-Q_{+}\phi(z)}\left[\theta_{\rm l}(z)+\theta_{\rm
r}(z)\right]=0.$

Next, we express Eq. (B) as a total derivative, i.e.

$\displaystyle-\frac{1}{8\pi\ell_{\rm B}}\frac{{\rm d}}{dz}\left(\frac{{\rm
d}\phi}{{\rm d}z}\right)^{2}+\frac{{\rm d}}{dz}\sum_{n=1}^{s}\rho_{{\rm
b}n}e^{-q_{n}\phi(z)}$ $\displaystyle+\sigma_{\rm s}Q_{-}\frac{{\rm
d}F(z)}{dz}-\frac{\sigma_{\rm s}}{2aI}\frac{{\rm d}G(z)}{dz}=0,$ (52)

where the functions $F(z)$ and $G(z)$ satisfy the Eqs.

$\displaystyle\frac{{\rm d}F(z)}{{\rm
d}z}=\phi^{\prime}(z)\left[\delta(z-p)+\delta(z-d+p)\right],$ (53)
$\displaystyle\frac{{\rm d}G(z)}{{\rm
d}z}=Q_{+}\phi^{\prime}(z)e^{-Q_{+}\phi(z)}\left[\theta_{\rm l}(z)+\theta_{\rm
r}(z)\right].$ (54)

Upon the formal integration of Eq. (B) from $z^{\prime}=0$ to $z^{\prime}=z$,
the first integral $C(z)$ of Eq. (B) becomes

$C(z)=-\frac{\phi^{\prime}(z)^{2}}{8\pi\ell_{\rm
B}}+\sum_{n=1}^{s}\left[\rho_{i}(z)-\rho_{i}(0)\right]+\sigma_{\rm
s}Q_{-}F(z)-\frac{\sigma_{\rm s}G(z)}{2aI},$ (55)

where the functions $F(z)$ and $G(z)$ follow from Eqs. (54)-(53) as

$\displaystyle F(z)=\int_{0}^{z}{\rm
d}z^{\prime}\phi^{\prime}(z^{\prime})\left[\delta(z^{\prime}-p)+\delta(z^{\prime}-d+p)\right],$
(56) $\displaystyle G(z)=Q_{+}\int_{0}^{z}{\rm
d}z^{\prime}\phi^{\prime}(z^{\prime})e^{-Q_{+}\phi(z^{\prime})}\left[\theta_{\rm
l}(z^{\prime})+\theta_{\rm r}(z^{\prime})\right].$ (57)

Taking into account the discontinuity of the field $\phi^{\prime}(z)$ at $z=p$
and $z=d-p$, the evaluation of the integral in Eq. (56) yields

$\displaystyle F(z)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\phi^{\prime}(p^{+})+\phi^{\prime}(p^{-})\right]\theta(z-p)$
$\displaystyle+\frac{1}{2}\left\\{\phi^{\prime}\left[(d-p)^{+}\right]+\phi^{\prime}\left[(d-p)^{-}\right]\right\\}\theta(z-d+p).$

Then, in order to obtain the function $G(z)$ from the integral (57), with the
theta functions defined in Eqs. (14)-(15), we carry out an integration by
parts, and use the identity $\partial_{x}\theta(x)=\delta(x)$. After some
algebra, one obtains

$\displaystyle G(z)$ $\displaystyle=$
$\displaystyle-e^{-Q_{+}\phi(z)}\left[\theta_{\rm l}(z)+\theta_{\rm
r}(z)\right]$ $\displaystyle+e^{-Q_{+}\phi(p-a)}\theta(z-p+a)$
$\displaystyle-e^{-Q_{+}\phi(p+a)}\theta(z-p-a)$
$\displaystyle+e^{-Q_{+}\phi(d-p-a)}\theta(z-d+p+a)$
$\displaystyle-e^{-Q_{+}\phi(d-p+a)}\theta(z-d+p-a).$

Finally, plugging Eqs. (B)-(B) into Eq. (55), the first integral of Eq. (II.2)
takes the form

$\displaystyle C(z)$ $\displaystyle=$
$\displaystyle-\frac{\phi^{\prime}(z)^{2}}{8\pi\ell_{\rm
B}}+\sum_{n=1}^{s}\left[\rho_{n}(z)-\rho_{n}(0)\right]$
$\displaystyle+\frac{\sigma_{\rm s}}{2aI}e^{-Q_{+}\phi(z)}\left[\theta_{\rm
l}(z)+\theta_{\rm r}(z)\right]$ $\displaystyle-\frac{\sigma_{\rm
s}}{2aI}\left\\{e^{-Q_{+}\phi(p-a)}\theta(z-p+a)\right.$ $\displaystyle\hskip
34.14322pt-e^{-Q_{+}\phi(p+a)}\theta(z-p-a)$ $\displaystyle\hskip
34.14322pt+e^{-Q_{+}\phi(d-p-a)}\theta(z-d+p+a)$ $\displaystyle\hskip
34.14322pt\left.-e^{-Q_{+}\phi(d-p+a)}\theta(z-d+p-a)\right\\}$
$\displaystyle+\frac{\sigma_{\rm
s}Q_{-}}{2}\left[\phi^{\prime}(p^{+})+\phi^{\prime}(p^{-})\right]\theta(z-p)$
$\displaystyle+\frac{\sigma_{\rm
s}Q_{-}}{2}\left\\{\phi^{\prime}\left[(d-p)^{+}\right]+\phi^{\prime}\left[(d-p)^{-}\right]\right\\}$
$\displaystyle\hskip 35.56593pt\times\theta(z-d+p).$

In order to derive the net pressure from the first integral (B), we will
normalize the latter with its value in the infinite membrane separation
regime. To this aim, we set $z=d/2$ and take the limit $d\to\infty$. This
yields

$\displaystyle C_{\rm b}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{s}\left[\rho_{n{\rm
b}}-\rho_{n}(0)\right]+\frac{\sigma_{\rm
s}Q_{-}}{2}\left[\phi^{\prime}(p^{+})+\phi^{\prime}(p^{-})\right]$ (61)
$\displaystyle+\frac{\sigma_{\rm
s}}{2aI}\left\\{e^{-Q_{+}\phi(p+a)}-e^{-Q_{+}\phi(p-a)}\right\\}.$

Subtracting from Eq. (B) the constant (61), and using the dipolar density
functions (13) together with the identity
$\theta(z-z^{\prime})-1=-\theta(z^{\prime}-z)$, the normalized first integral
corresponding to the net interplate pressure $\beta P=C(z)-C_{\rm b}$ follows
as

$\displaystyle\beta P$ $\displaystyle=$
$\displaystyle-\frac{\phi^{\prime}(z)^{2}}{8\pi\ell_{\rm
B}}+\sum_{n=1}^{s}\left[\rho_{n}(z)-\rho_{n{\rm b}}\right]+\rho^{+}_{\rm
r}(z)+\rho^{+}_{\rm l}(z)$ (62) $\displaystyle+\rho^{+}_{\rm
r}(d-p+a)\theta(z-d+p-a)$ $\displaystyle-\rho^{+}_{\rm
r}(d-p-a)\theta(z-d+p+a)$ $\displaystyle+\rho^{+}_{\rm
l}(p-a)\theta(p-a-z)-\rho^{+}_{\rm l}(p+a)\theta(p+a-z)$
$\displaystyle+\frac{\sigma_{\rm
s}Q_{-}}{2}\left\\{\phi^{\prime}\left[(d-p)^{+}\right]+\phi^{\prime}\left[(d-p)^{-}\right]\right\\}$
$\displaystyle\hskip 35.56593pt\times\theta(z-d+p)$
$\displaystyle-\frac{\sigma_{\rm
s}Q_{-}}{2}\left[\phi^{\prime}(p^{+})+\phi^{\prime}(p^{-})\right]\theta(p-z).$

We also note that as the mirror symmetry $\phi(d-z)=\phi(z)$ implies
$\rho^{+}_{\rm l}(p\pm a)=\rho^{+}_{\rm r}(d-p\mp a)$ and
$\phi^{\prime}(p^{\pm})=\mp\phi^{\prime}\left[(d-p)^{\mp}\right]$, the
pressure (62) can be further simplified as

$\displaystyle\beta P$ $\displaystyle=$
$\displaystyle-\frac{\phi^{\prime}(z)^{2}}{8\pi\ell_{\rm
B}}+\sum_{n=1}^{s}\left[\rho_{n}(z)-\rho_{n{\rm b}}\right]+\rho^{+}_{\rm
r}(z)+\rho^{+}_{\rm l}(z)$ $\displaystyle+\rho^{+}_{\rm
r}(d-p+a)\left[\theta(z-d+p-a)+\theta(p-a-z)\right]$
$\displaystyle-\rho^{+}_{\rm
r}(d-p-a)\left[\theta(z-d+p+a)+\theta(p+a-z)\right]$
$\displaystyle+\frac{\sigma_{\rm
s}Q_{-}}{2}\left\\{\phi^{\prime}\left[(d-p)^{+}\right]+\phi^{\prime}\left[(d-p)^{-}\right]\right\\}$
$\displaystyle\hskip 35.56593pt\times\left[\theta(z-d+p)+\theta(p-z)\right].$

In agreement with the mechanical equilibrium condition, the r.h.s. of the
first integral (62) is naturally independent of the coordinate $z$. In Sec.
II.3, we use Eq. (62) to relate the interaction pressure to the charge density
values in different zones of the intermembrane region.

##  Appendix C Deriving the orientational probability (28)

We derive here the orientational probability characterizing the rotational
configuration of the surface dipoles. To this aim, we redefine the dipolar on-
site potentials as $V_{\rm l,r}(\mathbf{r},\mathbf{\Omega})\to V_{\rm
l,r}(\mathbf{r})$ such that they act only on the dipolar anions fixed at
$\mathbf{r}$. This allows to obtain the fixed anion number densities as

$\rho^{-}_{\alpha}(\mathbf{r})=\left.\frac{\delta(\beta H)}{\delta
V^{\pm}_{\alpha}(\mathbf{r})}\right|_{V_{\alpha}=0}=:\sigma_{\rm
s}\delta(z-z_{\rm
c})\int\frac{\mathrm{d}^{2}\mathbf{\Omega}}{4\pi}n_{\alpha}(\mathbf{\Omega}),$
(64)

with the auxiliary parameter $z_{\rm c}=p$ if $\alpha=l$, and $z_{\rm c}=d-p$
for $\alpha=r$. The second equality of Eq. (64) defines the orientational
probability such that the latter yields the correct dipolar anion density
$\rho^{-}_{\alpha}(z)=\sigma_{\rm s}\delta(z-z_{\rm c})$ and normalization
condition
$\int\mathrm{d}^{2}\mathbf{\Omega}\;n_{\alpha}(\mathbf{\Omega})/(4\pi)=1$.
Evaluating the derivative in Eq. (64), and accounting for the symmetry with
respect to the mid-pore, the orientational probabilities follow as $n_{\rm
l}(a_{\rm z})=n_{\rm r}(a_{\rm z})=n(a_{\rm z})$, with

$n(a_{\rm z})=\frac{1}{I}e^{-Q_{+}\phi(p+a_{z})}.$ (65)

##  Appendix D Salt-free solution of the ZPB Eq. (II.2)

In this Appendix, we derive the exact solution of the ZPB Eq. (II.2) in a
salt-free solution. Setting $\rho_{n{\rm b}}=0$, Eq. (II.2) becomes

$\displaystyle\frac{k_{\rm
B}T}{e^{2}}\partial_{z}\varepsilon(z)\partial_{z}\phi(z)-\sigma_{\rm
s}Q_{-}\left[\delta(z-p)+\delta(z-d+p)\right]$
$\displaystyle+\frac{\sigma_{\rm
s}Q_{+}}{2aI}\left\\{e^{-Q_{+}\phi(z)}\theta_{\rm
l}(z)+e^{-Q_{+}\phi(z)}\theta_{\rm r}(z)\right\\}=0.$ (66)

We now note that according to Eqs. (14)-(15), the rotating dipole densities in
the curly bracket term of Eq. (D) overlap for $d\leq 2(p+a)$. This implies
that Eq. (D) should be solved in two different distance regimes.

###  D.1 Derivation of the solution (III.1.1) for $d\geq 2(p+a)$

We first consider the membrane separation regime $d>2(p+a)$ where the surface
dipoles on different walls do not overlap. Due to the mirror symmetry with
respect to the point $z=d/2$, and the theta functions (14)-(15) fixing the
regions scanned by the rotating dipoles, the solution to Eq. (D) has the
piecewise form

$\displaystyle\phi(z)$ $\displaystyle=$
$\displaystyle\phi_{0}\theta(z)\theta(p-a-z)+\phi_{1}(z)\theta(z-p+a)\theta(p-z)$
$\displaystyle+\phi_{2}(z)\theta(z-p)\theta(p+a-z)$
$\displaystyle+\phi_{3}\theta(z-p-a)\theta(d/2-z)$
$\displaystyle+\phi(d-z)\theta(z-d/2)\theta(d-z).$

In Eq. (D.1), $\phi_{0}$ and $\phi_{3}$ are constant potentials in the
unscreened slit zones inaccessible to the dipolar cations. The remaining
potential components with index $i=1,2$ solve the differential equation

$\frac{d^{2}\phi_{i}(z)}{dz^{2}}+\frac{1}{a\mu_{+}I}e^{-Q_{+}\phi_{i}(z)}=0.$
(68)

Multiplying Eq. (68) by $d\phi_{i}(z)/dz$, integrating the result, and
imposing the continuity of the potential and the field at $z=p\pm a$, the
first integrals of Eq. (68) follow as

$\displaystyle\frac{d\phi_{1}(z)}{dz}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}\kappa_{0}}{Q_{+}}\left\\{e^{-Q_{+}\left[\phi_{1}(z)-\phi_{0}\right]}-1\right\\}^{1/2},$
(69) $\displaystyle\frac{d\phi_{2}(z)}{dz}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}\kappa^{\prime}_{0}}{Q_{+}}\left\\{e^{-Q_{+}\left[\phi_{2}(z)-\phi_{3}\right]}-1\right\\}^{1/2},$
(70)

where $\kappa_{0}$ and $\kappa^{\prime}_{0}$ are integration constants.
Solving Eqs. (69)-(70), and imposing at $z=p$ the continuity of the potential
together with the Gauss’ law in Eq. (17), one finally obtains
$\phi_{3}=\phi_{0}$, $\kappa^{\prime}_{0}=\kappa_{0}$ with

$\frac{\kappa_{0}a}{\sqrt{2}}\tan\left(\frac{\kappa_{0}a}{\sqrt{2}}\right)=\frac{Q_{+}a}{2\mu_{-}},$
(71)

and the non-uniform potential components in Eq. (D.1),

$\displaystyle\phi_{1}(z)=\phi_{0}+\frac{1}{Q_{+}}\ln\left\\{\cos^{2}\left[\frac{\kappa_{0}}{\sqrt{2}}(z-p+a)\right]\right\\},$
(72)
$\displaystyle\phi_{2}(z)=\phi_{0}+\frac{1}{Q_{+}}\ln\left\\{\cos^{2}\left[\frac{\kappa_{0}}{\sqrt{2}}(z-p-a)\right]\right\\}.$
(73)

###  D.2 Derivation of the solution (30) for $2(p+a)\geq d\geq 2p+a$

####  D.2.1 Deriving the piecewise potential

In the slit size regime $2(p+a)\geq d\geq 2p+a$, the rotating dipoles on the
left and right walls overlap at $d-p-a\leq z\leq p+a$. As a result, the mid-
pore region is characterized by a non-uniform potential. Thus, the piecewise
solution to Eq. (D) reads

$\displaystyle\phi(z)$ $\displaystyle=$
$\displaystyle\phi_{0}\theta(z)\theta(p-a-z)+\phi_{1}(z)\theta(z-p+a)\theta(p-z)$
(74) $\displaystyle+\phi_{2}(z)\theta(z-p)\theta(d-p-a-z)$
$\displaystyle+\phi_{3}(z)\theta(z-d+p+a)\theta(d/2-z)$
$\displaystyle+\phi(d-z)\theta(z-d/2)\theta(d-z).$

In Eq. (74), $\phi_{0}$ is the uniform potential component in the pure solvent
region. Then, the non-uniform potential components solve the differential Eqs.

$\frac{d^{2}\phi_{i}(z)}{dz^{2}}+\frac{1}{a\mu_{+}I}e^{-Q_{+}\phi_{i}(z)}=0$
(75)

for $i=1,2$, and

$\frac{d^{2}\phi_{3}(z)}{dz^{2}}+\frac{2}{a\mu_{+}I}e^{-Q_{+}\phi_{3}(z)}=0$
(76)

in the dipolar overlap region.

Integrating Eqs. (75)-(76), and imposing the continuity of the solution
through the slit, the first integrals follow as

$\displaystyle\frac{d\phi_{1}(z)}{dz}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}\kappa_{1}}{Q_{+}}\left\\{e^{-Q_{+}\left[\phi_{1}(z)-\phi_{0}\right]}-1\right\\}^{1/2},$
(77) $\displaystyle\frac{d\phi_{2}(z)}{dz}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}\kappa_{2}}{Q_{+}}\left\\{\frac{\kappa_{1}^{2}}{\kappa_{2}^{2}}e^{-Q_{+}\left[\phi_{2}(z)-\phi_{0}\right]}-1\right\\}^{1/2},$
(78) $\displaystyle\frac{d\phi_{3}(z)}{dz}$ $\displaystyle=$
$\displaystyle\frac{2\kappa_{3}}{Q_{+}}\left\\{\frac{\kappa_{1}^{2}}{\kappa_{3}^{2}}e^{-Q_{+}\left[\phi_{3}(z)-\phi_{0}\right]}-1\right\\}^{1/2},$
(79)

where $\kappa_{i}$ with $i=\\{1,2,3\\}$ are integration constants. Solving
Eqs. (77)-(79), and imposing the continuity of the field and potential
together with the Gauss’ law in Eq. (17), the potential components finally
follow as

$\displaystyle\phi_{1}(z)=\phi_{0}+\frac{1}{Q_{+}}\ln\left\\{\cos^{2}\left[\frac{\kappa_{1}}{\sqrt{2}}(z-p+a)\right]\right\\},$
(80)
$\displaystyle\phi_{2}(z)=\phi_{0}+\frac{1}{Q_{+}}\ln\left\\{\frac{\kappa_{1}^{2}}{\kappa_{2}^{2}}\cos^{2}\left[\frac{\kappa_{2}}{\sqrt{2}}(z-z_{0})\right]\right\\},$
(81)
$\displaystyle\phi_{3}(z)=\phi_{0}+\frac{1}{Q_{+}}\ln\left\\{\frac{\kappa_{1}^{2}}{\kappa_{3}^{2}}\cos^{2}\left[\kappa_{3}\left(z-\frac{d}{2}\right)\right]\right\\},$
(82)

where the constants $\kappa_{i}$ and $z_{0}$ satisfy the non-linear Eqs.

$\displaystyle\cos\left[\frac{\kappa_{2}}{\sqrt{2}}(p-z_{0})\right]=\frac{\kappa_{2}}{\kappa_{1}}\cos\left(\frac{\kappa_{1}a}{\sqrt{2}}\right),$
(83)
$\displaystyle\tan\left[\frac{\kappa_{2}}{\sqrt{2}}(p-z_{0})\right]=\frac{\kappa_{1}}{\kappa_{2}}\tan\left(\frac{\kappa_{1}a}{\sqrt{2}}\right)-\frac{\sqrt{2}Q_{+}}{\kappa_{2}\mu_{-}},$
(84) $\displaystyle\cos\left[\frac{\kappa_{2}}{\sqrt{2}}(d-p-a-
z_{0})\right]=\frac{\kappa_{2}}{\kappa_{3}}\cos\left[\kappa_{3}\left(\frac{d}{2}-p-a\right)\right],$
(85) $\displaystyle\sin\left[\frac{\kappa_{2}}{\sqrt{2}}(d-p-a-
z_{0})\right]=\sqrt{2}\sin\left[\kappa_{3}\left(\frac{d}{2}-p-a\right)\right].$

We finally note that if one evaluates the integral $I$ in Eq. (16) with the
salt-free potential solutions (D.1) and (74) derived above, one obtains the
self-consistency condition $Q_{-}=Q_{+}$. The latter equality corresponds to
the interfacial electroneutrality constraint originating from the absence of
salt.

####  D.2.2 Deriving the analytical Eqs. (38) and (39) for pressure

Here, we derive closed form expressions for the intermembrane pressure in the
salt-free liquid. To this aim, we substitute the potential profile (74) into
Eq. (II.3). This yields

$\displaystyle\beta P$ $\displaystyle=$ $\displaystyle\frac{\mu_{+}\sigma_{\rm
s}\kappa_{1}^{2}}{2Q_{+}}\left\\{1-\frac{\kappa_{3}^{2}}{\kappa_{1}^{2}}\sec^{2}\left[\kappa_{3}\left(\frac{d}{2}-p-a\right)\right]\right\\}$
(87) $\displaystyle+\sigma_{\rm
s}\left\\{\sqrt{2}\kappa_{1}\tan\left(\frac{\kappa_{1}a}{\sqrt{2}}\right)-\frac{Q_{+}}{\mu_{-}}\right\\}.$

Next, we combine Eqs. (D.2.1) and (D.2.1) to get

$\kappa_{3}^{2}\sec^{2}\left[\kappa_{3}\left(\frac{d}{2}-p-a\right)\right]=2\kappa_{3}^{2}-\kappa_{2}^{2}.$
(88)

Then, combining Eqs. (83) and (84), one obtains

$\sqrt{2}\kappa_{1}\tan\left(\frac{\kappa_{1}a}{\sqrt{2}}\right)-\frac{Q_{+}}{\mu_{-}}=\frac{\kappa_{2}^{2}\mu_{-}}{\sqrt{2}Q_{+}}\left(1-\frac{\kappa_{1}^{2}}{\kappa_{2}^{2}}\right).$
(89)

Plugging the identities (88)-(89) into Eq. (87), the interplate pressure
finally simplifies to

$\beta P=\frac{\kappa_{2}^{2}-\kappa_{3}^{2}}{2\pi\ell_{\rm B}Q_{+}^{2}}.$
(90)

In order to obtain an insightful expression for the pressure (90), we carry-
out a perturbative expansion of the transcendental Eqs. (83)-(D.2.1) in terms
of the dipole density $\sigma_{\rm s}$. To this aim, we inject into these
equations the perturbative series $z_{0}\approx z_{01}+z_{02}\sqrt{\sigma_{\rm
s}}+z_{03}\sigma_{\rm s}+z_{04}\sigma_{\rm s}^{3/2}$, and
$\kappa_{i}\approx\kappa_{i1}\sqrt{\sigma_{\rm s}}+\kappa_{i2}\sigma_{\rm
s}+\kappa_{i3}\sigma_{\rm s}^{3/2}$ for $1\leq i\leq 3$. This yields

$\displaystyle z_{0}\approx p+a,$ (91)
$\displaystyle\kappa_{1}\approx\kappa_{2}\approx\sqrt{\frac{Q_{+}}{a\mu_{-}}}-\frac{a^{2}}{12}\left(\frac{Q_{+}}{a\mu_{-}}\right)^{3/2},$
$\displaystyle\kappa_{3}\approx\sqrt{\frac{Q_{+}}{a\mu_{-}}}+\frac{1}{24}\left(\frac{Q_{+}}{a\mu_{-}}\right)^{3/2}\left[3(2p+2a-d)^{2}-2a^{2}\right].$

Plugging the approximative solutions (D.2.2)-(D.2.2) into Eq. (90), the
dipolar membrane pressure follows as

$\beta P\approx-2\pi\ell_{\rm B}Q_{+}^{2}\sigma_{\rm
s}^{2}\frac{\left(p+a-d/2\right)^{2}}{a^{2}}.$ (94)

##  Appendix E DH-level solution of the ZPB Eq. (II.2) with added salt

In this appendix, we derive the DH-level analytical solution of the ZPB Eq.
(II.2) for a $1:1$ salt solution. Thus, we set $q_{\pm}=\pm 1$ and $\rho_{{\rm
b}n}=\rho_{\rm b}$, and linearize Eq. (II.2) in terms of the potential
$\phi(z)$. This yields

$\displaystyle\partial_{z}^{2}\phi(z)-\left\\{\kappa_{\rm
s}^{2}+(\xi_{1}^{2}-\kappa_{\rm s}^{2})\left[\theta_{\rm l}(z)+\theta_{\rm
r}(z)\right]\right\\}\phi(z)$ (95)
$\displaystyle=\frac{2}{\mu_{-}}\left[\delta(z-p)+\delta(z-d+p)\right]-\psi\left[\theta_{\rm
l}(z)+\theta_{\rm r}(z)\right],$

where we defined the auxiliary parameters

$\kappa_{\rm s}^{2}=8\pi\ell_{\rm B}\rho_{\rm b};\hskip
14.22636pt\xi_{1}^{2}=\kappa_{\rm s}^{2}+Q_{+}\psi;\hskip
14.22636pt\psi=\frac{1}{\mu_{+}aI}.$ (96)

###  E.1 Solution of Eq. (95) for $d\geq 2(p+a)$

In the separation regime $d\geq 2(p+a)$, the potential solution to Eq. (95)
follows in the piecewise form

$\displaystyle\phi(z)$ $\displaystyle=$
$\displaystyle\phi_{1}(z)\theta(z)\theta(p-a-z)+\phi_{2}(z)\theta(z-p+a)\theta(p-z)$
$\displaystyle+\phi_{3}(z)\theta(z-p)\theta(p+a-z)$
$\displaystyle+\phi_{4}(z)\theta(z-p-a)\theta(d/2-z)$
$\displaystyle+\phi(d-z)\theta(z-d/2)\theta(d-z),$

with the potential components

$\displaystyle\phi_{1}(z)$ $\displaystyle=$ $\displaystyle
b_{1}\cosh(\kappa_{\rm s}z),$ (98) $\displaystyle\phi_{2}(z)$ $\displaystyle=$
$\displaystyle b_{2}e^{\xi_{1}z}+b_{3}e^{-\xi_{1}z}+\frac{\psi}{\xi_{1}^{2}},$
(99) $\displaystyle\phi_{3}(z)$ $\displaystyle=$ $\displaystyle
b_{4}e^{\xi_{1}z}+b_{5}e^{-\xi_{1}z}+\frac{\psi}{\xi_{1}^{2}},$ (100)
$\displaystyle\phi_{4}(z)$ $\displaystyle=$ $\displaystyle
b_{6}\cosh\left[\kappa_{\rm s}\left(\frac{d}{2}-z\right)\right].$ (101)

Imposing the continuity of the potential (E.1) throughout the intermembrane
region together with the Gauss’ law (17), the coefficients in Eqs. (99)-(101)
follow as

$\displaystyle b_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{\gamma_{1}\sinh\left[\kappa_{\rm
s}(p-a)\right]}\left\\{(1-\Delta_{1})e^{\xi_{1}(p-a)}b_{2}+\frac{\Delta_{0}\psi}{\xi_{1}^{2}}\right\\},$
$\displaystyle b_{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{s_{-}}\frac{1+\Delta_{3}e^{-2\xi_{1}a}}{1-\Delta_{1}\Delta_{3}e^{-4\xi_{1}a}}e^{-\xi_{1}p}$
$\displaystyle-\frac{\psi}{\xi_{1}^{2}}\frac{\Delta_{2}+\Delta_{0}\Delta_{3}e^{-2\xi_{1}a}}{1-\Delta_{1}\Delta_{3}e^{-4\xi_{1}a}}e^{-\xi_{1}(p+a)},$
$\displaystyle b_{3}$ $\displaystyle=$ $\displaystyle
b_{2}\Delta_{1}e^{2\xi_{1}(p-a)}-\frac{\Delta_{0}\psi}{\xi_{1}^{2}}e^{\xi_{1}(p-a)};\;b_{4}=b_{2}+\frac{e^{-\xi_{1}p}}{s_{-}},$
$\displaystyle b_{5}$ $\displaystyle=$ $\displaystyle
b_{2}\Delta_{1}e^{2\xi_{1}(p-a)}-\frac{\Delta_{0}\psi}{\xi_{1}^{2}}e^{\xi_{1}(p-a)}-\frac{e^{\xi_{1}p}}{s_{-}},$
(105) $\displaystyle b_{6}$ $\displaystyle=$
$\displaystyle-\frac{1}{\gamma_{1}\sinh\left[\kappa_{\rm
s}\left(\frac{d}{2}-p-a\right)\right]}$ $\displaystyle\hskip
8.53581pt\times\left\\{b_{2}\left(1-\Delta_{1}e^{-4\xi_{1}a}\right)e^{\xi_{1}(p+a)}\right.$
$\displaystyle\hskip
25.60747pt\left.+\frac{2}{s_{-}}\cosh(\xi_{1}a)+\frac{\Delta_{0}\psi}{\xi_{1}^{2}}e^{-2\xi_{1}a}\right\\},$

where we defined the auxiliary parameters $s_{-}=\xi_{1}\mu_{-}$,
$\gamma_{1}=\kappa_{\rm s}/\xi_{1}$, and

$\displaystyle\Delta_{0}$ $\displaystyle=$
$\displaystyle\frac{\gamma_{1}\tanh\left[\kappa_{\rm
s}(p-a)\right]}{1+\gamma_{1}\tanh\left[\kappa_{\rm
s}(p-a)\right]};\;\Delta_{1}=1-2\Delta_{0},$ (107) $\displaystyle\Delta_{2}$
$\displaystyle=$ $\displaystyle\frac{\gamma_{1}\tanh\left[\kappa_{\rm
s}\left(\frac{d}{2}-p-a\right)\right]}{1+\gamma_{1}\tanh\left[\kappa_{\rm
s}\left(\frac{d}{2}-p-a\right)\right]},$ (108) $\displaystyle\Delta_{3}$
$\displaystyle=$ $\displaystyle 1-2\Delta_{2}.$ (109)

###  E.2 Solution of Eq. (95) for $2(p+a)\geq d\geq 2p+a$

In the distance regime $2(p+a)\geq d\geq 2p+a$, the piecewise potential
solving Eq. (95) has the form

$\displaystyle\phi(z)$ $\displaystyle=$
$\displaystyle\phi_{1}(z)\theta(z)\theta(p-a-z)$
$\displaystyle+\phi_{2}(z)\theta(z-p+a)\theta(p-z)$
$\displaystyle+\phi_{3}(z)\theta(z-p)\theta(d-p-a-z)$
$\displaystyle+\phi_{4}(z)\theta(z-d+p+a)\theta(d/2-z)$
$\displaystyle+\phi(d-z)\theta(z-d/2)\theta(d-z),$

with the potential components

$\displaystyle\phi_{1}(z)$ $\displaystyle=$ $\displaystyle
c_{1}\cosh(\kappa_{\rm s}z),$ (111) $\displaystyle\phi_{2}(z)$
$\displaystyle=$ $\displaystyle
c_{2}e^{\xi_{1}z}+c_{3}e^{-\xi_{1}z}+\frac{\psi}{\xi_{1}^{2}},$ (112)
$\displaystyle\phi_{3}(z)$ $\displaystyle=$ $\displaystyle
c_{4}e^{\xi_{1}z}+c_{5}e^{-\xi_{1}z}+\frac{\psi}{\xi_{1}^{2}},$ (113)
$\displaystyle\phi_{4}(z)$ $\displaystyle=$ $\displaystyle
c_{6}\cosh\left[\xi_{2}\left(\frac{d}{2}-z\right)\right]+\frac{2\psi}{\xi_{2}^{2}}.$
(114)

In Eqs. (111)-(114), we defined the additional parameters
$\xi_{2}^{2}=\kappa_{\rm s}^{2}+2Q_{+}\psi$ and $\gamma_{2}=\xi_{2}/\xi_{1}$,
and the coefficients

$\displaystyle c_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{\gamma_{1}\sinh\left[\kappa_{\rm
s}(p-a)\right]}\left\\{(1-\Delta_{1})e^{\xi_{1}(p-a)}c_{2}+\frac{\Delta_{0}\psi}{\xi_{1}^{2}}\right\\},$
$\displaystyle c_{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{s_{-}}\frac{e^{\xi_{1}(p+2a)}+\Delta_{5}e^{\xi_{1}(2d-3p)}}{\Delta_{5}e^{2\xi_{1}(d-p)}-\Delta_{1}e^{2\xi_{1}p}}$
$\displaystyle+\frac{\psi}{\Delta_{5}e^{2\xi_{1}(d-p)}-\Delta_{1}e^{2\xi_{1}p}}$
$\displaystyle\hskip
8.53581pt\times\left\\{\frac{1}{\xi_{1}^{2}}\left[\Delta_{4}e^{\xi_{1}(d-p+a)}-\Delta_{0}e^{\xi_{1}(p+a)}\right]\right.$
$\displaystyle\hskip
25.60747pt\left.-\frac{2}{\xi_{2}^{2}}\Delta_{4}e^{\xi_{1}(d-p+a)}\right\\},$
$\displaystyle c_{3}$ $\displaystyle=$ $\displaystyle
c_{2}e^{2\xi_{1}(p-a)}\Delta_{1}-\frac{\psi}{\xi_{1}^{2}}\Delta_{0}e^{\xi_{1}(p-a)};\;c_{4}=c_{2}+\frac{e^{-\xi_{1}p}}{s_{-}},$
$\displaystyle c_{5}$ $\displaystyle=$
$\displaystyle\Delta_{5}\left(c_{2}+\frac{1}{s_{-}}e^{-\xi_{1}p}\right)e^{2\xi_{1}(d-p-a)}$
$\displaystyle-\Delta_{4}\psi\left(\frac{1}{\xi_{1}^{2}}-\frac{2}{\xi_{2}^{2}}\right)e^{\xi_{1}(d-p-a)},$
$\displaystyle c_{6}$ $\displaystyle=$
$\displaystyle\frac{1}{\gamma_{2}\sinh\left[\xi_{2}\left(\frac{d}{2}-p-a\right)\right]}$
$\displaystyle\times\left\\{(1-\Delta_{5})\left(c_{2}+\frac{e^{-\xi_{1}p}}{s_{-}}\right)e^{\xi_{1}(d-p-a)}\right.$
$\displaystyle\hskip
14.22636pt\left.+\Delta_{4}\psi\left(\frac{1}{\xi_{1}^{2}}-\frac{2}{\xi_{2}^{2}}\right)\right\\},$

with

$\Delta_{4}=\frac{\gamma_{2}\tanh\left[\xi_{2}\left(\frac{d}{2}-p-a\right)\right]}{1+\gamma_{2}\tanh\left[\xi_{2}\left(\frac{d}{2}-p-a\right)\right]};\;\Delta_{5}=1-2\Delta_{4}.$
(120)

## References

  * (1) J. Israelachvili, Intermolecular and Surface Forces, (Academic Press) (1992).
  * (2) R. Podgornik, H.H. Strey and V.A. Parsegian, Molecular Interactions in Lipids, DNA and DNA-lipid Complexes, in Gene Therapy: Therapeutic Mechanisms and Strategies, 209-239 (Marcel Dekker, New York, 2000).
  * (3) Electrostatic Effects in Soft Matter and Biophysics, edited by C. Holm, P. Kekicheff, and R. Podgornik (Kluwer Academic, Dordrecht, 2001).
  * (4) G. L.-Caballero et al., Soft Matter 10, 2805 (2014).
  * (5) M. Delsanti, J. P. Dalbiez, O. Spalla, L. Belloni, and M. Drifford, ACS Symp. Ser. 548, 381 (1994).
  * (6) M. Olvera de la Cruz, L. Belloni, M. Delsanti, J.P. Dalbiez, O. Spalla, and M. Drifford, J. Chem. Phys. 103, 5781 (1995).
  * (7) F. H. J. van der Heyden, D. Stein, K. Besteman, S. G. Lemay, and C. Dekker, Phys. Rev. Lett. 96, 224502 (2006).
  * (8) S. Qiu, Y. Wang, B. Cao, Z. Guo, Y. Chen, and G. Yang, Soft Matter, 2015, 11, 4099-4105.
  * (9) V. A. Parsegian, N. Fuller and R. P. Rand, Proc. Natl. Acad. Sci. USA, 76, 2750 (1979).
  * (10) R. P. Rand and V.A. Parsegian, Biochimica et biophysica acta, 988, 351 (1989).
  * (11) T.J. McIntosh, A.D. Magid, S.A. Simon, Biophys J, 57, 1187 (1990).
  * (12) V.A. Parsegian, R.P. Rand, N.L. Fuller, J. Phys. Chem., 95, 4777 (1991).
  * (13) L. Herrmann, A. Johner, P. Kékicheff, Phys. Rev. Lett., 113, 268302 (2014).
  * (14) M. M. Koerner, L.A. Palacio, J. W. Wright, K. S. Schweitzer, B. D. Ray, and H. I. Petrache, Biophys. J., 101, 362 (2011).
  * (15) A. Averbakha, V.I. Lobyshev, J. Biochem. Biophys. Methods, 45, 23 (2000).
  * (16) L. Fink, C. Allolio, J. Feitelson, C. Tamburu, D. Harries, and U. Raviv, Langmuir 36, 10715 (2020).
  * (17) Marcelja, S. and Radic, N., Chem. Phys. Lett. 42, 129-130 (1976).
  * (18) Leikin, S and Parsegian, V A and Rau, D C and Rand, R P, Annual Review of Physical Chemistry 44, 369-395 (1993).
  * (19) B. Jönsson and H. Wennerström, J.Chem. Soc., Faraday Trans. 2, 79, 19 (1983).
  * (20) R. Kjellander, J. Chem. Soc., Faraday Trans. 2, 80, 1323 (1984).
  * (21) M. L. Belaya, M. V. Feigel’man, and V. G. Levadny, Chem. Phys. Lett. 126, 361 (1986).
  * (22) M. L. Belaya, M. V. Feigel’man, and V. G. Levadny, Langmuir 3, 648 (1987).
  * (23) K. Bohinc, G. V. Bossa, and S. May, Adv. Colloid Interface Sci. 249, 220 (2017).
  * (24) A. Pertsin, D. Platonov, and M. Grunze, Langmuir 23, 1388 (2007).
  * (25) E. Schneck, F. Sedlmeier, and R. R. Netz, Proc. Natl. Acad. Sci. USA 109, 14405 (2012).
  * (26) M. Kanduč, E. Schneck, and R. R. Netz, Langmuir, 29, 9126 (2013).
  * (27) S. Buyukdagli and R. Podgornik, Phys. Rev. E 102, 012806 (2020).
  * (28) G.L. Gouy, J. Phys. 9, 457 (1910).
  * (29) D.L. Chapman, Phil. Mag. 25, 475 (1913).
  * (30) Henderson D., Blum L., Lebowitz J.L., J. Electroanal. Chem. 102, 315 (1979).
  * (31) Henderson D., Blum L., J. Chem. Phys. 75, 2025 (1981).
  * (32) In the case of a pure solvent liquid, evaluating the integral $I$ defined in Eq. (16) with the salt-free potential solutions (D.1) or (74), one obtains the self-consistency condition $Q_{-}=Q_{+}$. The latter equality corresponds to the interfacial electroneutrality constraint originating from the absence of salt.
  * (33) We note that a similar effect has been previously observed in membrane models with fixed interfacial dipoles DipFix .