Unbounded eigenfunctions in the stability problem

for a three-layer flow in porous media

Gelu Paşa

Simion Stoilow Institute of Mathematics, Calea Grivitei 21, Bucharest, Romania

E-mail<EMAIL_ADDRESS>and<EMAIL_ADDRESS>
We study the linear stability of the displacement of three Stokes fluids with
constant viscosity in a porous medium when the middle fluid is contained in a
bounded region. We use the Hele-Shaw model. The eigenfunctions of the
stability system are the amplitudes of the linear perturbations. These
amplitudes must be small. We get unbounded eigenfunctions. So the stability
problem has no physical sense.

AMS 2010 Subject Classification: 34B09; 34D20; 35C09; 35J20; 76S05.

Keywords: Hele-Shaw displacements; Constant viscosity fluids; Hydrodynamic
stability.

1\. Introduction

We consider the linear stability problem for the displacement of three Stokes
fluids in a Hele-Shaw cell, parallel with the horizontal fix plane $x_{1}Oy$
\- see [1], [21]. The three constant viscosities $\mu_{L},\mu,\mu_{R}$ are (in
fact) divided by the constant permeability of the equivalent porous medium.
The fluid $\mu_{L}$ is displacing the middle fluid $\mu$, which in turn is
pushing a fluid with viscosity $\mu_{R}$, in the positive direction $Ox_{1}$.
The middle liquid $\mu$ is contained in a bounded region. We suppose

$\mu_{L}<\mu<\mu_{R}.$ (1)

A large number of successive intermediate liquids can be considered between
the initial displacing fluids with viscosities $\mu_{L},\mu_{R}$. On this way
some authors try to minimize the Saffman-Taylor instability (see [5] \- [12]),
by using the seminal paper [4]. The obtained results have some weak points,
explained in some works of the present author (see [15] -[20]). We give here
only a selective list of related papers. Some have been published recently and
contain more detailed references.

The linear stability system was obtained in [4] and also in [19]. An existence
condition for the non-zero solution of the stability system was used in [4] to
obtain a quadratic equation for the the growth rates; some numerical results
were given. A compatibility condition was derived in [19]; we proved that this
condition is not fulfilled. Thus the stability problem has no solution.
However, some (natural) additional hypothesis were used in [19].

In this paper we use the above existence condition for obtain the behavior of
the eigenfunctions of the stability system in the range of very large
wavenumbers. The new element of this paper is following. We get unbounded
eigenfunctions for increasing wavenumbers. Therefore the perturbations are
also unbounded for large wavenumbers. Thus the stability problem has no
physical sense. This time, we not use any additional hypothesis.

2\. The stability problem

The displacement is given by the velocity $(U,0)$ of the fluid $\mu_{L}$ far
upstream. We use the moving reference $x=x_{1}-Ut$, where $t$ is the time. The
intermediate liquid is contained in the segment $(a,b),b\leq 0$. The relevant
(basic) flow equations are quite similar with Darcy’s law in porous media. A
stationary solution exists, with two straight interfaces $x=a,\,\,x=b$ between
the fluid layers - see [13]. The Laplace-Young law acts on both interfaces,
where two positive surface tensions $T(a),T(b)$ exist.

We perform a linearized stability analysis by normal modes. We use the Fourier
decomposition of the velocity perturbation in the $Ox$ direction:

$u^{\prime}={\color[rgb]{1,0,0}\epsilon}f(x)\exp(iky+\sigma t).$ (2)

Here $k,\sigma,f$ are the wave numbers, growth rates and amplitudes.
$\epsilon$ is a small positive number. The function $f$ must be bounded.
Moreover, $f$ must decay to zero in the far field and is continuous. The
derivative of $f$ is not continuous in $a,b$. In [21] are imposed the
conditions

$f(x)=f(a)e^{k(x-a)},\,\,\forall x\leq a,\,\,k\geq 0;$
$f(x)=f(b)e^{-k(x-b)},\,\,\forall x\geq b,\,\,k\geq 0.$ (3)

We insert (2) in the linearized disturbance equations obtained from the
Darcy’s law. We follow the procedure given in [3], [4] and we obtain the
problem (2a), (2b), (3) of [4]:

$f_{xx}(x,k)-k^{2}f(x,k)=0,\,\,x\in(a,b),\quad b\leq 0,\quad\forall k\geq 0;$
(4) $f_{x}^{+}(a,k)=[-kE_{a}(k)/\sigma(k)+\mu_{L}k/\mu]f(a,k);$ (5)
$f_{x}^{-}(b,k)=[kE_{b}(k)/\sigma(k)-\mu_{R}k/\mu]f(b,k);$ (6)
$E_{a}(k):=\frac{(\mu-\mu_{L})Uk-T(a)k^{3}}{\mu},$
$E_{b}(k):=\frac{(\mu_{R}-\mu)Uk-T(b)k^{3}}{\mu};$ (7)

$f_{x}$ is the derivative; $f_{x}^{-},f_{x}^{+}$ are the lateral derivatives;
$f$ is depending also on $k$.

3\. The non-existence result

The general solution of (4) is

$f(x,k)=A(k)e^{kx}+B(k)e^{-kx}.$ (8)

The coefficients $A,B$ were considered constant in [4]. We use the relations
(5)-(6) and obtain a stability system for $A(k),B(k)$, with variable
coefficients $c(k),d(k),g(k),h(k)$ (see formulas (2b), (3), (5), pag. 112101-5
in [4] and the system (12) of [19]):

$A(k)e^{ka}c(k)+B(k)e^{-ka}d(k)=0,\quad A(k)e^{kb}g(k)+B(k)e^{-kb}h(k)=0;$ (9)
$c(k)=[\mu_{L}-\mu-E_{a}(k)/\sigma(k)],\quad d(k)=[\mu_{L}+\mu-
E_{a}(k)/\sigma(k)],$ (10) $g(k)=[\mu_{R}+\mu-E_{b}(k)/\sigma(k)],\quad
h(k)=[\mu_{R}-\mu-E_{b}(k)/\sigma(k)].$ (11)

A nontrivial solution $(A(k),B(k))$ exists if the following condition is
verified

$e^{2k(a-b)}c(k)h(k)-g(k)d(k)=0,\quad\forall k\geq 0.$ (12)

Our main goal is to find $\lim_{k\to\infty}A(k)$ and $\lim_{k\to\infty}B(k)$
by using (9) - (12).

Remark 1. We have

$d(k),c(k),g(k),h(k)\neq 0,\forall k>0.$ (13)

Indeed, suppose $\exists\,\,m\,\,s.t.\,\,d(m)=0,c(m)=-2\mu$. From
$(\ref{B3})_{1}$ and (8) we get

$A(m)=0,\,f(x,m)=B(m)e^{-mx},\,f_{x}(x,m)=-mB(m)e^{-mx}.$

We emphasize that $B(m)\neq 0$, otherwise $f(x,m)\equiv 0$; the eigenfunctions
can not be identically zero. We simplify with $B(m)$, then (5)-(6) in $k=m$
lead us to

$\mu_{L}e^{-ma}+\mu
e^{-ma}=\frac{E_{a}(m)}{\sigma(m)}e^{-ma};\quad\mu_{R}e^{-mb}+\mu
e^{-mb}=\frac{E_{b}(m)}{\sigma(m)}e^{-mb}.$ (14)

In both relation we have the same $\sigma(m)$ then we get
$E_{a}(m)/(\mu_{L}+\mu)=E_{b}(m)/(\mu_{R}+\mu).$ This is an unexpected
restriction on $\mu$, not considered as a hypothesis.

We get $c(k),g(k),h(k)\neq 0,\forall\,k>0$ by using the same arguments. For
example, suppose $\exists\,\,n\,\,s.t.\,\,c(n)=0$. Then (8), $(\ref{B3})_{1}$
and (5)-(6) in $k=n$ give us $B(n)=0,\,\,f(x,n)=A(n)e^{nx},\,\,A(n)\neq
0,\,\,f_{x}(x,n)=nA(n)e^{nx}$, thus
$E_{a}(n)/(\mu_{L}-\mu)=E_{b}(n)/(\mu_{R}+\mu)$.

$\bullet$

Remark 2. The following properties hold:

$\lim_{k\to\infty}d(k)=0\Rightarrow\lim_{k\to\infty}g(k),\lim_{k\to\infty}c(k),\lim_{k\to\infty}h(k)\quad\mbox{are
real bounded (not zero) quantities}.$ (15)
$\lim_{k\to\infty}g(k)=0\Rightarrow\lim_{k\to\infty}d(k),\lim_{k\to\infty}c(k),\lim_{k\to\infty}h(k)\quad\mbox{are
real bounded (not zero) quantities}.$ (16)

Indeed, suppose $\lim_{k\to\infty}d(k)=0$ . Thus
$\lim_{k\to\infty}E_{a}/\sigma=(\mu_{L}+\mu)$ and we have

$\lim_{k\to\infty}g(k)=\mu_{R}+\mu-(T_{b}/T_{a})(\mu_{L}+\mu),$
$\lim_{k\to\infty}h(k)=\mu_{R}-\mu-(T_{b}/T_{a})(\mu_{L}+\mu),\quad\lim_{k\to\infty}c(k)=-2\mu.$

We use the same arguments and get (16).

$\bullet$

For a function $H(k)$, we use also the notation $H(k)\rightarrow\alpha$ when
$\lim_{k\to\infty}H(k)=\alpha$.

$\bullet$

From Remark 2 it follows that $d(k)\rightarrow 0$ is a possible “solution” of
(12) when $k\rightarrow\infty$. Indeed, the limits of $c(k),h(k),g(k)$ are
bounded, thus $e^{2k(a-b)}c(k)h(k)\rightarrow 0$ and $d(k)g(k)\rightarrow 0$.
Next we prove that $d(k)\rightarrow 0$ or $g(k)\rightarrow 0$ are necessary
conditions.

$\bullet$

Proposition 1. We have the following two possibilities

$or\quad d(k)\rightarrow 0\quad or\quad g(k)\rightarrow 0.$ (17)

Proof. We recall (12). If $c(k)\rightarrow\infty$ or $h(k)\rightarrow\infty$,
then $\lim_{k\to\infty}e^{2k(a-b)}c(k)h(k)$ it might not exist. First, we
prove here that both $c(k),h(k)$ have bounded limits limits to infinity.

Indeed, suppose $c(k)\rightarrow\infty$. That means, in fact,
$E_{a}/\sigma(k)\rightarrow\infty$, because only this term of $c(k)$ is
depending on $k$. We use the notation $I(k):=E_{a}(k)/\sigma$ and get

$e^{2k(a-b)}I^{2}(k)[\frac{\mu_{L}-\mu}{I(k)}-1][\frac{\mu_{R}-\mu}{I(k)}-\frac{E_{b}(k)}{E_{a}(k)}]-I^{2}(k)[\frac{\mu_{L}+\mu}{I(k)}-1][\frac{\mu_{R}+\mu}{I(k)}-\frac{E_{b}(k)}{E_{a}(k)}]=0.$
(18)

We simplify with $I^{2}(k)$. As $E_{b}(k)/E_{a}(k)\rightarrow T_{b}/T_{a}$
(which is a bounded quantity), when $k\rightarrow\infty$ we get
$0-T_{a}/T_{b}=0$. Therefore the equation (12) is not verified when
$c(k)\rightarrow\infty$.

Suppose $c(k),h(k)\rightarrow\infty$. Thus
$I(k):=E_{a}(k)\sigma(k)\rightarrow\infty,\quad
J(k):=E_{b}(k)\sigma(k)\rightarrow\infty$ and from (12) we get

$e^{2k(a-b)}I(k)J(k)[\frac{\mu_{L}-\mu}{I(k)}-1][\frac{\mu_{R}-\mu}{J(k)}-1]-I(k)J(k)[\frac{\mu_{L}+\mu}{I(k)}-1][\frac{\mu_{R}+\mu}{J(k)}-1]=0.$
(19)

When $k\rightarrow\infty$ we get $0-1=0$. Thus (12) is not verified when
$c(k),h(k)\rightarrow\infty$.

The second step is following. Both $c(k),h(k)$ have bounded limits when
$k\rightarrow\infty$, therefore $e^{2k(a-b)}c(k)h(k)\rightarrow 0.$ From (12)
it follows $d(k)g(k)\rightarrow 0$ and we obtain the result (17).

$\bullet$

Remark 3. The above proposition gives us two possible eigenvalues, only for
large $k$:

$k\rightarrow\infty\Rightarrow\sigma_{1}(k)\approx
E_{a}(k)/(\mu_{L}+\mu),\quad\sigma_{2}(k)\approx E_{b}(k)/(\mu_{R}+\mu).$ (20)

This was also obtained in [19] by direct calculations from (12), with the
natural hypothesis

$e^{-2k}k^{3}\approx 0,\,\,e^{-2k}<<k^{3};\quad e^{-2k}k^{6}\approx
0,\,\,e^{-2k}k^{6}<<k^{6}\quad\mbox{ for large enough}\quad k.$

$\bullet$

Proposition 2. There exists unbounded eigenfunctions $f(x,k)$ with $A(k),B(k)$
given by (9).

Proof. Consider $d(k)\rightarrow 0$ and $f(a,k)$ bounded. Then from
$(\ref{B3})_{1}$ we obtain

$(-2\mu)A(k)e^{ka}\rightarrow 0.$ (21)

Let $Y(k)=-g(k)/h(k)$. From (13), (15), (16) it follows that $Y(k)$ is bounded
for large $k$. We use $(\ref{B3})_{2}$ and obtain

$B(k)e^{-ka}=Y(k)A(k)e^{k(2b-a)}.$ (22)

Here is the main point of our paper. We use (22) to get the limit of
$B(k)e^{-ka}$. There exist a lot of possibilities, related with $A(k)$ and the
magnitude of the interval $(a,b)$. For example:

i) if $A(k)=D=constant$ and $a=-10,b=-1$, then $2b-a=8>0$ and
$|B(k)e^{-ka}|\rightarrow\infty$. This contradicts the hypothesis $f(a,k)$
bounded;

ii) if $A=D$ and $a=-2,b=-1,5$ then $2b-a=-1<0$ and both terms of $f(a,k)$ are
convergent to zero for increasing $k$. This is a “good” case.

iii) if $A(k)$ is a polynomial, $a=-10,b=-1$, we still have
$A(k)e^{ka}\rightarrow 0$ and $|B(k)e^{-ka}|\rightarrow\infty.$ This
contradicts again the hypothesis $f(a,k)$ bounded.

Moreover, the formulas (21)-(22) give us

$\lim_{k\to\infty}B(k)e^{-ka}=\lim_{k\to\infty}\\{Y(k)[A(k)e^{ka}]e^{2k(b-a)}\\}.$
(23)

Thus $\lim_{k\to\infty}f(a,k)$ might even not exist.

$\bullet$

Remark 4. Consider the case iii) with $b=0$ and formula (23) with
$A(k)=e^{-ka}$. We get a very strange conclusion. The eigenfunctions (so the
perturbations amplitudes) could be exponentially increasing, for fixed $k$, as
functions of $(b-a)$.  This is in total contradiction with the results of [2],
[3] and [13]. In these three papers, a variable profile was considered in the
intermediate region - denoted by (I.R). On the water-liquid interface, it was
considered a continuous viscosity. The main conclusion was: a large enough
(I.R.) (and a suitable viscosity profile) can almost suppress the instability
Saffman-Taylor. This effect was confirmed by experimental results. On the
contrary, there are no experimental results for the models studied in [4] and
[5] \- [12].

$\bullet$

4\. Conclusions

In the seminal paper [4] were given the following results. From the relations
(5)-(6) was obtained the system (9) with $A,B$ constant. The condition (12)
was used to get a quadratic equation for $\sigma$. Numerical values and some
estimates of $\sigma$ were also given. All these results were used later in
[5] \- [12]. Some weak points exist in these cited papers, which were
explained in [15] \- [19].

In this paper we consider the general case when $A,B$ are depending on $k$.
The eigenfunctions $f(x,k)$ are not bounded (in general) as functions of $k$,
so the amplitudes of perturbations could be not small. Thus the solution of
the stability system has no physical meaning, even if $A=A(k),B=B(k)$. The
equation (12) is still true. But we can say that the corresponding eigenvalues
also do not make physical sense.

## References

  * [1] J. Bear, Dynamics of Fluids in Porous Media. Elsevier, New York, 1972.
  * [2] C. Caraso and G. Paşa, An optimal viscosity profile in the secondary oil recovery, RAIRO M2AN (Modelisation Mathematique et Analyse Numerique), 32, 2, (1998), 211-221.
  * [3] P. Daripa and G. Paşa, An optimal viscosity profile in enhanced oil recovery by polymer flooding. Int. J. Eng. Sci. 42 (2004), 2029-2039.
  * [4] P. Daripa, Studies on stability in three-layer Hele-Shaw flows. Physics of Fluids 20 (2008), 112101.
  * [5] P. Daripa, Hydrodynamic stability of multi-layer Hele-Shaw flows. J. Stat. Mech., Published 10 December 2008.
  * [6] P. Daripa,On estimates for short wave stability and long wave instability in 3-layer Hele-Shaw flows. Physica A: Stat. Mech. Appl., 390 (2011) pp. 3069-3076.
  * [7] P. Daripa and X. Ding, A Numerical Study of Instability Control for the Design of an Optimal Policy of Enhanced Oil Recovery by Tertiary Displacement Processes. Transp. In Porous Media, 93 (2012), 675-703.
  * [8] P. Daripa, On Stabilization of Multi-layer Hele-Shaw and Porous Media Flows in the Presence of Gravity. Transp. In Porous Media, 95 (2012), 349-371.
  * [9] P. Daripa and X. Ding, Universal Stability Properties for Multi-Layer Hele-Shaw Flows and Application to Instability Control. SIAM J. Appl. Math. 72 (2012), 1667-1685.
  * [10] P. Daripa, Some useful upper bounds for the selection of optimal profiles. Physica A: Statistical Mechanics and its Applications, 391 (2012), 4065–4069.
  * [11] P. Daripa and X. Ding, Selection principle of optimal profiles for multi-layer Hele-Shaw flows and stabilization. Transport In Porous Media, 96 (2013), 353-367.
  * [12] P. Daripa and C. Gin, Studies on dispersive stabilization of porous media flows. Physics of Fluids, 28 (2016), 082105\.
  * [13] S.B. Gorell and G.M. Homsy, A theory of the optimal policy of oil recovery by secondary displacement process. SIAM J. Appl. Math. 43 (1983), 79-98.
  * [14] H. S. Hele-Shaw, Investigations of the nature of surface resistence of water and of streamline motion under certain experimental conditions. Inst. Naval Architects Transactions 40 (1898), 21-46.
  * [15] G. Pasa, Some Contradictions in the Multi-Layer Hele-Shaw Flow. Int. Journal of Petroleum Technology, 6 (2019), 41-48.
  * [16] G. Pasa, A paradox in Hele-Shaw displacements. Annali Univ. Ferrara, 66 (2020), 99–111.
  * [17] G. Paşa, The multi-layer Hele-Shaw model with constant viscosity fluids can not minimize the Saffman-Taylor instability. Int. Journal of Applied Mathematics, 33 (2020), 697-708.
  * [18] G. Paşa, Stability analysis for a multi-layer Hele-Shaw displacement. Annali Univ. Ferrara, 68 (2022), 1–9.
  * [19] G. Paşa, On the growth rates for a three-layer flow in porous media. arXiv:2203.11150v1 [math.AP] 18 March 2022.
  * [20] G. Paşa, On the Linear Stability Problem for a Three-Layer Displacement in Porous Media. International Journal of Petroleum Technology, 9 (2022), 8-13
  * [21] P.G. Saffman and G.I. Taylor, The penetration of a fluid in a porous medium or Helle-Shaw cell containing a more viscous fluid. Proc. Roy. Soc. A, 245 (1958), 312-329.


# Learned Compression of Encoding Distributions

Mateen Ulhaq  Ivan V. Bajić

###### Abstract

The entropy bottleneck introduced by Ballé  _et al._ [1] is a common component
used in many learned compression models. It encodes a transformed latent
representation using a static distribution whose parameters are learned during
training. However, the actual distribution of the latent data may vary wildly
across different inputs. The static distribution attempts to encompass all
possible input distributions, thus fitting none of them particularly well.
This unfortunate phenomenon, sometimes known as the amortization gap, results
in suboptimal compression. To address this issue, we propose a method that
dynamically adapts the encoding distribution to match the latent data
distribution for a specific input. First, our model estimates a better
encoding distribution for a given input. This distribution is then compressed
and transmitted as an additional side-information bitstream. Finally, the
decoder reconstructs the encoding distribution and uses it to decompress the
corresponding latent data. Our method achieves a Bjøntegaard-Delta (BD)-rate
gain of −7.10% on the Kodak test dataset when applied to the standard fully-
factorized architecture. Furthermore, considering computational complexity,
the transform used by our method is an order of magnitude cheaper in terms of
Multiply-Accumulate (MAC) operations compared to related side-information
methods such as the scale hyperprior.

Index Terms—  Learned compression, entropy bottleneck, entropy model,
probabilistic modeling, side information

##  1 Introduction

Fig. 1: Visualization of the suboptimality of using a single static encoding
distribution. This distribution is optimal, on average, among all static
distributions, but it is suboptimal for any specific data instance.

In this paper, we explore the compression of probability distributions that
are used in compression. Specifically, we will focus on probability
distributions that are used to encode the latent representations produced by
an image compression model. Entropy models compress latents using _encoding
distributions_, which are discrete probability mass functions (pmfs) used for
encoding. For example, Ballé  _et al._ [1] encode a latent representation
${\boldsymbol{y}}$ using the following entropy models:

  1. 1.

A “fully-factorized” _entropy bottleneck_. Each element ${y}_{c,i}$ within the
$c$-th channel is encoded using a single channel-specific non-parametric
encoding distribution ${p}_{c}(y)$, which is typically fixed once training
completes.

  2. 2.

A _Gaussian conditional_. For each element ${y}_{i}$, the parameters
${\mu}_{i}$ and ${\sigma}_{i}$ of its encoding distribution are computed,
conditioned on the latent
${\boldsymbol{\hat{z}}}=\operatorname{Quantize}[h_{a}({\boldsymbol{y}})]$. In
later works [2], a Gaussian mixture model (GMM) is sometimes used instead,
with parameters ${w}_{i}^{(k)}$, ${\mu}_{i}^{(k)}$, and ${\sigma}_{i}^{(k)}$
for each Gaussian in the mixture.

Of these two, the entropy bottleneck is more flexible at modeling arbitrary
distributions. However, because all elements within a given channel are
modelled using the same distribution, it is quite poor at specializing and
adapting to specific input distributions. Instead, it models a “conservative”
distribution that tries to encompass all possible inputs, which is often
suboptimal, as visualized in fig. 1. This is sometimes known as the
_amortization gap_ [3, 4]. In contrast, the Gaussian conditional is better at
adapting to a specific input, though it is limited to modeling Gaussian
distributions. Furthermore, it must also pay an additional cost for
transmitting the distribution parameters, though this is often worth the
trade-off, e.g., a savings of roughly 40% for a 10% increase in rate in [1].

In this paper, we propose a method to address this amortization gap via the
compression of an input-specific encoding distribution. This allows our
adaptive entropy bottleneck to adapt to the input distribution, rather than
using a static dataset-optimized encoding distribution which suffers from the
amortization gap. We demonstrate that our low-cost input-adaptive method can
be used by models containing the entropy bottleneck in order to achieve better
rate-distortion performance. Our code is available
online111https://github.com/multimedialabsfu/learned-compression-of-encoding-
distributions.

##  2 Related works

In their “scale hyperprior” architecture, Ballé  _et al._ [1] utilize an
entropy bottleneck to compress the latent representation
${\boldsymbol{z}}=h_{a}({\boldsymbol{y}})$. This is known as “side-
information,” and its compressed bitstream is transmitted before the bitstream
for ${\boldsymbol{y}}$. The reconstructed latent ${\boldsymbol{\hat{z}}}$ is
then used to determine the parameters222[1] fixes ${\boldsymbol{\mu}}=0$,
though later papers [5] allow ${\boldsymbol{\mu}}$ to vary.
${\boldsymbol{\mu}}$ and ${\boldsymbol{\sigma}}$ of a Gaussian conditional
distribution which is used to reconstruct ${\boldsymbol{\hat{y}}}$. In many
such models that utilize the scale hyperprior, the rate usage of the entropy
bottleneck often comprises roughly 10% of the total rate. Evidently, the
entropy bottleneck is a key component of many state-of-the-art (SOTA) image
compression models, and so it is important to address the amortization gap
suffered by the entropy bottleneck.

Balcilar  _et al._ [3] propose non-learned methods for the compression and
transmission of encoding distributions. For the factorized entropy model, they
first measure a discrete histogram of the $c$-th latent channel. Then, they
construct an approximation of the histogram using a $K$-Gaussian mixture model
parameterized by $w^{(k)}$, $\mu^{(k)}$, and $\sigma^{(k)}$. The
aforementioned parameters are compactly encoded into 8 bits each, and
transmitted. That is, for a model using $C$ activated channels, $(3K-1)\cdot
C\cdot 8$ bits are required to transmit the parameters. (Typically,
$K\in[1,3]$ is used, and $C\in[16,192]$ for low-rate models.) The
approximation model is then used as the encoding distribution
$p_{{\boldsymbol{\hat{y}}}}({\boldsymbol{\hat{y}}}_{c})=\sum_{k=1}^{K}w_{k}\;\mathcal{\hat{N}}({\boldsymbol{\hat{y}}}_{c};\,\mu_{k},\,{\sigma_{k}}^{2})$.
For the Gaussian conditional entropy model, they instead transmit a factor
$\beta^{(c)}$ that corrects the error in the center (and most likely) bin’s
probability. A key difference between our work and this prior work is that we
formulate a learnable method for the compression of the encoding distribution.

Galpin  _et al._ [6] introduce a lightweight non-learned context model for
raster-scan ordered encoding, similar to Minnen  _et al._ [5]. They compare
their method in conjunction with [7, 3] in various configurations on both the
standard factorized model and its GDN $\to$ ReLU replacement variant. Our work
is also broadly related to source instance overfitting for compression [8, 9].
However, unlike these approaches, we do not overfit the neural network, but
rather the latent distribution of the source instance, which is
computationally more feasible.

##  3 Proposed method

###  3.1 Compression of encoding distributions

Let ${\boldsymbol{x}}$ be an input image, and
${\boldsymbol{y}}=g_{a}({\boldsymbol{x}})$ be its transformed latent
representation, and
${\boldsymbol{\hat{y}}}=\operatorname{Quantize}[{\boldsymbol{y}}]$ be its
quantized version. Consider a fully-factorized entropy model [1], which models
the $j$-th channel of ${\boldsymbol{\hat{y}}}$ (denoted as
${\boldsymbol{\hat{y}}}_{j}$) using a single channel-specific encoding
distribution ${\boldsymbol{\hat{p}}}_{j}$. In this setting, all elements
within a given channel are modelled using the same encoding distribution.
Then, the “true” distribution ${\boldsymbol{p}}_{j}$ can be determined exactly
by computing the normalized histogram of ${\boldsymbol{\hat{y}}}_{j}$.
Notably, ${\boldsymbol{p}}_{j}$ is also the distribution with the minimum rate
cost for encoding ${\boldsymbol{\hat{y}}}_{j}$ using the entropy model. That
is, the ideal encoding distribution ${\boldsymbol{\hat{p}}}^{*}_{j}$ is given
by

${\boldsymbol{\hat{p}}}^{*}_{j}=\operatorname*{arg\,min}_{{\boldsymbol{\hat{p}}}_{j}}\left[H({\boldsymbol{p}}_{j})+D_{\mathrm{KL}}({\boldsymbol{p}}_{j}\parallel{\boldsymbol{\hat{p}}}_{j})\right]\implies{\boldsymbol{\hat{p}}}^{*}_{j}={\boldsymbol{p}}_{j}.$

Unfortunately, using the true distribution ${\boldsymbol{p}}_{j}$ for encoding
is infeasible, since the decoder does not have access to it. Instead, we
propose a method that generates a reconstructed approximation
${\boldsymbol{\hat{p}}}$ that targets ${\boldsymbol{p}}$, subject to a rate
trade-off.

fig. 2 visualizes various collections of encoding distributions used for
encoding different input images using a fully-factorized entropy model. The
left column shows various input images — except for the “(Default)” image,
which is an abstract amortized representation of randomly sampled images from
the training dataset. The middle and right columns visualize the negative log-
likelihoods (in bits) of the true (${\boldsymbol{p}}$) and reconstructed
(${\boldsymbol{\hat{p}}}$) encoding distributions, respectively. Each
collection of distributions is visualized as a color plot, where

  1. (i)

the $x$-coordinate is the channel index $j$ of the latent ${\boldsymbol{y}}$,

  2. (ii)

the $y$-coordinate is the bin index $i$ of the encoding distribution, and

  3. (iii)

the $z$-coordinate (i.e., color intensity) is the negative log-likelihood of
$({\boldsymbol{p}}_{j})_{i}$ or $({\boldsymbol{\hat{p}}}_{j})_{i}$, clipped
within $[0,10]$.

Fig. 2: Visualization of target (${\boldsymbol{p}}$) and reconstructed
(${\boldsymbol{\hat{p}}}$) encoding distributions. Our proposed method
reconstructs ${\boldsymbol{\hat{p}}}$, which is then used by the fully-
factorized entropy model to encode the latent derived from a given input
image. Each collection of distributions is visualized as a color plot, with
channels varying along the $x$-axis, bins varying along the $y$-axis, and
negative log-likelihoods represented by the $z$-axis (i.e., color).

As shown, the default amortized encoding distributions account for a wide
range of possible input distributions. However, the ideal encoding
distributions for any given image often have much narrower ranges. The average
cost of encoding a randomly drawn sample from a distribution
${\boldsymbol{p}}$ is equal to the cross-entropy between that distribution and
the one used for encoding, ${\boldsymbol{\hat{p}}}$. Thus, a mismatch between
the two leads to excessive rate usage. In contrast, with our proposed method,
the encoding distributions more closely approximate the ideal distributions
for encoding a specific image. Our method balances the trade-off between the
rate cost of encoding the encoding distributions ${\boldsymbol{\hat{p}}}$ and
the rate cost of encoding the latent representation ${\boldsymbol{\hat{y}}}$
using the encoding distributions. Indeed, our method accurately reconstructs
regions of high probability, which contribute the most towards saving rate.

###  3.2 Architecture overview

Fig. 3: Adaptive encoding distribution architecture.

Our proposed method applied to the traditional “fully-factorized” architecture
[1] is visualized in fig. 3. In this configuration, the latent representation
${\boldsymbol{y}}$ is first passed through the histogram layer described in
section 3.3, which estimates the target distribution ${\boldsymbol{p}}$. Then,
${\boldsymbol{p}}$ is passed through a distribution compression architecture
consisting of the analysis transform $h_{a,q}$ which generates a latent
representation ${\boldsymbol{q}}$, followed by an entropy model (e.g., entropy
bottleneck), and then finally a synthesis transform $h_{s,q}$ which
reconstructs the distribution ${\boldsymbol{\hat{p}}}$. The rate cost of
encoding ${\boldsymbol{q}}$ is labelled $R_{q}$, and the rate cost of encoding
${\boldsymbol{y}}$ using ${\boldsymbol{\hat{p}}}$ is labelled $R_{y}$.

Like the scale hyperprior from [1], our proposed method also requires the
transmission of side-information from which the encoding distributions are
reconstructed. However, the scale hyperprior applies 2D convolutions to
downsample the 3D input tensor, and later reconstructs encoding distributions
that vary along all 3 dimensions. In contrast, our method initially collapses
the spatial dimensions into a 1D histogram, to which it applies 1D
convolutions, and later reconstructs encoding distributions that vary only
along the channel dimension. This allows our method to be more computationally
efficient when the number of bins and encoding distributions is small.

###  3.3 Histogram estimation

To adapt the encoding distribution to the input, we must first determine the
ideal target encoding distribution. We do this by estimating the histogram for
a given latent channel. To ensure differentiability backwards through the
histogram module, we use the method of kernel density estimation (KDE) [10].

Let $[y_{\mathrm{min}},y_{\mathrm{max}}]$ be a supporting domain that
encompasses all values $y_{1},y_{2},\ldots,y_{N}$ in a specific channel. Let
$\\{b_{1},b_{2},\ldots,b_{B}\\}$ denote the bin centers of the $B$-bin
histogram, where $b_{i}=y_{\mathrm{min}}+i-1$. In the method of kernel density
estimation, which is used to estimate the probability density function from a
set of observations, a kernel function is placed centered at each observed
sample, and those functions are summed to create the overall density function.
To construct the histogram, we evaluate the kernel density estimate at each
$b_{i}$. Or rather, we use an equivalent formulation where we instead place a
kernel function at each bin center. Then, the probability mass in the $i$-th
bin can be estimated as

$(\operatorname{Hist}({\boldsymbol{y}}))_{i}=\frac{1}{\mathrm{Normalization}}\sum_{n=1}^{N}K\left(\frac{y_{n}-b_{i}}{\Delta
b}\right),$

where $\Delta b=b_{i+1}-b_{i}$ is the bin width, and $K(u)$ is the kernel
function, which we choose to be the triangular function

$K_{\mathrm{soft}}(u)=\max\\{0,1-|u|\\}.$

Conveniently, the triangular kernel function ensures that the total mass of a
single observation sums to $1$. Thus, to ensure that the total mass of the
histogram under all $N$ observations sums to $1$, we may simply choose
$\mathrm{Normalization}=N$.

While the above provides a piecewise differentiable estimate of the histogram,
it is also possible to determine the exact _hard_ histogram by using a
rectangular kernel function, which allocates $y_{n}$ to its nearest bin:

$K_{\mathrm{hard}}(u)=\operatorname{rect}(u).$

We can then use the “straight-through” estimator (STE) [11] trick to get the
benefits of both:

$\operatorname{Hist}({\boldsymbol{y}})=\operatorname{Hist}_{\mathrm{soft}}({\boldsymbol{y}})+\operatorname{detach}(\operatorname{Hist}_{\mathrm{hard}}({\boldsymbol{y}})-\operatorname{Hist}_{\mathrm{soft}}({\boldsymbol{y}})).$

Thus, the hard histogram is used for the forward pass, and the soft histogram
for the backward pass (i.e., backpropagation).

###  3.4 Loss function

Our model follows the variational formulation based on [1]; further details on
the construction may be found in earlier works such as [12, 13, 14]. Briefly,
for a variable ${\boldsymbol{y}}$, its counterparts $\tilde{{\boldsymbol{y}}}$
and $\hat{{\boldsymbol{y}}}$ denote the noise-modelled (for training) and
quantized (for inference) versions of ${\boldsymbol{y}}$, respectively.
Furthermore, the variational encoder
$q_{{\boldsymbol{\phi}}}({\boldsymbol{\tilde{y}}},{\boldsymbol{\tilde{q}}}\mid{\boldsymbol{x}})$
with trainable parameters $\phi$ is a probabilistic model of the joint
distribution for ${\boldsymbol{\tilde{y}}}$ and ${\boldsymbol{\tilde{q}}}$,
for a given input ${\boldsymbol{x}}$. Then, the loss function that we seek to
minimize over the input data distribution
$p_{{\boldsymbol{x}}}({\boldsymbol{x}})$ is

$\begin{split}\mathcal{L}&=\operatorname{\mathbb{E}}_{{\boldsymbol{x}}\sim
p_{{\boldsymbol{x}}}({\boldsymbol{x}})}\operatorname{\mathbb{E}}_{{\boldsymbol{\tilde{y}}},{\boldsymbol{\tilde{q}}}\sim
q_{{\boldsymbol{\phi}}}({\boldsymbol{\tilde{y}}},{\boldsymbol{\tilde{q}}}\mid{\boldsymbol{x}})}\left[\begin{multlined}-\log
p_{{\boldsymbol{\tilde{y}}}\mid{\boldsymbol{\tilde{q}}}}({\boldsymbol{\tilde{y}}}\mid{\boldsymbol{\tilde{q}}})\\\
-\lambda_{q}\log p_{{\boldsymbol{\tilde{q}}}}({\boldsymbol{\tilde{q}}})\\\
+\lambda_{x}D({\boldsymbol{x}},{\boldsymbol{\tilde{x}}})\end{multlined}-\log
p_{{\boldsymbol{\tilde{y}}}\mid{\boldsymbol{\tilde{q}}}}({\boldsymbol{\tilde{y}}}\mid{\boldsymbol{\tilde{q}}})\\\
-\lambda_{q}\log p_{{\boldsymbol{\tilde{q}}}}({\boldsymbol{\tilde{q}}})\\\
+\lambda_{x}D({\boldsymbol{x}},{\boldsymbol{\tilde{x}}})\right]\\\
&=\operatorname{\mathbb{E}}_{{\boldsymbol{x}}\sim
p_{{\boldsymbol{x}}}({\boldsymbol{x}})}\left[R_{{\boldsymbol{\tilde{y}}}}({\boldsymbol{x}})+\lambda_{q}R_{{\boldsymbol{\tilde{q}}}}({\boldsymbol{x}})+\lambda_{x}D({\boldsymbol{x}},{\boldsymbol{\tilde{x}}})\right].\end{split}$

For a given input ${\boldsymbol{x}}$,
$R_{{\boldsymbol{\tilde{y}}}}({\boldsymbol{x}})$ and
$R_{{\boldsymbol{\tilde{q}}}}({\boldsymbol{x}})$ are the rate costs of
${\boldsymbol{\tilde{y}}}$ and ${\boldsymbol{\tilde{q}}}$, respectively.
Furthermore, $D({\boldsymbol{x}},{\boldsymbol{\tilde{x}}})$ is a distortion
metric — typically mean-squared error (MSE) — between the input
${\boldsymbol{x}}$ and its reconstruction ${\boldsymbol{\tilde{x}}}$.
Additionally, $\lambda_{x}$ is a trade-off hyperparameter between the rate and
distortion. And lastly, $\lambda_{q}$ is a trade-off hyperparameter between
the rate cost of encoding ${\boldsymbol{\tilde{q}}}$ against the amount of
rate saved when using ${\boldsymbol{\tilde{q}}}$ to encode
${\boldsymbol{\tilde{y}}}$.

The choice for $\lambda_{q}$ depends upon the ratio between the target and
trained-upon input dimensions:

$\lambda_{q}=\frac{H_{{\boldsymbol{x}},\mathrm{trained}}W_{{\boldsymbol{x}},\mathrm{trained}}}{H_{{\boldsymbol{x}},\mathrm{target}}W_{{\boldsymbol{x}},\mathrm{target}}}.$

For instance, for the target dimension of $768\times 512$ (i.e., images from
the Kodak test dataset) and a trained-upon dimension of $256\times 256$ (i.e.,
image patches from the Vimeo-90K dataset), we set $\lambda_{q}=\frac{1}{6}$.

Let ${\boldsymbol{p}}_{j}=(p_{j1},\ldots,p_{jB})$ refer to the discrete
$B$-bin probability distribution used to encode the $j$-th channel of
${\boldsymbol{\tilde{y}}}$. Then concretely, the rate cost of
${\boldsymbol{\tilde{y}}}$ is the cross-entropy between the true and
reconstructed333Since the reconstructed distribution is determined directly
from ${\boldsymbol{y}}$ rather than from ${\boldsymbol{\tilde{y}}}$, we denote
it by ${\boldsymbol{\hat{p}}}$, not ${\boldsymbol{\tilde{p}}}$. discretized
distributions,

$R_{{\boldsymbol{\tilde{y}}}}=\sum_{j=1}^{M}\sum_{i=1}^{B}-p_{ji}\log\hat{p}_{ji}.$

Similarly, let ${\boldsymbol{\tilde{q}}}_{j}$ denote the $j$-th channel of
${\boldsymbol{\tilde{q}}}$. Then, the rate cost of ${\boldsymbol{\tilde{q}}}$
is measured by the typical calculation,

$R_{{\boldsymbol{\tilde{q}}}}=\sum_{j=1}^{M}\sum_{i=1}^{\operatorname{length}({\boldsymbol{l}}_{j})}-{l_{ji}}\log{l_{ji}},\quad\text{where
}{{\boldsymbol{l}}_{j}}=p_{{\boldsymbol{\tilde{q}}}_{j}}({\boldsymbol{\tilde{q}}}_{j}).$

###  3.5 Optimization

We will now compute the derivatives for our proposed method. To simplify the
notation, we will focus on a single channel — the $j$-th channel — and
directly denote ${\boldsymbol{\tilde{y}}}_{j}$ as ${\boldsymbol{y}}$, and
${\boldsymbol{p}}_{j}$ as ${\boldsymbol{p}}$. From the perspective of the
backward pass,

$\begin{split}p_{i}=(\operatorname{Hist}({\boldsymbol{y}}))_{i}&=\frac{s^{2}}{HW}\sum_{n}K_{\mathrm{soft}}{\left(\frac{y_{n}-b_{i}}{\Delta
b}\right)}\\\
&=\frac{s^{2}}{HW}\sum_{n}\max\left\\{0,1-\left|\frac{y_{n}-b_{i}}{\Delta
b}\right|\right\\},\end{split}$

where $s$ is the downscale factor (e.g., $s=2^{4}$), and the dimension of the
input image ${\boldsymbol{x}}$ is $H\times W$. And so,

$\begin{aligned} \frac{\partial p_{i}}{\partial
y_{k}}&=\frac{s^{2}}{HW}\sum_{n}\frac{\partial}{\partial
y_{k}}\left[\max\left\\{0,1-\left|\frac{y_{n}-b_{i}}{\Delta
b}\right|\right\\}\right]\\\ &=\frac{s^{2}}{HW}\frac{\partial}{\partial
y_{k}}\left[\max\left\\{0,1-\left|\frac{y_{k}-b_{i}}{\Delta
b}\right|\right\\}\right]\\\ &=\frac{s^{2}}{HW}\frac{-1}{\Delta
b}\cdot\underbrace{H{\left(1-\left|\frac{y_{k}-b_{i}}{\Delta
b}\right|\right)}}_{1\text{ if }|y_{k}-b_{i}|<\Delta b\text{ else
}0}\cdot\underbrace{\left(2\cdot H{\left(\frac{y_{k}-b_{i}}{\Delta
b}\right)}-1\right)}_{-1\text{ if }y_{k}<b_{i}\text{ else }1}.\end{aligned}$

The corresponding derivative for the reconstructed distribution
${\boldsymbol{\hat{p}}}$ may be determined as

$\begin{split}\frac{\partial\hat{p}_{i}}{\partial
y_{k}}&=\sum_{l}\frac{\partial\hat{p}_{i}}{\partial p_{l}}\frac{\partial
p_{l}}{\partial y_{k}}\approx\sum_{l}\delta_{il}\frac{\partial p_{l}}{\partial
y_{k}}=\frac{\partial p_{i}}{\partial y_{k}},\end{split}$

where we have assumed $\frac{\partial\hat{p}_{i}}{\partial
p_{l}}\approx\delta_{il}$, and that the dominant term in the sum is the term
where $l=i$.

The rate cost of encoding the channel ${\boldsymbol{y}}$ with its associated
encoding distribution ${\boldsymbol{p}}$ is given by

$R_{y}=\frac{HW}{s^{2}}\sum_{i}-p_{i}\log\hat{p}_{i}.$

Then, we may compute the derivative as follows:

$\begin{aligned} \frac{\partial R_{y}}{\partial
y_{k}}&=\frac{-HW}{s^{2}}\sum_{i}\frac{\partial}{\partial
y_{k}}\left[p_{i}\log\hat{p}_{i}\right]\\\
&=\frac{-HW}{s^{2}}\sum_{i}\left[\frac{p_{i}}{\hat{p}_{i}}\frac{\partial\hat{p}_{i}}{\partial
y_{k}}+\frac{\partial p_{i}}{\partial y_{k}}\log\hat{p}_{i}\right]\\\
&\approx\frac{-HW}{s^{2}}\sum_{i}\left[\frac{p_{i}}{\hat{p}_{i}}\frac{\partial
p_{i}}{\partial y_{k}}+\frac{\partial p_{i}}{\partial
y_{k}}\log\hat{p}_{i}\right]\quad\text{assume
}\frac{\partial\hat{p}_{i}}{\partial y_{k}}\approx\frac{\partial
p_{i}}{\partial y_{k}}\\\
&=\frac{-HW}{s^{2}}\sum_{i}\left[\frac{p_{i}}{\hat{p}_{i}}+\log\hat{p}_{i}\right]\frac{\partial
p_{i}}{\partial y_{k}}\\\ &=\frac{-HW}{s^{2}}\sum_{\scriptscriptstyle
i\in\\{\lfloor{y_{k}-y_{\mathrm{min}}}\rfloor+1,\lceil{y_{k}-y_{\mathrm{min}}}\rceil+1\\}}\left[\frac{p_{i}}{\hat{p}_{i}}+\log\hat{p}_{i}\right]\frac{\partial
p_{i}}{\partial y_{k}}\\\ &=\frac{-1}{\Delta
b}\left(\begin{multlined}\left[{\frac{p_{i}}{\hat{p}_{i}}+\log\hat{p}_{i}}\right]_{i=\lceil{y_{k}-y_{\mathrm{min}}}\rceil+1}\\\
\qquad-\left[{\frac{p_{i}}{\hat{p}_{i}}+\log\hat{p}_{i}}\right]_{i=\lfloor{y_{k}-y_{\mathrm{min}}}\rfloor+1}\end{multlined}\left[{\frac{p_{i}}{\hat{p}_{i}}+\log\hat{p}_{i}}\right]_{i=\lceil{y_{k}-y_{\mathrm{min}}}\rceil+1}\\\
\qquad-\left[{\frac{p_{i}}{\hat{p}_{i}}+\log\hat{p}_{i}}\right]_{i=\lfloor{y_{k}-y_{\mathrm{min}}}\rfloor+1}\right)\\\
&\approx\frac{-1}{\Delta
b}\left(\log\hat{p}_{\scriptscriptstyle\lceil{y_{k}-y_{\mathrm{min}}}\rceil+1}-\log\hat{p}_{\scriptscriptstyle\lfloor{y_{k}-y_{\mathrm{min}}}\rfloor+1}\right)\enspace\text{assume
}\hat{p}_{i}\approx p_{i}.\end{aligned}$

(1)

Thus, assuming that $\hat{p}_{i}\approx p_{i}$ and
$\frac{\partial\hat{p}_{i}}{\partial y_{k}}\approx\frac{\partial
p_{i}}{\partial y_{k}}$, the gradient of the rate cost $R_{y}$ is directly
proportional to the difference in code lengths of the two bins whose centers
are nearest to the value $y_{k}$. Intuitively, this makes sense, since the
rate cost of encoding the $k$-th element is directly proportional to the
linear interpolation between the code lengths of the two nearest bins:

$\begin{split}R_{y_{k}}=-\alpha\cdot\log\hat{p}_{\lfloor{y_{k}-y_{\mathrm{min}}}\rfloor+1}-(1-\alpha)\cdot\log\hat{p}_{\lceil{y_{k}-y_{\mathrm{min}}}\rceil+1},\end{split}$

where $\alpha=(y_{k}-b_{\lfloor{y_{k}-y_{\mathrm{min}}}\rfloor+1})/\Delta b$.
The linear interpolation may be justified by the fact that $\hat{y}_{k}$
inhabits only a single bin at a time. The probability of $\hat{y}_{k}$
inhabiting the left bin is $\alpha$, and the probability of inhabiting the
right bin is $1-\alpha$. This aligns with Shannon’s measure of entropy, which
is the expected value of the code length.

For comparison, the standard “entropy bottleneck” represents the likelihood of
a symbol by

$p(y)=c{\left(y+\nicefrac{{1}}{{2}}\right)}-c{\left(y-\nicefrac{{1}}{{2}}\right)}=\int_{y-\nicefrac{{1}}{{2}}}^{y+\nicefrac{{1}}{{2}}}f(t)\,dt,$

where $c$ is a cumulative distribution function and $f(y)=\frac{d}{dy}c(y)$ is
a probability density function, from which the encoding distribution can be
determined by integrating over the bin. Then, the rate cost for the “entropy
bottleneck” is merely the sum of the negative log-likelihoods (i.e., code
lengths),

$R_{y}=\sum_{i}-\log p(y_{i}).$

We may then compute its derivative as follows:

$\begin{split}\frac{\partial R_{y}}{\partial
y_{k}}&=-\sum_{i}\frac{\partial}{\partial y_{k}}\left[\log
p(y_{i})\right]=-\frac{\partial}{\partial y_{k}}\left[\log p(y_{k})\right]\\\
&=-[p(y_{k})]^{-1}\cdot\left[f{\left(y_{k}+\nicefrac{{1}}{{2}}\right)}-f{\left(y_{k}-\nicefrac{{1}}{{2}}\right)}\right].\end{split}$
(2)

Interestingly, whereas the derivative for the proposed method (under the
assumption that $\hat{p}_{i}\approx p_{i}$) computed in 1 contains the
difference between the “right” and “left” log-likelihoods, the derivative for
the entropy bottleneck computed in 2 contains the difference between the
“right” and “left” evaluations of the probability density function.

##  4 Experimental setup

###  4.1 Architecture details

As shown in fig. 4, our $h_{a,q}$ and $h_{s,q}$ are implemented using a simple
five-layer convolutional neural network. In between each of the convolutional
layers shown is a ReLU activation function, as well as a channel shuffle
operation, as is done in ShuffleNet [15]. There are two strides of length $2$,
resulting in a total downscaling factor of $s_{q}=4$. For all our models, we
set $K=15$ to control the kernel sizes, and $G=8$ to control the number of
channel groups. Furthermore, we set $(N_{q},M_{q})=(32,16)$ for low-resolution
models, and $(N_{q},M_{q})=(64,32)$ for high-resolution models. The number of
bins $B$ is set between $128$ and $1024$ across different models. We have
elected to use an entropy bottleneck design for simplicity, though one can
likely further improve the compression performance of the encoding
distribution compression architecture by using more powerful entropy modeling
techniques (e.g., a scale hyperprior). Since $\lambda_{q}$ is a parameter that
depends on the ratio between the target and trained-upon input dimensions, it
may also be advisable to train a single model that supports $\lambda$-rate
control strategies (e.g., G-VAE [16, 17] and QVRF [18]) for both latent
representations $y,q$. However, we have not yet explored this possibility.

Fig. 4: Architecture layer diagram for $h_{a,q}$ and $h_{s,q}$ transforms. $k$
denotes kernel size, $g$ denotes number of channel groups, and
$\downarrow,\uparrow$ denote stride.

###  4.2 Training details

Our models are trained on $256\times 256$ image patches from the Vimeo-90K
triplet dataset [19]. A training batch size of $16$ was used, along with the
Adam optimizer [20] with an initial learning rate of $10^{-4}$ that was
decayed by a factor of $0.1$ whenever the validation loss plateaued.
Specifically, we loaded the weights of the pretrained models from CompressAI
[21], which were also trained using the same setup as above. We replaced the
static distributions of the EntropyBottleneck module with the dynamically
generated adaptive distributions from our proposed encoding distribution
compression module. Then, we froze the weights for $g_{a}$ and $g_{s}$, and
trained only the weights for our encoding distribution compression model
(i.e., for $h_{a,q}$ and $h_{s,q}$, and the entropy model for $q$). Finally,
we evaluated our models on the standard Kodak test dataset [22] containing 24
images of size $768\times 512$.

##  5 Experimental results

fig. 5 shows the rate-distortion (RD) curves comparing a given base model
against the same model enhanced with an encoding distribution compression
method.

Fig. 5: RD curves for the Kodak dataset.

In table 1, we compare the rate savings of our proposed method against the
base model at various quality levels. Additionally, we also list the maximum
rate savings that can be theoretically achieved by perfectly eliminating the
amortization gap between the true and reconstructed distributions, i.e., by
using the true distribution directly as the encoding distribution with zero
additional transmission cost. As shown, our approach achieves a 7.10% BD-rate
reduction over the base model, in comparison to a maximum possible reduction
of 8.21% BD-rate reduction achievable by perfectly eliminating the
amortization gap.

Table 1: Potential and achieved rate savings for the bmshj2018-factorized model [1] when equipped with our proposed distribution compression method. Quality |  Original (bpp) |  Potential gain (bpp) |  Potential gain (%) |  Our (bpp) |  Our gain (bpp) |  Our gain (%)
---|---|---|---|---|---|---
1 | 0.122 | -0.012 | -9.45 | 0.113 | -0.010 | -8.06
2 | 0.188 | -0.016 | -8.66 | 0.174 | -0.014 | -7.53
3 | 0.287 | -0.024 | -8.20 | 0.267 | -0.020 | -7.13
4 | 0.440 | -0.035 | -7.90 | 0.409 | -0.031 | -6.95
5 | 0.647 | -0.051 | -7.83 | 0.602 | -0.045 | -6.98
6 | 0.965 | -0.080 | -8.25 | 0.898 | -0.067 | -6.95
7 | 1.349 | -0.111 | -8.24 | 1.254 | -0.095 | -7.04
8 | 1.830 | -0.149 | -8.14 | 1.705 | -0.126 | -6.88
* |  |  | -8.21 |  |  | -7.10

In order to compute the potential savings, we first ran the base model on each
image from the Kodak test dataset [22], giving us a collection of true
distributions ${\boldsymbol{p}}$ for each image. Furthermore, the base model
provides a static encoding distribution
${\boldsymbol{p_{\textnormal{default}}}}$. Then, we computed the potential
savings (in bpp) for each image as $(\Delta
R)_{\textnormal{max}}=\frac{1}{s^{2}}\sum_{j}D_{\mathrm{KL}}({\boldsymbol{p}}_{j}\parallel({\boldsymbol{p_{\textnormal{default}}}})_{j})$.
We then averaged these values across the dataset.

In table 2, we compare the rate savings for various base models equipped with
a distribution compression method. “Ratio” is the fraction of the total rate
occupied by the original unmodified model. “Gap” is the maximum potential rate
gain with perfect distribution reconstruction and zero additional transmission
cost. “Gain” is the actual rate gain achieved by the model. All models are
evaluated on the Kodak dataset [22].

Table 2: Comparison of rate savings for various models. Model |  Quality |  Factorized |  Total
---|---|---|---
Ratio (%) |  Gap (%) |  Gain (%) |  Gain (%)
bmshj2018-factorized [1] \+ Balcilar2022 [3] | 1 | 100 | -9.45 | -6.79 | -6.79
bmshj2018-factorized [1] \+ ours | 1 | 100 | -9.45 | -8.06 | -8.06
bmshj2018-factorized [1] \+ Balcilar2022 [3] | * | 100 | -8.21 | -5.74 | -5.74
bmshj2018-factorized [1] \+ ours | * | 100 | -8.21 | -7.10 | -7.10

table 3 reports the number of trainable parameters and the number of multiply-
accumulate operations (MACs) per pixel for various model configurations. The
results are calculated assuming an input image size of $768\times 512$, and a
latent representation size of $48\times 32$. As shown, our model requires far
fewer parameters and MACs/pixel than the comparable scale hyperprior [1]
model. In particular, for comparable configurations, our method’s $h_{a,q}$
and $h_{s,q}$ transforms require $96\text{--}97\%$ fewer parameters and
$96\text{--}99\%$ fewer MACs/pixel than the scale hyperprior’s $h_{a}$ and
$h_{s}$ transforms.

Table 3: Trainable parameter counts and number of multiply-accumulate operations (MACs) per pixel.  Model configuration |  Params |  MAC/px |  Params |  MAC/px
---|---|---|---|---
$(M_{y},N_{q},M_{q},K,G,B)$ |  $h_{a,q}$ |  $h_{s,q}$
Ours $(192,32,16,15,8,256)$ | 0.029M | 10 | 0.029M | 10
Ours $(320,64,32,15,8,1024)$ | 0.097M | 126 | 0.097M | 126
$(N,M)$ |  $h_{a}$ |  $h_{s}$
bmshj2018-hyperprior [1] $(128,192)$ | 1.040M | 1364 | 1.040M | 1364
bmshj2018-hyperprior [1] $(192,320)$ | 2.396M | 3285 | 2.396M | 3285

##  6 Conclusion

In this paper, we proposed a learned method for the compression of encoding
distributions. Our method effectively measures and compresses the encoding
distributions used by the entropy bottleneck. The experiments we performed
where we only trained the distribution compression component show that this
method is effective at significantly reducing the amortization gap. Since many
learned compression models use the entropy bottleneck component, our method
provides them with a low-cost improvement in bitrate. Furthermore, our work
opens up the possibility of using learned distribution compression as a
paradigm for correcting encoding distributions.

###  6.1 Future work

A few steps remain towards making encoding distribution corrective methods
viable parts of more advanced entropy models such as [1, 2, 23]. These include
the following:

  * •

The proposed adaptive entropy bottleneck needs to be formulated in such a way
that it can be trained fully end-to-end. When training a pretrained base model
equipped with our adaptive distribution compression module, we found that
unfreezing the pretrained transform weights led to a degradation in RD
performance.

  * •

Application of our adaptive distribution method to the Gaussian conditional
component of entropy models. This component (along with the entropy
bottleneck) is used to construct entropy modeling methods such as those used
in [1, 2, 23]. Most SOTA learned image compression models predict the location
and scale parameters of the Gaussian encoding distributions. There are
potential rate improvements to be made by adapting the shapes of the encoding
distributions to those that more closely match the data distribution.

## References

  * [1] Johannes Ballé, David Minnen, Saurabh Singh, Sung Jin Hwang, and Nick Johnston,  “Variational image compression with a scale hyperprior,”  in Proc. ICLR, 2018.
  * [2] Zhengxue Cheng, Heming Sun, Masaru Takeuchi, and Jiro Katto,  “Learned image compression with discretized gaussian mixture likelihoods and attention modules,”  in Proc. IEEE/CVF CVPR, 2020, pp. 7939–7948.
  * [3] Muhammet Balcilar, Bharath Damodaran, and Pierre Hellier,  “Reducing the amortization gap of entropy bottleneck in end-to-end image compression,”  in Proc. Picture Coding Symposium (PCS). IEEE, 2022, pp. 115–119.
  * [4] Chris Cremer, Xuechen Li, and David Duvenaud,  “Inference suboptimality in variational autoencoders,”  in Int. Conf. on Machine Learning (ICML). PMLR, 2018, pp. 1078–1086.
  * [5] David Minnen, Johannes Ballé, and George D Toderici,  “Joint autoregressive and hierarchical priors for learned image compression,”  Advances in neural information processing systems, vol. 31, 2018.
  * [6] Franck Galpin, Muhammet Balcilar, Frédéric Lefebvre, Fabien Racapé, and Pierre Hellier,  “Entropy coding improvement for low-complexity compressive auto-encoders,”  in Proc. IEEE DCC, 2023, pp. 338–338.
  * [7] Joaquim Campos, Simon Meierhans, Abdelaziz Djelouah, and Christopher Schroers,  “Content adaptive optimization for neural image compression,”  in Proc. IEEE/CVF CVPR Workshops, 2019.
  * [8] Ties van Rozendaal, Iris A. M. Huijben, and Taco S. Cohen,  “Overfitting for fun and profit: Instance-adaptive data compression,”  in Proc. ICLR, 2021, pp. 1–18.
  * [9] Honglei Zhang, Francesco Cricri, Hamed Rezazadegan Tavakoli, Maria Santamaria, Yat-Hong Lam, and Miska M. Hannuksela,  “Learn to overfit better: finding the important parameters for learned image compression,”  in Proc. IEEE VCIP, 2021, pp. 1–5.
  * [10] Emanuel Parzen,  “On Estimation of a Probability Density Function and Mode,”  The Annals of Mathematical Statistics, vol. 33, no. 3, pp. 1065–1076, 1962.
  * [11] Yoshua Bengio, Nicholas Léonard, and Aaron Courville,  “Estimating or propagating gradients through stochastic neurons for conditional computation,” 2013.
  * [12] Diederik P Kingma and Max Welling,  “Auto-encoding variational bayes,” 2013.
  * [13] Johannes Ballé, Valero Laparra, and Eero P. Simoncelli,  “End-to-end optimized image compression,”  in Proc. ICLR, 2017.
  * [14] Lucas Theis, Wenzhe Shi, Andrew Cunningham, and Ferenc Huszár,  “Lossy image compression with compressive autoencoders,”  in Proc. ICLR, 2017.
  * [15] Xiangyu Zhang, Xinyu Zhou, Mengxiao Lin, and Jian Sun,  “ShuffleNet: An extremely efficient convolutional neural network for mobile devices,”  in Proc. IEEE/CVF CVPR, 2018, pp. 6848–6856.
  * [16] Ze Cui, Jing Wang, Shangyin Gao, Bo Bai, Tiansheng Guo, and Yihui Feng,  “G-VAE: A continuously variable rate deep image compression framework,” 2020.
  * [17] Ze Cui, Jing Wang, Shangyin Gao, Tiansheng Guo, Yihui Feng, and Bo Bai,  “Asymmetric gained deep image compression with continuous rate adaptation,”  in Proc. IEEE/CVF CVPR, 2022, pp. 10532–10541.
  * [18] Kedeng Tong, Yaojun Wu, Yue Li, Kai Zhang, Li Zhang, and Xin Jin,  “QVRF: A quantization-error-aware variable rate framework for learned image compression,”  in Proc. IEEE ICIP, 2023.
  * [19] Tianfan Xue, Baian Chen, Jiajun Wu, Donglai Wei, and William T Freeman,  “Video enhancement with task-oriented flow,”  Int. J. Comput. Vis. (IJCV), vol. 127, pp. 1106–1125, 2019.
  * [20] Diederik P. Kingma and Jimmy Ba,  “Adam: A method for stochastic optimization,” 2015.
  * [21] Jean Bégaint, Fabien Racapé, Simon Feltman, and Akshay Pushparaja,  “CompressAI: a PyTorch library and evaluation platform for end-to-end compression research,” 2020.
  * [22] Eastman Kodak,  “Kodak lossless true color image suite (PhotoCD PCD0992),” .
  * [23] Dailan He, Ziming Yang, Weikun Peng, Rui Ma, Hongwei Qin, and Yan Wang,  “ELIC: Efficient learned image compression with unevenly grouped space-channel contextual adaptive coding,”  in Proc. IEEE/CVF CVPR, 2022, pp. 5718–5727.


# Water inside charged nanoslits: Structure and Dielectric study with a novel
water model FAB/$\epsilon$

Raúl Fuentes-Azcatl<EMAIL_ADDRESS>José Rafael Bordin
Departamento de Física, Instituto de Física e Matemática, Universidade Federal
de Pelotas. Caixa Postal 354, 96001-970, Pelotas, Brazil.
<EMAIL_ADDRESS>Marcia C. Barbosa<EMAIL_ADDRESS>Instituto de
Física, Universidade Federal do Rio Grande do Sul, Caixa Postal 15051, CEP
91501-970, Porto Alegre, RS, Brazil

###### Abstract

In this work, the dielectric behavior of water inside charged nanoslit of
graphene is studied to analized the water molecules under electrical
confinement; through polarizing the nanoslit of graphene, creating an electric
field inside the nanopore. How the water molecules are structured under this
type of electrical confinement is studied with two force fields of water, the
three-site water models here used are the SCP/$\epsilon$ and the
FAB/$\epsilon$, the first is a rigid model that improves the SPC model and the
second is a flexible model that improves all the force fields of three sites
non-polarizables and flexibles.

Departamento SIC. Universidad Autónoma Metropolitana, Lerma. Av. de las Garzas
10, El panteón, 52005 Lerma de Villada, Edomex. México.

##  1 Introduction

Water is unique not only because life, as we known, happens in water but also
due to its relevance in industry. It exhibits more than 70 known anomalous
behavior 1, and the water behavior becomes even more stranger when confined in
pores with few nanometers 2, 3, 4. Specially, the comprehension of confined
water is important from many aspects. Water confined inside biological
proteins 5, in the DNA grooves 6 or inside rocks nanopores, as it is in the
shale formation process 7 presents additional anomalies and modify in a non
trivial way the confined media. Also, many new technologies have been proposed
based in the water behavior at this scale 8. Among the many nanomaterials
employed in these new technologies 9, single layers of graphene (SLG) emerge
as one of the most prominent material 10. It was successfully used in a wide
range of applications such as water filtration 11, 12, single-molecule sensors
13, 14, ion selectivity 15, 16, 17, and energy conversion and storage 18, 19,
20, 21.

Water ordering and orientation at the nanopore interface play a major role in
the nanofluidic research and all the possible applications. For instance,
recent experimental 22, 23 and computational 24, 25, 26, 27, 28 studies have
indicated an confinement induced reduction in the water dielectric constant
that can be related to the water structure at the interface or to long-ranged
anisotropic dipole correlations combined with an excluded-volume effect of the
low-dielectric pore material. Also, the surface layer structure can
drastically affect the dynamical behavior of nanoconfined water 29, 30, in
fact, the dynamics and the static dielectric constant are related as well 31.
Therefore, understand the water behavior at the solid boundaries under extreme
confinement is of great importance (and a great challenge) to fundamental
research 32 and technological applications.

Despite the huge advances in nanotechnology, experiments of fluids confined at
the length scale of few nanometers remain as a considerable challenge 33, 34,
35. Is this sense, the use of Molecular Dynamics (MD) simulations arises as a
suited approach to study nanofluidic systems 36. Simulations are used to
interpret experiments, but also to study properties that can not be accessed
directly by them. The proper depiction by MD simulations relies on the
definition of the Classical Force Fields based in additive pairwise potentials
for the particle particles interaction 37, 38.

In the specific case of water, a large number of atomistic models have been
proposed. Most of them are rigid nonpolarizable models, with 3 or more sites.
39, 40, 41 Even this simple models can reproduce most of the bulk water
behavior, there are recent models that can reproduce the structure and
infrared spectrum as the FBA/$\epsilon$ and TIP4P/$\epsilon$ flexible 42, 43,
improving all current non-polarizable and flexible models 43. Properly include
polarization in classical simulations remain as an huge task, but efforts are
being done in this direction in the recent years 44, 45, 46, but the
flexibility offers a type of polarization to the molecule giving freedom of
movement to the charge, as in the model FBA/$\epsilon$ 42.

Also, the addition of the molecule flexibillity is important particularly in
the study of water solutions. For instance, the H–O–H angle is 104,474o for
water molecules in the gas phase, with a distance of 0.095718 nm between
oxygen and each hydrogen. However, in the liquid phase at 298 K and 1 bar, the
average angle changes to 106o 47. Once each water molecule, in bulk, form up
to four hydrogen bonds, small variations in the H–O–H angle can be relevant to
the molecules cluster structure. Under extreme confinement, where the hydrogen
bond is suppressed by the geometrical constraints, variations in the angle can
be even more relevant for the confined water structural, dynamic,
thermodynamic and dielectric properties.

A drawback of the classical atomistic models fitted using bulk properties is
that distinct classical models, even though agree in some measures in the
confined geometries, lead to different quantitative results in other
quantities 48, 49, 50. For instance, Celebi and co-workers 30 recent analyzed
the slip length of water inside graphene nanopores employing rigid non-
polarizable models and found considerable differences and deviations in the
viscosity and slip length when compared with experimental results. With this,
a question that naturally arises is how the molecule flexibillity can improve
the reproduction of the behavior of nanoconfined water when compared with the
classical atomistic approach.

In order to answer this question we here study the behavior of two water
models confinement between two parallel single layers graphene (SLG) sheets.
The SLG walls can be neutral, modeling a nanopore made by pristine graphene,
or charged, as a minimal model for a nanocapacitor. This charged system is not
only interesting from the practical applications of nanocapacitors, but to
understand how the charged wall and electric field, modify the structure of
the water and how it will reflect in the fluid behavior.
We compare the behavior under confinement of two recent three site water
models that improved the traditional SPC water model 51. The rigid
SPC/$\epsilon$ model 40, was constructed based on the SPC model, but show a
with better depiction of the dielectric behavior. More recently, the flexible
FAB/$\epsilon$ 42 model was constructuded to enhance the thermodynamic
properties at different pressure and temperature conditions in comparison to
all three-site models by adding extra degrees of fredom to the sites. With
this, we can analyze how variations in the internal angle and bond lengths
will affects the behavior of the interfacial and bulk-like water inside
neutral and charged nanopores.

Our paper is organized as follows. In Section 2 we introduce the water and
nanopore models and the details about the simulation method. In Section 3 we
show and discuss our results, and the conclusions are presented in Section 4.

##  2 The Models

###  2.1 The Water Models

In this work we compare two water models, one rigid and one fully flexible.
The rigid model, namely SPC/$\epsilon$ 40 model, is an improved SPC model with
better depiction of the dielectric behavior which helped to have better force
fields of monovalent ions 52, 53. The another water model is the recently
proposed FAB/$\epsilon$ 42, flexibility helps to study changes in structure
due to other substances such as ionic liquids 54, fully flexibility of the
model: it has harmonic potentials at the O-Hbond and at the angle formed by
the three molecule site,which helps to improved the force field as in the CO2
molecule 55.

While in the SPC/$\epsilon$ model the atoms in each molecule are maintained
connected by rigid constraints, in the FAB/$\epsilon$ flexible model there is
the addition of the harmonic potential in the O-H bond,

$U_{k}(r)=\frac{k_{r}}{2}(r-r_{0})^{2}$ (1)

and in the H-O-H angle,

$U_{\theta}(\theta)=\frac{k_{\theta}}{2}(\theta-\theta_{0})^{2},$ (2)

where $r$ is the bond distance and $\theta$ is the bond angle. The subscript
$0$ denotes their equilibrium values, and $k_{r}$ and $k_{\theta}$ are the
corresponding spring constants.

For the intermolecular potential between two molecules the LJ and Coulomb
interactions are used for non-polarizable models,

$u(r)=4\epsilon_{\alpha\beta}\left[\left(\frac{\sigma_{\alpha\beta}}{r}\right)^{12}-\left(\frac{\sigma_{\alpha\beta}}{r}\right)^{6}\right]+\frac{1}{4\pi\epsilon_{0}}\frac{q_{\alpha}q_{\beta}}{r}$
(3)

where $r$ is the distance between sites $\alpha$ and $\beta$, $q_{\alpha}$ is
the electric charge of site $\alpha$, $\epsilon_{0}$ is the permitivity of
vacuum, $\epsilon_{\alpha\beta}$ is the LJ energy scale and
$\sigma_{\alpha\beta}$ the repulsive diameter for an $\alpha-\beta$ pair. The
cross interactions between unlike atoms are obtained using the Lorentz-
Berthelot mixing rules,

$\sigma_{\alpha\beta}=\left(\frac{\sigma_{\alpha\alpha}+\sigma_{\beta\beta}}{2}\right);\hskip
28.45274pt\epsilon_{\alpha\beta}=\left(\epsilon_{\alpha\alpha}\epsilon_{\beta\beta}\right)^{1/2}$
(4)

The table 1 shows the force field values of the models mentioned above.

Table 1: Parameters of the three-site water models considered in this work.  model | $k_{b}$ | $r_{OH}$ | $k_{a}$ | $\Theta$ | $\varepsilon_{OO}$ | $\sigma_{OO}$ | $q_{O}$ | $q_{H}$
---|---|---|---|---|---|---|---|---
| kJ/ $mol$ Å2 | Å | kJ/ $mol$ rad2 | deg | kJ/mol | Å | e | e
SPC/$\epsilon$ 40 | - | 1.000 | - | 109.45 | 0.705859 | 3.1785 | -0.8900 | 0.4450
FBA/$\epsilon$ 42 | 3000 | 1.027 | 383 | 114.70 | 0.792324 | 3.1776 | -0.8450 | 0.4225

###  2.2 The Graphene Model

The graphene walls were modeled as a highly oriented pyrolytic graphite (HOPG)
sheet. It is a highly pure and ordered form of synthetic graphite. It is
characterized by a low mosaic spread angle, meaning that the individual
graphite crystallites are well aligned with each other. The original force
field for the carbons in this graphene model is composed only by a Lennard-
Jones interaction 56. Here, charges are included in every carbon atom to have
an uniform surface charge densities, $\sigma_{S}$, as given by Shim et al. 57.
The charges are increased from zero (neutral confinement) to create the
desired voltage between the walls.

Table 2: Force field parameters of grafene model.

Model | $q_{C}/e$ |  $\sigma$/Å |  $(\epsilon/k_{B})$/K
---|---|---|---
C | 0.0143† | 3.401 | 35.46

$\dagger$-example of charge for a uniform surface charge densities of
$\sigma_{S}$= 0.86e/nm2

##  3 Simulation Details

We employed Molecular Dynamics (MD) simulation in the $NVT$ ensemble to study
the behavior of the confined water inside a nano-capacitor. The nano-capacitor
was constructed with two parallel SLG with 240 carbon atoms placed in an area
Lx =2.39315 Ly = 2.39315 nm separated by a distance of Lz = 2.85 nm. Then, 458
H2O molecules were initially placed randomly inside the nano-capacitor space.
Simulations with electric potential ranging from 0 to 150 V were performed.
The electrical potential was created by inducing a surface charge density in
the SLG parallel plates by including point charges in each carbon atom in
order to get the desired value. The left wall is the positive plate and the
right wall the negative one. The net charge of the system is zero.

All simulations were conducted using the molecular dynamics package GROMACS
201658. Periodic boundary conditions were imposed in the $xy$-plane. The
equations of motion were solved using the leap-frog algorithm 37, 58 with a
time step of 1 fs. The temperature was set to $T=$ 298 K with the Nosé-Hoover
thermostat 59 with coupling parameter 0.6 ps. The LINCS algorithm was employed
for the ridig bons. The Lennard-Jones potential was truncated at 10 Åwithout
any long-range correction in energy or pressure. The real part of the Coulomb
potential was also truncated at 10 Å. The Fourier component of the Ewald sums
was evaluated with the particle mesh Ewald (PME) method 60 using a grid
spacing of 0.35 Å and a three degree polynomial for the interpolation. All
simulations were run for 50 ns, were the first 5 ns are the equilibration time
and the next 45 ns the production time, with a output frequency for positions
and velocities at each 1000 steps.

##  4 Results

(a)  (b)

Figure 1: Density profile along the Z direction for the FAB/$\varepsilon$
model (a) and for the SPC/$\varepsilon$ model (b) for various potentials
between the plates. Three regions are defines: region I, near the left wall,
region II, the bulk-like region at the center, and region III, near the right
wall.

Water (and fluids in general) under confinement assumes a layered structure
near the walls and have a bulk-like region in the nanopore center 61, 62, 3,
63, 64, 65. We observe the same behavior for both water models as illustrated
by the density profiles along the $z$-direction for the distinct applied
voltages are shown in figure 1(a) for the flexible FAB/$\varepsilon$ model and
in graph 1(b) for the rigid the SPC/$\varepsilon$ model. The water central
density is the value expected by the bulk H2O density of each model for the
specific temperature and pressure, and the density profiles are similar to the
observed in previous works for water confined by neutral and charged
nanoplates 3, 66, 67, 68, 69. The same trends are observed for both H2O
models. To understand the distinct water behaviors inside the nanopore, we
divided it in three regions: the contact layers in the region I - near the
positive wall - and in the region III - at the negative wall and the bulk-like
region II.

(a)  (b)

(c)  (d)

(e)

Figure 2: (a) Dipole moment (b)H–O–H average angle (c) O-H average length for
confined H2O versus voltage in the distinct regions between the plates for the
SPC/$\varepsilon$ rigid model (square) ) and FAB/$\varepsilon$ flexible
(circles) water models. (d) Schematic view of the water near the neutral
nanopore walls and close to the charged walls.

The distinct occupancy and layering in the regions defined by the slabs I, II
and III can lead to distinct physical properties for the water molecules. In
order to understand the role of model flexibility in water properties in each
region we evaluated the molecule dipole moment $\mu$ and variation in the
internal bonds and angles. These quantities are constant for rigid models and
varies for flexible models. The dipole moment versus voltage for the different
regions I, II and III, are shown in figure 2(a), while the H–O–H internal
angles and the H–O bond length are shown in the figures 2(b) and (c),
respectively. The results indicates that for the FAB/$\varepsilon$ model $\mu$
changes as the voltage increases in a non monotonic way. In the region II,
$\mu$ presents values slightly above the bulk with a small increase with the
voltage. In the region I, $\mu$ increases with the voltage. In the region II,
$\mu$ decreases with the voltage below a threshold, $V<V^{*}$ ($V^{*}$ is a
threshold voltage), and increases above it.

The behavior of $\mu$ correlates with changes in the H-O-H angle, $\theta$ and
the O-Hbond length as shown in figure 2(b) and (c). Near the cathode, in the
slab region I, $\theta$ increases as $V$ increases for $V<V^{*}$, a
consequence of the repulsion between the H’s atoms and the positive wall.
Likewise, this repulsion compress the hydrogens in direction to the O atom,
leading to a decrease of the O-Hbond length for $V<V^{*}$, as shown in the
figure 2(c). These structural changes lead to the decrease in $\mu$ observed
in the region I, $V<V^{*}$ in figure 2(a). Above this voltage threshold, $\mu$
increases with $V$. This change in the slope is followed by the structural
changes. As $V$ grows above this threshold, $\theta$ decreases and the O-Hbond
length increases.

The inflection point in $V=V^{*}$ coincides with the voltage where the peak in
the density profile moves closer to the charged wall, as we can see in figure
1(a) for $V=75$ V and 150 V. This indicates that the negatively charged O
atoms are now closer to the wall, increasing the wall repulsion to the H atoms
and increasing the bond length and the bond angle.

Near the anode, in the region III, we do not see the monotonic behavior of
$\mu$, $\theta$ and the H–O bond length. $\mu$ increases with $V$, as the bond
length also increases and $\theta$ decreases, reinforcing the relations
observed near the positive plate. Here, also the increase in $\mu$ is
consequence of the molecules going closer to the negative wall. However, here
the H atoms are attracted. As the graphene wall pull them in her direction,
the bond length increases and $\theta$ decrease. The elongation in the bond
length lead to the increase in the dipole moment. Also, is interesting to
mention that the confinement changes the bond length (and $\mu$) even for the
bulk-like region II and neutral walls case. As reported by Fuentes and Barbosa
42, the bond length for the FAB/$\varepsilon$ in the liquid phase at 298 K and
1 bar is 0.9495 Å. However, for zero volts, the average bond length in the
region II is 1.058 Å, this indicates that the molecule internal structure in
the bulk-like regime is affected by the confinement.

Figure 3: Infrared spectrum of FAB/$\varepsilon$ flexible model around
different voltages is compared with the experimental value of bulk liquid at
298K and 1 bar of temperature and pressure respectively.

The changes in the intermolecular bonds length and angle can also be observed
in the Infrared (IR) spectrum of the flexible model. Analyzing the behavior of
the IR properties with respect to the increase of the electrical potential, we
can see in figure 3 that there are small differences that corresponds to the
frequency at the maximum for the angle bending, the bond stretching and the
bands. It’s corresponds to the librational motion with symmetries A2, B2. As
$V$ increases, the intermolecular stretch definition is lost as shown in Table
3. This leads to the molecule becoming less flexible and the cage vibrations,
that are located at 50 cm-1, becomes more defined as the applied tension
increases.

Table 3: Experimental and simulation IR data results of FAB/$\varepsilon$
water model42. In liquid phase at 298K and 1bar of temperature and pressure
respectively. Wave numbers(cm -1 ) at the peak of the bands of the power
spectrum

| Exp. | FBA/$\varepsilon$
---|---|---
Cage vibrations | 50 | 50
Intermolecular stretch | 183.4 | 278
Librations A2 | 430 | 536
Librations B2 | 650 | 685
Bending (H-O-H) | 1643.5 | 1613
Stretching (O-H) | 3404 | 3060

The polarization factor Gk measures the equilibrium fluctuations of the
collective dipole moment of the system and is related to the orientational
correlation function. Kirkwood postulate that it should be possible to express
$\varepsilon$ in terms of a short-range orientational correlation function 70.
In this way, here we evaluate the polarization factor GK 71,

$G_{K}=<{\bf M}^{2}>/N\mu^{2}$ (5)

Here, $N$ is the total number of molecules, $M$ is the sum total of dipoles
${\mu}$ in the system. Local orientational correlations are averaged out by
thermal motion after the first few coordination shells.

Analyzing the polarization values Gk and the average dielectric constant
$\epsilon$ inside the pore as seen in the figure 4, the polarization have a
great increment up the V* and from there it decreases, as the dielectric
constant.

(a)  (b)

Figure 4: (a) Kirkwood factor Gk and (b) Total dielectric constant as function
of the electrical potential for the two H2O models. Black circles stands for
the FAB/$\varepsilon$ flexible model and red squares for the SPC/$\varepsilon$
rigid model.

The changes in the water distinct dipole moment and the suppression of the
molecule flexibility observed in the IR influences how the water behaves as a
dielectric medium. The static dielectric constant $\varepsilon$ of a polar
liquid is related to the thermal equilibrium fluctuations of the polarization
at zero field 72. Polarization fluctuations are long-range and vary with the
shape of the dielectric body therefore, the confinement strongly affects the
long-range interaction. We can see in the figure 4 that both flexible and
rigid water have a similar behavior. For low voltages the dielectric constant
decreases with voltage and reaches a minimum value. For voltages above
$V_{min}$, the dielectric constant increases with voltage to a maximum,
$V=V*$. For $V>V*$ the dielectric constant decreases with voltage. For the
flexible model the maximum is located at $V^{*}$. The maximum values in
$\varepsilon$ for the FAB/$\varepsilon$ model is higher than the observed for
the rigid SPC/$\varepsilon$ water. This is a consequence from the fact that
for rigid models, variations in the dielectric permitivity are related to
distinct orientations of the molecules, changing the molecular arrangement.
For flexible models, besides the orientations, the variation in the internal
bonds length and angle leads to stronger voltage effects in the dielectric
behavior.

Is interesting to compare the total dielectric constant in bulk and confined
for both water models. For the rigid model SPC/$\varepsilon$, the bulk
dielectric constant is $\varepsilon_{\rm bulk}=78$ 40. The comparison between
this value and the dielectric contant of confined SPC/$\varepsilon$ as shown
in figure 2 indicates that the confinement decreases the dielectric constant.
This result agrees with previous computational and experimental studies 73,
23, 74, 25, 75, 76 that indicates that water dipoles tends to orientate
parallel to the confining surfaces, suppressing the parallel contribution to
the total dielectric constant and the total dielectric constant itself.

For FAB/$\varepsilon$ water we observe a distinct scenario. This model have,
in bulk, a dielectric constant of $\varepsilon=76$. As we can see in the
figure 4, the total dielectric constant in the confined regime is smaller than
the bulk value for low applied voltages. However, in the voltage regime where
the dipole moment of the interfacial water molecules in regions I and III
sharply increases, above 50V, the total dielectric constant increases,
becoming higher than the bulk value. This shows that not only the orientation
of the water molecule at the nanopore interface affects the dielectric
constant, but at high voltages the internal bonds can play a major role.

Finally, we analyze the dynamical behavior. Recent works have explored the
relation between the dynamics of the contact layer and their electric
properties. Experiments by Fumagalli and co-workers 23 have shown that the
contact layer - with 3 molecules thickness - is electrically dead, with out-
of-plane $\varepsilon$ of only 2. Using simulations of the rigid SPC/E water
model, Mondal et al observed that, despite electrically dead, this layers is
not dynamically dead near a hydrophobic surface 31. For our hydrophobic case,
when the voltage is 0V, the interfacial layer for the rigid SPC/$\varepsilon$
water has a slightly higher diffusion than the flexible FBA/$\varepsilon$, as
we can see in the figure 5, despite of have the same dielectric constant.
Interesting, the diffusion of the contact layers for the SPC/$\varepsilon$
model reaches a minima at 22.5 V - like the dielectric constant. However, for
the flexible model $D$ increases, diffusing faster than the non-confined case,
whose diffusion constant is 2.39 $cm^{2}s^{-1}$ at 300 K and 1 bar. This can
be related also to the dielectric constant: the FBA/$\varepsilon$ water
$\varepsilon$ decreases less in comparison to the rigid model. As the voltage
increases even more, both models achieves higher diffusion than the non-
confined case for the SPC/$\varepsilon$ model is 2.4 $cm^{2}s^{-1}$. This
increase in the diffusion constant is followed by the increase in
$\varepsilon$, reinforcing the relation between this two quantities for water
under extreme confinement.

Figure 5: Diffusion coefficients for the H2O confined between two charged
walls calculated in the different regions defined in figure 1 and for
different electric potentials. Black circles FAB/$\varepsilon$ flexible model
and red squares with the SPC/$\varepsilon$ rigid model.

##  5 Conclusions

Molecular Dynamics simulations were carried out to study confinement of H2O
between two parallel graphene walls with a surface charge distribution, i.e a
electric potential. Two different H2O models were tested, the FAB/$\epsilon$
flexible and SPC/$\epsilon$ rigid models. At zero electric potential both
models present similar structure along the nano-capacitor and when the
electric potential is increased a strong structure is formed close to the
walls. In particular the FAB/$\epsilon$ model can adapt their structure and
give us more information about this changes and if they are representative in
the dielectric properties as the electric field increases.

Due to the dipole moment induced at very strong electrical potentials, the
dielectric constant changes within the nanopore and, therefore, the
electrostatic properties are expected to change. We see a non-monotonic
behavior of the total dielectric constant as function of the applied voltage,
with a minima at 22.5 V. This minima is more pronounced for the rigid
molecule, and is followed by a minima in the diffusion at the interfacial
layer for this water model.The flexible molecule has a less pronounced minima
in $\varepsilon$ and did not show a minima in $D$. This findings helps to
elucidate the role of the water model to understand the behavior of fluids
confined at extreme conditions.

RFA acknowledge support from CONACYT, SECTIE, UAM Supercomputing Center
(Yotla) for the computer time allocated for some of our calculations. MCB and
JRB thanks the Brazilian science agencies - Conselho Nacional de
Desenvolvimento Cientıfico e Tecnologico (CNPq) (INCT-Fc) and Coordenacāo de
Aperfeiçoamento de Pessoal de Nível Superior (CAPES Print Program)

## References

  * Chaplin 2020 Chaplin, M. Anomalous properties of water. http://www.lsbu.ac.uk/water/anmlies.html, 2020
  * Zangi 2004 Zangi, R. Water confined to a slab geometry: a review of recent computer simulation studies. _Journal of Physics: Condensed Matter_ 2004, _16_ , S5371
  * Gao et al. 2016 Gao, Z.; Giovambattista, N.; Sahin, O. Phase Diagram of Water Confined by Graphene. _Schientific Reports_ 2016, _8_ , 6228
  * Muñoz-Santiburcio and Marx 2021 Muñoz-Santiburcio, D.; Marx, D. Confinement-Controlled Aqueous Chemistry within Nanometric Slit Pores. _Chemical Reviews_ 2021, _121_ , 6293–6320, PMID: 34006106
  * Lynch et al. 2020 Lynch, C. I.; Rao, S.; Sansom, M. S. P. Water in Nanopores and Biological Channels: A Molecular Simulation Perspective. _Chemical Reviews_ 2020, _120_ , 10298–10335, PMID: 32841020
  * Passos et al. 2014 Passos, C.; Kuhn, P.; Barbosa, M. DNA denaturation in the rodlike polyelectrolyte model. _Physica A: Statistical Mechanics and its Applications_ 2014, _413_ , 481–488
  * Sun et al. 2019 Sun, F.; Yao, Y.; Li, G.; et al., Effect of critical thickness on nanoconfined water fluidity: review, communication, and inspiration. _J. Petrol. Explor. Prod. Technol._ 2019, _9_ , 1149–1159
  * Bampoulis et al. 2018 Bampoulis, P.; Sotthewes, K.; Dollekamp, E.; Poelsema, B. Water confined in two-dimensions: Fundamentals and applications. _Surface Science Reports_ 2018, _73_ , 233–264
  * Shen et al. 2021 Shen, J.; Liu, G.; Han, Y.; et al., Artificial channels for confined mass transport at the sub-nanometre scale. _Nat. Rev. Mater._ 2021, _6_ , 294–312
  * Naushad 1989 Naushad, M. _A New Generation Material Graphene: Applications in Water Technology_ ; Springer, Cham, 1989; p 471
  * Wei et al. 2018 Wei, Y.; Zhang, Y.; Gao, X.; Ma, Z.; Wang, X.; Gao, C. Multilayered graphene oxide membranes for water treatment: A review. _Carbon_ 2018, _139_ , 964–981
  * Xu and Wang 2017 Xu, L.; Wang, J. The application of graphene-based materials for the removal of heavy metals and radionuclides from water and wastewater. _Critical Reviews in Environmental Science and Technology_ 2017, _47_ , 1042–1105
  * Geng et al. 2014 Geng, J.; Kim, K.; Zhang, J.; Escalada, A.; Tunuguntla, R.; Comolli, L. R.; Allen, F. I.; Shnyrova, A. V.; Cho, K.; Munoz, D.; Wang, Y. M.; Grigoropoulos, C. P.; Ajo-Franklin, C. M.; Frolov, V. A.; Noy, A. Stochastic transport through carbon nanotubes in lipid bilayers and live cell membranes. _Nature_ 2014, _514_ , 612–615
  * Valenga et al. 2020 Valenga, M. G. P.; Felsner, M. L.; de Matos, C. F.; de Castro, E. G.; Galli, A. Development and validation of voltammetric method for determination of amoxicillin in river water. _Analytica Chimica Acta_ 2020, _1138_ , 79–88
  * Thomas et al. 2014 Thomas, M.; Corry, B.; Hilder, T. A. What Have We Learnt About the Mechanisms of Rapid Water Transport, Ion Rejection and Selectivity in Nanopores from Molecular Simulation? _Small_ 2014, _10_ , 1453–1465
  * Bordin et al. 2012 Bordin, J. R.; Diehl, A.; Barbosa, M. C.; Levin, Y. Ion fluxes through nanopores and transmembrane channels. _Phys. Rev. E_ 2012, _85_ , 031914
  * Andreeva et al. 2021 Andreeva, D.; Trushin, M.; Nikitina, A.; et al, Two-dimensional adaptive membranes with programmable water and ionic channels. _Nat. Nanotechnol._ 2021, _16_ , 174–180
  * Park and Jung 2014 Park, H. G.; Jung, Y. Carbon nanofluidics of rapid water transport for energy applications. _Chem. Soc. Rev._ 2014, _43_ , 565–576
  * Liu et al. 2019 Liu, Z.; Zhang, H.; Eredia, M.; Qiu, H.; Baaziz, W.; Ersen, O.; Ciesielski, A.; Bonn, M.; Wang, H. I.; Samorì, P. Water-Dispersed High-Quality Graphene: A Green Solution for Efficient Energy Storage Applications. _ACS Nano_ 2019, _13_ , 9431–9441, PMID: 31386338
  * Olabi et al. 2021 Olabi, A.; Abdelkareem, M. A.; Wilberforce, T.; Sayed, E. T. Application of graphene in energy storage device – A review. _Renewable and Sustainable Energy Reviews_ 2021, _135_ , 110026
  * Fuentes-Azcatl and Domínguez 2019 Fuentes-Azcatl, R.; Domínguez, H. Carbon Dioxide Confined between Two Charged Single Layers of Graphene: Molecular Dynamics Studies. _The Journal of Physical Chemistry C_ 2019, _123_ , 23705–23710
  * Sugahara et al. 2019 Sugahara, A.; Ando, Y.; Kajiyama, S.; Yazawa, K.; Gotoh, K.; Otani, M.; Okubo, M.; Yamada, A. Negative dielectric constant of water confined in nanosheets. _Nature Communications_ 2019, _10_ , 850
  * Fumagalli et al. 2018 Fumagalli, L.; Esfandiar, A.; Fabregas, R.; Hu, S.; Ares, P.; Janardanan, A.; Yang, Q.; Radha, B.; Taniguchi, T.; Watanabe, K.; Gomila, G.; Novoselov, K. S.; Geim, A. K. Anomalously low dielectric constant of confined water. _Science_ 2018, _360_ , 1339–1342
  * Zhang et al. 2013 Zhang, C.; Gygi, F.; Galli, G. Strongly Anisotropic Dielectric Relaxation of Water at the Nanoscale. _The Journal of Physical Chemistry Letters_ 2013, _4_ , 2477–2481
  * Zhu et al. 2021 Zhu, H.; Hu, B.; Hu, H.; He, W.; Huang, J.; Li, G. Tunable dielectric constant of water confined in graphene oxide nanochannels. _Journal of Molecular Liquids_ 2021, _324_ , 115139
  * Motevaselian and Aluru 2020 Motevaselian, M. H.; Aluru, N. R. Confinement-Induced Enhancement of Parallel Dielectric Permittivity: Super Permittivity Under Extreme Confinement. _The Journal of Physical Chemistry Letters_ 2020, _11_ , 10532–10537, PMID: 33290076
  * Motevaselian and Aluru 2020 Motevaselian, M. H.; Aluru, N. R. Universal Reduction in Dielectric Response of Confined Fluids. _ACS Nano_ 2020, _14_ , 12761–12770, PMID: 32966055
  * Olivieri et al. 2021 Olivieri, J.-F.; Hynes, J. T.; Laage, D. Confined Water’s Dielectric Constant Reduction Is Due to the Surrounding Low Dielectric Media and Not to Interfacial Molecular Ordering. _The Journal of Physical Chemistry Letters_ 2021, _12_ , 4319–4326, PMID: 33914550
  * Köhler et al. 2017 Köhler, M. H.; Bordin, J. R.; da Silva, L. B.; Barbosa, M. C. Breakdown of the Stokes–Einstein water transport through narrow hydrophobic nanotubes. _Physical Chemistry Chemical Physics_ 2017, _19_ , 12921–12927
  * Celebi et al. 2019 Celebi, A. T.; Nguyen, C. T.; Hartkamp, R.; Beskok, A. The role of water models on the prediction of slip length of water in graphene nanochannels. _The Journal of Chemical Physics_ 2019, _151_ , 174705
  * Mondal and Bagchi 2020 Mondal, S.; Bagchi, B. Water Layer at Hydrophobic Surface: Electrically Dead but Dynamically Alive? _Nano Letters_ 2020, _20_ , 8959–8964, PMID: 33253589
  * Cerveny et al. 2016 Cerveny, S.; Mallamace, F.; Swenson, J.; Vogel, M.; Xu, L. Confined Water as Model of Supercooled Water. _Chemical Reviews_ 2016, _116_ , 7608–7625, PMID: 26940794
  * Sam et al. 2020 Sam, A.; Hartkamp, R.; Kannam, S. K.; Babu, J. S.; Sathian, S. P.; Daivis, P. J.; Todd, B. D. Fast transport of water in carbon nanotubes: a review of current accomplishments and challenges. _Molecular Simulation_ 2020, _0_ , 1–20
  * Xu 2020 Xu, H. Probing nanopore structure and confined fluid behavior in shale matrix: A review on small-angle neutron scattering studies. _International Journal of Coal Geology_ 2020, _217_ , 103325
  * Yue Hong 2021 Yue Hong, Q. L. Progress of Scanning Probe Microscopy in the Study of Confined Water in Two-Dimensional Materials. _Chinese Science Bulletin_ 2021, _66_ , 1689–1702
  * Gravelle 2016 Gravelle, S. https://hal.archives-ouvertes.fr/hal-02375018v1
  * Allen and Tildesley 1989 Allen, M. P.; Tildesley, D. J. _Computer simulation of liquids_ ; Clarendon Press, 1989; p 385
  * Underwood and Greenwell 2018 Underwood, T. R.; Greenwell, H. C. The Water-Alkane Interface at Various NaCl Salt Concentrations: A Molecular Dynamics Study of the Readily Available Force Fields. _Scientific Reports_ 2018, _8_ , 352
  * Fuentes-Azcatl and Alejandre 2014 Fuentes-Azcatl, R.; Alejandre, J. Non-Polarizable Force Field of Water Based on the Dielectric Constant: TIP4P/$\epsilon$. _The Journal of Physical Chemistry B_ 2014, _118_ , 1263–1272, PMID: 24422512
  * Fuentes-Azcatl et al. 2015 Fuentes-Azcatl, R.; Mendoza, N.; Alejandre, J. Improved SPC force field of water based on the dielectric constant: SPC/$\varepsilon$. _Physica A: Statistical Mechanics and its Applications_ 2015, _420_ , 116 – 123
  * Fuentes-Azcatl and Barbosa 2016 Fuentes-Azcatl, R.; Barbosa, M. C. Thermodynamic and dynamic anomalous behavior in the TIP4P/$\epsilon$ water model. _Physica A: Statistical Mechanics and its Applications_ 2016, _444_ , 86–94
  * Fuentes-Azcatl and Barbosa 2020 Fuentes-Azcatl, R.; Barbosa, M. C. Flexible bond and angle, FBA/$\varepsilon$ model of water. _Journal of Molecular Liquids_ 2020, _303_ , 112598
  * Fuentes-Azcatl 2021 Fuentes-Azcatl, R. Flexible model of water based on the dielectric and electromagnetic spectrum properties: TIP4P/$\epsilon_{Flex}$. _Journal of Molecular Liquids_ 2021, _338_ , 116770
  * Yu et al. 2010 Yu, H.; Whitfield, T. W.; Harder, E.; Lamoureux, G.; Vorobyov, I.; Anisimov, V. M.; MacKerell, A. D.; Roux, B. Simulating Monovalent and Divalent Ions in Aqueous Solution Using a Drude Polarizable Force Field. _Journal of Chemical Theory and Computation_ 2010, _6_ , 774–786, PMID: 20300554
  * Halgren and Damm 2001 Halgren, T. A.; Damm, W. Polarizable force fields. _Current Opinion in Structural Biology_ 2001, _11_ , 236 – 242
  * Dang et al. 1991 Dang, L. X.; Rice, J. E.; Caldwell, J.; Kollman, P. A. Ion solvation in polarizable water: molecular dynamics simulations. _Journal of the American Chemical Society_ 1991, _113_ , 2481–2486
  * Hasted 1972 Hasted, J. B. In _The Physics and Physical Chemistry of Water._ ; Franks, F., Ed.; Water (A Comprehensive Treatise); Springer, Boston, 1972
  * Tsimpanogiannis et al. 2019 Tsimpanogiannis, I. N.; Moultos, O. A.; Franco, L. F. M.; de M. Spera, M. B.; Erdős, M.; Economou, I. G. Self-diffusion coefficient of bulk and confined water: a critical review of classical molecular simulation studies. _Molecular Simulation_ 2019, _45_ , 425–453
  * Tsimpanogiannis et al. 2020 Tsimpanogiannis, I. N.; Jamali, S. H.; Economou, I. G.; Vlugt, T. J. H.; Moultos, O. A. On the validity of the Stokes–Einstein relation for various water force fields. _Molecular Physics_ 2020, _118_ , e1702729
  * Abal et al. 2020 Abal, J. P. K.; Bordin, J. R.; Barbosa, M. C. Salt parameterization can drastically affect the results from classical atomistic simulations of water desalination by MoS2 nanopores. _Phys. Chem. Chem. Phys._ 2020, _22_ , 11053–11061
  * Berendsen et al. 1987 Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The missing term in effective pair potentials. _The Journal of Physical Chemistry_ 1987, _91_ , 6269–6271
  * Fuentes-Azcatl and Barbosa 2016 Fuentes-Azcatl, R.; Barbosa, M. C. Sodium Chloride, NaCl/$\varepsilon$: New Force Field. _The Journal of Physical Chemistry B_ 2016, _120_ , 2460–2470, PMID: 26890321
  * Fuentes-Azcatl and Barbosa 2018 Fuentes-Azcatl, R.; Barbosa, M. C. Potassium bromide, KBr/$\epsilon$: New Force Field. _Physica A: Statistical Mechanics and its Applications_ 2018, _491_ , 480–489
  * Fuentes-Azcatl et al. 2021 Fuentes-Azcatl, R.; Araujo, G. J. C.; Lourenço, T. C.; Costa, C. T. O. G.; Carneiro, J. W. d. M.; Costa, L. T. Dielectric behavior of water in [bmim] [Tf2N] room-temperature ionic liquid: molecular dynamic study. _Theoretical Chemistry Accounts_ 2021, _140_ , 125
  * Fuentes-Azcatl and Domínguez 2019 Fuentes-Azcatl, R.; Domínguez, H. Prediction of experimental properties of CO2: improving actual force fields. _Journal of Molecular Modeling_ 2019, _25_
  * Bonthuis et al. 2011 Bonthuis, D. J.; Gekle, S.; Netz, R. R. Dielectric Profile of Interfacial Water and its Effect on Double-Layer Capacitance. _Phys. Rev. Lett._ 2011, _107_ , 166102
  * Shim et al. 2011 Shim, Y.; Jung, Y.; Kim, H. J. Graphene-Based Supercapacitors: A Computer Simulation Study. _The Journal of Physical Chemistry C_ 2011, _115_ , 23574–23583
  * Páll et al. 2015 Páll, S.; Abraham, M. J.; Kutzner, C.; Hess, B.; Lindahl, E. Tackling Exascale Software Challenges in Molecular Dynamics Simulations with GROMACS. Solving Software Challenges for Exascale. Cham, 2015; pp 3–27
  * Tuckerman et al. 2001 Tuckerman, M. E.; Liu, Y.; Ciccotti, G.; Martyna, G. J. Non-Hamiltonian molecular dynamics: Generalizing Hamiltonian phase space principles to non-Hamiltonian systems. _The Journal of Chemical Physics_ 2001, _115_ , 1678–1702
  * Essmann et al. 1995 Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. _The Journal of Chemical Physics_ 1995, _103_ , 8577–8593
  * Sobrino Fernández et al. 2016 Sobrino Fernández, M.; Peeters, F. M.; Neek-Amal, M. Electric-field-induced structural changes in water confined between two graphene layers. _Phys. Rev. B_ 2016, _94_ , 045436
  * Zokaie and Foroutan 2015 Zokaie, M.; Foroutan, M. Comparative study on confinement effects of graphene and graphene oxide on structure and dynamics of water. _RSC Adv._ 2015, _5_ , 39330–39341
  * Krott and Bordin 2013 Krott, L. B.; Bordin, J. R. Distinct dynamical and structural properties of a core-softened fluid when confined between fluctuating and fixed walls. _The Journal of Chemical Physics_ 2013, _139_ , 154502
  * Bordin et al. 2014 Bordin, J. R.; Krott, L. B.; Barbosa, M. C. High pressure induced phase transition and superdiffusion in anomalous fluid confined in flexible nanopores. _The Journal of Chemical Physics_ 2014, _141_ , 144502
  * Krott et al. 2015 Krott, L. B.; Bordin, J. R.; Barraz, N. M.; Barbosa, M. C. Effects of confinement on anomalies and phase transitions of core-softened fluids. _The Journal of Chemical Physics_ 2015, _142_ , 134502
  * Xia and Berkowitz 1995 Xia, X.; Berkowitz, M. L. Electric-Field Induced Restructuring of Water at a Platinum-Water Interface: A Molecular Dynamics Computer Simulation. _Phys. Rev. Lett._ 1995, _74_ , 3193–3196
  * Harrach and Drossel 2014 Harrach, M. F.; Drossel, B. Structure and dynamics of TIP3P, TIP4P, and TIP5P water near smooth and atomistic walls of different hydroaffinity. _The Journal of Chemical Physics_ 2014, _140_ , 174501
  * Celebi et al. 2018 Celebi, A. T.; Barisik, M.; Beskok, A. Surface charge-dependent transport of water in graphene nano-channels. _Microfluidics and Nanofluidics_ 2018, _22_ , 7
  * Zhao et al. 2020 Zhao, Z.; Zhou, R.; Sun, C. Molecular Dynamics Study of Water Diffusivity in Graphene Nanochannels. _International Journal of Thermophysics_ 2020, _41_ , 79
  * Kirkwood 1939 Kirkwood, J. G. The Dielectric Polarization of Polar Liquids. _The Journal of Chemical Physics_ 1939, _7_ , 911–919
  * Glättli et al. 2002 Glättli, A.; Daura, X.; van Gunsteren, W. F. Derivation of an improved simple point charge model for liquid water: SPC/A and SPC/L. _The Journal of Chemical Physics_ 2002, _116_ , 9811–9828
  * Fröhlich 1958 Fröhlich, H. _Theory of Dielectrics: Dielectric Constant and Dielectric Loss, 2nd ed._ ; Oxford University Press, 1958; p 385
  * Lee et al. 1984 Lee, C.; McCammon, J. A.; Rossky, P. J. The structure of liquid water at an extended hydrophobic surface. _The Journal of Chemical Physics_ 1984, _80_ , 4448–4455
  * Tocci et al. 2014 Tocci, G.; Joly, L.; Michaelides, A. Friction of Water on Graphene and Hexagonal Boron Nitride from Ab Initio Methods: Very Different Slippage Despite Very Similar Interface Structures. _Nano Letters_ 2014, _14_ , 6872–6877, PMID: 25394228
  * Sato et al. 2018 Sato, T.; Sasaki, T.; Ohnuki, J.; Umezawa, K.; Takano, M. Hydrophobic Surface Enhances Electrostatic Interaction in Water. _Phys. Rev. Lett._ 2018, _121_ , 206002
  * Varghese et al. 2019 Varghese, S.; Kannam, S. K.; Hansen, J. S.; P. Sathian, S. Effect of Hydrogen Bonds on the Dielectric Properties of Interfacial Water. _Langmuir_ 2019, _35_ , 8159–8166, PMID: 31121091

(a)  Figure 6: TOC Figure