firstofone

# Fourier coefficients and cuspidality of modular forms: a new approach

Soumya Das  Department of Mathematics
Indian Institute of Science
Bangalore – 560012, India
<EMAIL_ADDRESS>
###### Abstract.

We provide a simple and new induction based treatment of the problem of
distinguishing cusp forms from the growth of the Fourier coefficients of
modular forms. Our approach also gives the best possible ranges of the weights
for this problem, and has wide adaptability.

###### Key words and phrases:

Bounds on Fourier coefficients, Siegel modular forms, cuspidality, cusp forms

###### 2000 Mathematics Subject Classification:

Primary 11F30, 11F46, Secondary 11F50

##  1\. Introduction

The main objective of this paper is to provide a short, simple and new proof
of the fact that suitable growth of Fourier coefficients of modular forms (at
any given cusp) serve to distinguish between cusp forms and non-cusp forms.
Whereas this statement is not very difficult to prove for elliptic modular
forms (cf. [Koh10], [BD14], [BK19], or Proposition 3.1) it becomes more
interesting – and non-trivial for higher rank groups. Here we would be working
with Siegel’s modular group acting on Siegel’s upper half-space, but most of
the paper generalizes without much issue to other kinds of tube domains.

In [Koh10], Kohnen asked whether the Hecke bound distinguishes cusp forms – of
course this question is interesting only for congruence subgroups. Over time,
this question was investigated thoroughly, and more generally the following
conjecture was formulated. Essentially it has two parameters to deal with. For
the notation used below, see Section 2.

######  Conjecture 1.1.

Given a real number $c$, a weight $k_{0}$ and $F\in M^{n}_{k}(\Gamma)$
($\Gamma\subset\mathrm{Sp}_{n}(\mathbf{Z})$ is a congruence subgroup), suppose
that for all $k\geq k_{0}$, the bound $a_{F}(T)\ll_{F}\det(T)^{c}$ holds for
all $T\in\Lambda_{n}^{+}$ (see 2 for the definition) at some cusp and for some
$c<k-\frac{n+1}{2}$. Then $F\in S^{n}_{k}(\Gamma)$.

It is clear that without loss one can assume that the cusp is $\infty$. The
two parameters mentioned above are $c$ and $k_{0}$. Part of the question was
also to find the largest value of $c$ and the smallest $k_{0}$ such that for
all $k\geq k_{0}$ such that the above growth condition implies that $F$ is a
cusp form. For small weights, the situation can be rather delicate.

This question is quite natural since many times a modular form is defined via
its Fourier expansion. Usually the theoretical tool to check cuspidality is
the so-called Siegel’s $\Phi$-operator. One needs to check that $F$ is in the
kernel of such $\Phi$-operators at all the one-dimensional cusps, by which we
mean the elements of $\Gamma\backslash\mathrm{Sp}_{n}(\mathbf{Z})/C_{n,n-1}$,
where $C_{n,n-1}$ denotes the Klingen parabolic subgroup with last row
$0^{(1,2n-1)}$. The strategy to answer the above question directly by
decomposing the space in question in terms of Eisenstein series and cusp forms
quickly turns very complicated in higher degrees owing to the presence of
various Klingen-type Eisenstein series. Their Fourier expansion is mostly not
known in general, and from whatever is known, their Fourier coefficients grow
almost similar to each other. This renders the direct approach described above
as infeasible. The reader can refer to the introduction to [BK19, BD14, BD18]
more on this.

Thus all investigations on this problem in higher degrees hinged on finding
avenues avoiding the use of Klingen-Eisenstein series. To put things into
perspective, we briefly sketch the evolution of the techniques that were
employed on this problem, over the years.

– (i) In [Koh10, KM14, Lin14], covering the cases $n=1,2$, the direct approach
of using the newform theory and explicit Fourier expansions of (elliptic and
Jacobi respectively) Eisenstein series was used, focusing on Hecke’e bound.

– (ii) In [BD14] the problem was reduced to eigenforms and Andrianov’s
identity was used to obtain information about the poles of the standard
$L$-function $D_{F}(s)$ followed by Zahrkovskaya relations (which relate
$D_{F}(s)$ with $D_{\Phi F}(s)$) to reduce to cusp forms of lower degrees and
finally Shimura’s results on the analytic properties of $D_{F}(s)$. In fact
there was another idea of writing $F$ as linear combination of theta series
(basis problem of Eichler) for square-free levels and use special properties
of the theta series under the $\Phi$-operator. The cases covered (essentially)
were $M^{n}_{k}(\Gamma_{0}(N))$, $N$ square-free, $k>2n$, $c<k-n$.

– (iii) In [Miz15] a new idea, via Imai’s converse theorem (in degree $2$,
full level), was introduced into the problem. One proves that various Maaß
twists of the Koecher-Maaß series of $F$ have good analytic properties from
the hypothesis, and applies Imai’s converse theorem, focusing on Hecke’s
bound.

– (iv) In [BD16] special techniques involving ‘diagonal’ representatives of
$0,1$-dimensional cusps of $\Gamma^{(n)}_{0}(N)$ were employed, reducing the
problem to degrees $1$ or $2$ via the Witt-operators or the Fourier-Jacobi
expansion, and appealing to the results from (iv) in degrees $1,2$. The cases
covered (essentially) were $M^{n}_{k}(\Gamma_{0}(N))$, $k\geq\frac{n}{2}+1$,
$c<k-\frac{n+3}{2}$, $N\geq 1$; and also $k\geq\max\\{2,\frac{n-1}{2}\\}$,
$c<k-\frac{n+1}{2}$, $N=1$ or square-free.

– (v) In [BD18], the above procedure (cf. (ii)) was refined and a local
version of the approach was used, leading to the comparison of the sizes of
Hecke eigenvalues on a suitable set of primes. This allows one to ignore the
convergence properties of $D_{F}(s)$. The cases covered (essentially) were
$M_{k}(\Gamma)$, $k>2n$, $c<k-n$;

– (vi) In [BK19], a nice approach adapted from [Miy89], which transfers the
problem to the behaviour of $F(Z)$ as $\det(Y)\to 0$, was successfully carried
out. The approach works smoothly when the Fourier expansion of $F$ is
supported on positive definite matrices (at a cusp), but otherwise requires
some additional tricks to reduce to this situation. The cases covered
(essentially) were $M_{k}(\Gamma)$, $k>n+1$, $c<k-\frac{n+1}{2}$.

We prove the following result, whose proof is direct, and which improves upon
all of the the results described above in some aspect. Moreover, its proof is
different from all of the above (cf. (i) – (vi)), and it can be easily adapted
to various other situations.

######  Theorem 1.2.

Let $k\geq\frac{n}{2}$ be half-integral, $\Gamma$ any congruence subgroup, and
$F\in M^{n}_{k}(\Gamma)$. Suppose for any $T\in\Lambda_{n}^{+}$ and for some
$c<k-\frac{n+1}{2}$, one has the bounds $a_{F}(T)\ll\det(T)^{c}$. Then $F$ is
a cusp form.

The first and the main idea behind our proof is a suitable use of the Fourier-
Jacobi expansion. It is perhaps well-known that this serves as a very useful –
and effective – tool in tackling a variety of questions about Siegel modular
forms of higher degrees. In some sense our proof may be in retrospect somewhat
similar in spirit to the arguments in [BD22]; however the situations are
completely different.

It is clear (cf. [BD14, Remark 5.3]) that the supremum of the possible values
of $c$ in Theorem 1.2 is indeed $k-\frac{n+1}{2}$. The second crucial
observation is that the maximal exponent in the growth condition, viz.
$k-\frac{n+1}{2}$, is consistent along each Fourier-Jacobi expansion of $F$
(cf. Section 4). For instance, our inductive argument would not work if we
just request Hecke’s bound – it doesn’t fit into the inductive set-up – and it
works precisely because we have considered the all possible values of ‘$c$’ in
Theorem 1.2.

Now we talk about the quantity $k_{0}$ in Question 1.1. Whereas the large
values of $k_{0}$ in [BD14] occurred from the poles of $D_{F}(s)$, in [BK19]
it was partly due to the use of Lipschitz formula. Since our argument will
avoid all of these, we do not face any such issue. So our result is thus valid
for the minimum possible value of $k_{0}$, viz., $k_{0}=n/2$ – beyond which
there do not exist non-zero cusp forms, by the theory of singular weights (cf.
[Fre83]).

Another technical point is that in order for our argument to succeed (even
only for integral weights) we have to crucially work with half-integral
weights (and general congruence subgroups) from the beginning, unlike all of
the other works quoted above, where such weights were not considered as they
were not necessary in the proof. This is because via the Fourier-Jacobi
expansion, one encounters Jacobi forms whose theta components may have half-
integral weights. It is however noted in many places that some of the earlier
works generalize to half-integral weights. Another amusing feature is that we
need a statement which transfers cuspidality from the Fourier-Jacobi
coefficients of $F$ to $F$ itself. Whereas this was proven in [BD14, Theorem
5.7] as an ancillary result – for us it will be crucial.

We end the introduction by remarking that the method seemingly can not handle
vector-valued Siegel modular forms because there is no theta expansion
available for the Fourier-Jacobi coefficients in this setting. This is at
least known to the experts working in this area. One might try other
approaches, see subsection 5.4 for a discussion on this, and some other
related results, e.g. how the poles of the Rankin-Selberg $L$-series control
cuspidality (cf. subsection 5.2). We also discuss a refinement of 1.1 on
certain ‘thin-sets’, see subsection 5.1.

### Acknowledgements

The author thanks IISc. Bangalore, UGC Centre for Advanced Studies, DST India
for financial support.

##  2\. Notation and setting

We will use standard notation throughout the paper.

Let $M^{n}_{k}(\Gamma)$ denote the space of holomorphic Siegel cusp forms on a
congruence subgroup $\Gamma$ of scalar weight $k$. Every $F\in
M^{n}_{k}(\Gamma)$ can be expressed by a Fourier series of the form

(2.1) $F(Z)=\sum\nolimits_{T\in\Lambda_{n}}a_{F}(T)e(TZ/M);$

for some integer $M\geq 1$. Here $\Lambda_{n}$ denotes the set of half-
integral, positive semi-definite matrices; and $\Lambda_{n}^{+}$ its subset
consisting of positive-definite matrices. If $F$ is a cusp form, (2.1) is
supported on $\Lambda_{n}^{+}$ at all the cusps of $\Gamma$. We put
$\Gamma_{n}=\mathrm{Sp}_{n}(\mathbf{Z})$, and put
$M^{n}_{k}=M^{n}_{k}(\Gamma_{n})$, omit the superscript $n$ if $n=1$, etc.

For any $0\leq r\leq n-1$, $F\in M^{n}_{k}(\Gamma)$ has a Fourier-Jacobi
expansion for $F$ (of type $(n-r,r)$):

(2.2) $\displaystyle
F(Z)=\sum_{m\in\Lambda_{r}}\phi_{F,m}(\tau,z)e(m\tau^{\prime}/M)$

where $Z=\left(\begin{smallmatrix}\tau&z\\\
z^{t}&\tau^{\prime}\end{smallmatrix}\right)$, with
$\tau\in\mathbf{H}_{n-r},\tau^{\prime}\in\mathbf{H}_{r}$; $H_{n}$ being the
Siegel’s upper half-space of degree $n$.

We further let (with $A^{t}$ denoting the transpose of $A$)

(2.3) $\mathscr{U}(\Gamma)=\\{\left(\begin{smallmatrix}U&0\\\
0&U^{-t}\end{smallmatrix}\right)\in\Gamma\mid
U\in\mathrm{GL}_{n}(\mathbf{Z})\\}$

be the subgroup of the embedded copy of $\mathrm{GL}_{n}(\mathbf{Z})$ inside
$\Gamma$.

##  3\. Half-integral weights when $n=1$

We let $\kappa\in\tfrac{1}{2}\mathbf{Z}$ and $\Gamma$ be any congruence
subgroup of $\mathrm{SL}_{2}(\mathbf{Z})$ if $\kappa\in\mathbf{Z}$, or of
$\Gamma_{0}(4)$ if $\kappa\in\tfrac{1}{2}\mathbf{Z}\setminus\mathbf{Z}$,
unless stated otherwise.

In the Proposition below we prove our main result of the paper for $n=1$, for
all possible half-integral weights and all congruence subgroups. We use
Miyake’s characterisation of cusp forms (stated in [Miy89] for integral
weights), but also offer a different way to tackle the half-integral weights.
We could not find this result in the literature.

######  Proposition 3.1.

Let $\kappa\geq 1/2$ and $\mathfrak{f}\in M_{k}(\Gamma)$ be non-zero. If
$|a_{\mathfrak{f}}(n)|\ll_{\mathfrak{f}}n^{c}$ for any fixed $c<\kappa-1$,
then $\mathfrak{f}$ must be a cusp form.

###### Proof.

Our proof is based on the argument in [Miy89]. First of all from the Fourier
expansion of $\mathfrak{f}$, and the bound in the lemma, we find that (with
$\epsilon>0$ fixed, such that $c+\epsilon<\kappa-1$)

(3.1) $\displaystyle\mathfrak{f}(z)=a_{\mathfrak{f}}(0)+\sum\nolimits_{n\geq
1}n^{c}\exp(-\alpha\cdot ny)\ll|a_{\mathfrak{f}}(0)|+\sum\nolimits_{n\geq
1}n^{c}\cdot(ny)^{-(1+c+\epsilon)}$

for some constant $\alpha>0$ depending only on $\mathfrak{f}$. We avoided
using Lipschitz’s formula here.

We first assume that $c\geq-1$. Simply using that $e^{x}\gg x^{1+c+\epsilon}$
($x>0$), (3.1) shows that, as $y\to 0^{+}$, the bound

(3.2) $\displaystyle\mathfrak{f}(z)\ll y^{-(1+c+\epsilon)}.$

Note that $c$ may be negative, but we can assume without loss that $c\geq-1$
simply by using $n^{c}\leq n^{-1}$ in (3.1) otherwise. In this case one would
get $\displaystyle\mathfrak{f}(z)\ll y^{-\epsilon}$.

We now recall the statement necessary statement from Miyake:

Suppose $\mathfrak{f}\in M_{k}(\Gamma)$ be such that
$\mathfrak{f}(z)=O(y^{-\nu})$ as $y\to 0^{+}$ for some $0<\nu<k$. Then $f$ is
a cusp form.  (*)

Even though Miyake states the above result for $k$ integral, one can note that
(*) cf. [Miy89, Theorem 2.1.4] is valid for any real weight $k>0$ with
suitable modification of the setting. We however choose to invoke (*) only for
integral weights (clearly valid for $k\geq 1$), and reduce the case of
$1/2$-integral weights to it.

Thus when $\kappa\geq 1$ is integral, for all values of $c<\kappa-1$, (*)
implies that $\mathfrak{f}$ is cuspidal.

So let $\kappa=k+1/2\in\frac{1}{2}\mathbf{Z}\setminus\mathbf{Z}$ and without
loss assume that $\mathfrak{f}\in M_{k}(\Gamma(N))$. Consider the modular form
(recall that $4|N$ in this case)

$\mathfrak{g}(\tau):=\mathfrak{f}(4N\tau)\cdot\theta(\tau)\in
M_{k+1}(\Gamma_{1}(16N^{2})),$

where $\theta(\tau)=\sum_{n\in\mathbf{Z}}q^{n^{2}}$. Thus if
$\mathfrak{f}(\tau)=\sum_{n\geq 1}a_{\mathfrak{f}}(n)q^{n/4N}$, then from the
definition of $\mathfrak{g}$,

$a_{\mathfrak{g}}(n)=\sum\nolimits_{\ell+m^{2}=n}a_{\mathfrak{f}}(\ell).$

Since $\\{(\ell,m)|\ell+m^{2}=n\\}\ll n^{1/2}$, by using the bound
$|a_{\mathfrak{f}}(n)|\ll_{f}n^{c}$ we easily get that

$|a_{\mathfrak{g}}(n)|\ll_{\mathfrak{g}}n^{c+1/2}.$

That $\mathfrak{g}$ is a cusp form now follows from the integral weight case,
noting that $k+1\geq 1$. We claim this implies that
${\mathfrak{f}}_{4N}(\tau):=\mathfrak{f}(4N\tau)\in
M_{k+1/2}(\Gamma_{1}(16N^{2}))$ is also a cusp form. It is enough to show that
the order of vanishing $\mathrm{ord}_{s}({\mathfrak{f}}_{4N})$ at any cusp $s$
of $\Gamma_{1}(16N^{2})$ is positive.

To see this, first note that at the cusps $\infty,0,1/2$ of $\Gamma_{0}(4)$,
the theta function looks like $\sum_{n\in\mathbf{Z}}q^{n^{2}}$,
$\sum_{n\in\mathbf{Z}}q^{n^{2}/4}$, $q^{1/4}\sum_{n\in\mathbf{Z}}q^{n^{2}+n}$
respectively, possibly up to some non-zero constants. Since any cusp of
$\Gamma_{1}(16N^{2})$ is $\Gamma_{0}(4)$ equivalent to one of the cusps of
$\Gamma_{0}(4)$, and $\theta$ is automorphic under $\Gamma_{0}(4)$, we can
write

$1\leq\mathrm{ord}_{s}(\mathfrak{g})=\mathrm{ord}_{s}({\mathfrak{f}}_{4N})+\mathrm{ord}_{l(s)}(\theta),$

where $l(s)\in\\{\infty,0,1/2\\}$. When $l(s)=\infty,0$, we see that
$\mathrm{ord}_{l(s)}(\theta)=0$ and the assertion about $\mathfrak{f}_{4N}$ is
clear. When $l(s)=1/2$, we see that

$\mathrm{ord}_{s}({\mathfrak{f}}_{4N})\geq 1-\mathrm{ord}_{1/2}(\theta)=3/4.$

This implies that $\mathfrak{f}_{4N}$ is cuspidal and hence so is
$\mathfrak{f}$. The lemma is proved. ∎

##  4\. Higher degrees

We would now want to generalise the previous results to higher degrees. The
main tool will be the Fourier-Jacobi expansion along with induction on the
degree. As remarked in the introduction, the condition $c<k-\frac{n+1}{2}$
will fit into the inductive framework naturally, yielding Theorem 1.2.

###### Proof of Theorem 1.2.

We induct on $n$. For $n=1$, this is Lemma 3.1. Next, we assume the theorem
for all congruence subgroups in degree $n-1$. Our idea is to pass to degree
$n$ by using the Fourier-Jacobi expansion for $F$ (of type $(n-1,1)$) as given
below. We have

(4.1) $\displaystyle F(Z)=\sum\nolimits_{m\geq
0}\phi_{m}(\tau,z)e(m\tau^{\prime}),$

where
$\tau\in\mathbf{H},z\in\mathbf{C}^{n-1},\tau^{\prime}\in\mathbf{H}_{n-1}$. The
theta decomposition for each $\phi_{m}$ is denoted by

(4.2)
$\displaystyle\phi_{m}(\tau,z)=\sum\nolimits_{\mu\in\mathbf{Z}^{n-1}/(2m\,\mathbf{Z}^{n-1})}h_{\mu}(\tau)\theta_{m,\mu}(\tau,z).$

Then it is known (cf. [Zie89]) that $h_{\mu}\in
M_{k-\frac{1}{2}}(\Gamma^{\prime})$ for some congruence subgroup
$\Gamma^{\prime}$. Note that the theta series have weight $\frac{1}{2}$ in
this case. The Fourier expansion of $h_{\mu}$ is given by

(4.3) $\displaystyle
h_{\mu}(\tau)=e(-\frac{1}{4}m^{-1}[\mu]\tau)\sum_{l\in\Lambda_{n-1},\,l\geq\frac{1}{4}m^{-1}[\mu^{t}]}c(l,\mu)e(l\tau),$

where

(4.4) $\displaystyle c(l,\mu):=a_{F}(\begin{pmatrix}l&\mu/2\\\
\mu^{t}/2&m\end{pmatrix}).$

We next would prove that for all $m\geq 1$, $\phi_{m}$ are cuspidal. An
application of [BD14, Theorem 5.7] with $\mathcal{S}=\\{0\\}$ would then show
that $F$ is indeed cuspidal, finishing the proof.

To prove that $\phi_{m}$ is cuspidal, we can consider $m$ to be fixed. We
would prove that for each $\mu$, $h_{\mu}$ is cuspidal. This being granted,
looking at the Fourier expansion of $\phi_{m}$, it immediately follows that
$\phi_{m}$ is cuspidal.

It remains to deal with the $h_{\mu}$’s. Note that $\mu$ is also fixed here.
We consider $\displaystyle g(\tau):=h_{\mu}(4m\tau)$ so that its Fourier
expansion may be written as

(4.5) $\displaystyle
g(\tau)=\sum_{L\in\mathbf{\Lambda}_{n-1},\,L\equiv-\mu\mu^{t}\bmod{4m}}c(\frac{L+\mu\mu^{t}}{4m},\mu)e(L\tau)$

We look at the positive definite Fourier coefficients of $g$. Upon using the
Jacobi decomposition for positive definite matrices, from (4.4) and (4.5) we
can write by our hypothesis on $F$ that,

(4.6) $\displaystyle
a_{g}(L)\ll\det\left(\begin{pmatrix}\frac{L+\mu\mu^{t}}{4m}&\mu/2\\\
\mu^{t}/2&m\end{pmatrix}\right)^{c}=m^{c}\det(L)^{c}\ll\det(L)^{c}.$

since $m$ is fixed, and $c$ is as in the theorem.

Now notice that $c<k-\frac{n+1}{2}=(k-\frac{1}{2})-\frac{n}{2}$. Since the
weight of $g$ is $k-\frac{1}{2}\geq\frac{n-1}{2}$ (cf. [Zie89, §3]), and $g$
is of degree $n-1$ on $\Gamma(4mN)$ ($N=$ level of $\Gamma$, cf. [BD22, §4,
Part B]), we can therefore invoke the induction hypothesis on $g$ to conclude
that all the $h_{\mu}$ are cusp forms, as desired. ∎

######  Remark 4.1.

We note that one could have used any of the Fourier-Jacobi expansions of $F$
(w.r.t.) the different parabolic subgroups, and the same argument as above
would have worked, showing that $\phi_{M}$ to be cuspidal for all $M>0$ of
size $1\leq r\leq n$. This is because of the identity

$k-\frac{n+1}{2}=k-\frac{r}{2}-\frac{n-r+1}{2}=k_{r}-\frac{(n-r)+1}{2};$

where $\phi_{M}$ has weight $k_{r}=k-\frac{r}{2}$ and its theta components
have degree $n-r$ (cf. [Zie89, §3]). However, we could not have invoked the
[BD14, Theorem 5.7], since the it requires that for all but finitely many
$M\in\mathrm{GL}_{r}(\mathbf{Z})\backslash\Lambda_{r}$, $\phi_{M}$ be
cuspidal. We do not have this information, as we consider only positive-
definite matrices.

######  Remark 4.2.

We note that the growth of the Fourier coefficients and that of $F$ as
$\det(Y)\to 0$ are actually equivalent. One way is the content of [BK19] (cf.
[Miy89]), and the other way follows immediately from the formula of
$\displaystyle a_{F}(T)=\int_{X\bmod 1}F(Z)e(-TZ)dX$, where $Y=T^{-1}$. This
suggests a different approach towards bounding Fourier coefficients of cusp
forms. In this aspect, one has a (weak) result from [Das24, Corollary 5.2].

##  5\. Further remarks

###  5.1. Question 1.1 on ‘thin sets’

In this subsection we want to discuss Question 1.1 on ‘thin sets’ – which are
in some sense "minimal" sets
$S\subset\mathrm{GL}_{n}(\mathbf{Z})\backslash\Lambda^{+}_{n}$ on which the
bound on the Fourier coefficients will guarantee cuspidality. Observe that in
almost all previously known results, say in level $1$, the assumption
$S=\mathrm{GL}_{n}(\mathbf{Z})\backslash\Lambda^{+}_{n}$ is present – and
crucially used in various ways – e.g. to invoke Lipschitz’s formula or
application of the Siegel’s $\Phi$ operator. When $S$ is significantly "thin",
the problem clearly becomes very hard.

– Example: As an example, in the extreme case that $S$ is just allowed to be
infinite, we get counter-examples, as follows. By a result of Kitaoka [Kit86,
Theorem 1.4.5, p. 73], if $G\in M^{n}_{k}(\Gamma)$ be such that its constant
terms are all zero at all the zero dimensional cusps of $\Gamma$, then say for
scalar $T=mI_{n}$ ($m\to\infty$) one actually has the bound

(5.1) $\displaystyle a_{G}(T)\ll
m^{n(k-(n+1)/2)-k+n+1}=\det(T)^{k-(n+1)/2-k/n+(n+1)/n}.$

The above bound is even better than what is needed for Question 1.1, yet $G$
is not cuspidal.

######  Conjecture 5.1.

A positive answer to Question 1.1 holds, provided the growth condition is
given on a set $S\subset\mathscr{U}(\Gamma)\backslash\Lambda^{+}_{n}$ with
positive upper density – i.e.,

$\delta_{+}(S):=\limsup_{X\to\infty}\frac{\\#\\{T\in S\mid\det(T)\leq
X\\}}{\\#\\{T\in\mathscr{U}(\Gamma)\backslash\Lambda_{n}^{+}\mid\det(T)\leq
X\\}}>0.$

###  5.2. Information from poles

We prove some results relating the cuspidality of a modular form with the
analytic behaviour of its Rankin-Selberg $L$-series. For a cusp $\omega_{j}$
of $\Gamma$, put

$R_{j}(f,f;s):=R(f|_{\omega_{j}},f|_{\omega_{j}};s)=\sum\nolimits_{n\geq
1}\frac{|a_{j}(n)|^{2}}{n^{s}},$

where $a_{j}(n)$ are the Fourier coefficients of $f|_{\omega_{j}}$. Similarly,
for $F\in M^{n}_{k}$, for $\mathrm{Re}(s)$ large enough, define the Rankin-
Selberg $L$-series of $F$ by

(5.2)
$\mathcal{R}(F,F;s):=\sum_{T\in\mathrm{GL}_{n}(\mathbf{Z})\backslash\Lambda_{n}^{+}}\frac{|a_{F}(T)|^{2}}{\epsilon(T)\,\det(T)^{s}},$

where $\epsilon(T)$ denotes the number of units of $T$.

######  Proposition 5.2.

Let $R_{j}(f,f;s)$ denote the Rankin-Selberg $L$-series for $f$ at the cusp
$\omega_{j}$ of $\Gamma$. Then $f\in M_{k}(\Gamma)$ is a cusp form if and only
if none of the $R_{j}(f,f;s)$ has a pole at $2k-1$.

For instance for the level $1$ Eisenstein series $E_{k}$, it is well known
that (see e.g. [BC08, p. 821] or [Zag81])

(5.3) $\displaystyle
R(E_{k},E_{k};s)=\frac{\zeta(s)\zeta(s-k+1)^{2}\zeta(s-2k+2)}{\zeta(2s-2k+2)},$

from which the pole at $s=2k-1$ is visible.

###### Proof.

For the proof one considers the $\Gamma$-invariant function
$g(z)=y^{k}|f(z)|^{2}$ and applies the theory of Rankin-Selberg series for
functions which are not of rapid decay via renormalization techniques. In
particular from the approach of truncation of domain (see [Dut97]) it is clear
that for a non-cusp form one will have at some cusp

(5.4) $\displaystyle g(z)=\psi(y)+O(y^{-A})\quad\quad\text{for all }A>0\text{
as }y\to\infty,$

with $\psi(y)=|a_{f}(0)|^{2}y^{k}$ and $a_{f}(0)\neq 0$; and for a cusp form
$\psi(y)=0$ at all cusps. The result now follows directly from [Dut97,
THEOREM] – in fact the residue at $s=2k-1$ is proportional to
$|a_{f}(0)|^{2}$. Note that our $R(f,f;s)$ differs from that of [Dut97] by the
transformation $s\to s-k+1$. ∎

The following corollary is immediate from [CN62], if we notice that $2k-1$ is
the rightmost pole of $R_{j}(f,f;s)$, where $\omega_{j}$ is a cusp of $\Gamma$
where $f$ has a non-zero constant term.

######  Corollary 5.3.

Let $f\in M_{k}(\Gamma)$ be non-zero and that $\kappa\geq 3/2$. Then, $f$ is
not a cusp form if and only if at some cusp $\omega_{j}$ (see above), one has

(5.5) $\sum\nolimits_{n\leq
x}|a_{j}(n)|^{2}=c_{\mathfrak{f},\Gamma}x^{2\kappa-1}+O_{\mathfrak{f},\Gamma}(x^{\eta})$

for some constant $c_{\mathfrak{f},\Gamma}>0$ not depending on $x$, and some
$\eta<2\kappa-1$.

It is natural to wonder what happens in the higher degrees. Following [BC08],
for $s\in\mathbf{C}$ we put (for a certain ‘growth-killing’ differential
operator $\mathcal{D}_{n}$)

(5.6) $\mathcal{R}^{*}(F,F;s):=\langle
E^{n}(s+\frac{n+1}{2}-k,\cdot),\,\mathcal{D}_{n}(\det(Y)^{k}|F|^{2})\rangle.$

so that
$\mathcal{R}^{*}(F,F;s)=c\Gamma_{n}(s)\phi_{n}(s)\phi_{n}(s+n-2k)\mathcal{R}(F,F;s)$,
where $\Gamma_{n}(s)$ is the generalized Gamma function, and $\phi_{n}(T)$ is
a polynomial defined by $\phi_{n}(T)=\prod_{j=0}^{n-1}(T-j/2)$. Here $c$ is a
function of $s$, which is not important. We observe, from (5.6), it is clear
that has possible poles at the points

(5.7) $s_{j}:=2k-\frac{n+j}{2},\quad\quad(1\leq j\leq n).$

At least if $k$ is suitably large, note that none of this poles come from the
Eisenstein series, its poles can only be at $j/4$, $0\leq j\leq 2n+2$.

######  Proposition 5.4.

Let $k>n+1$. Then $F\in M^{n}_{k}$ is a cusp form if and only if
$\mathcal{R}(F,F;s)$ does not have a pole at any of the points $s_{j}$.

For congruence subgroups, probably the above condition on the poles at all the
cusps should be enough, but one has to be more careful. Also it would be
interesting to know whether for each point $s_{j}$, there exists an Eisenstein
series (or some other modular form) whose Rankin-Selberg series has a pole at
$s_{j}$.

###### Proof.

One way is clear because for a cusp form $F$, we know that
$\mathcal{R}(F,F;s)$ has no poles at any of the $s_{j}$. This follows from the
integral representation of $\mathcal{R}(F,F;s)$, as can be read from [Kal84,
Theorem 2]. Namely, in this case the poles can at most be at

(5.8) $t_{j}:=k-j/4,\quad\quad 0\leq j\leq 2n+2,$

coming from the Eisenstein series (with a simple pole at $t_{0}=k$), and at
most at

(5.9) $\rho/2+k-(n+1)/2\quad\text{ and }\quad\rho/4+k/2-(n+1-j)/2,\quad\quad
1\leq j\leq\lfloor n/2\rfloor,$

where $\rho$ varies over the set of non-trivial zeros of the Riemann zeta
function. From this it is easy to check that none of points $s_{j}$ from (5.7)
can occur among these.

Conversely, let $F\in M^{n}_{k}$ be such that $\mathcal{R}(F,F;s)$ has no
poles at any of the $s_{j}$. Then it may have poles at most at one of the
points listed in (5.8) and (5.9). The largest of these is $s=k$. By [CN62],
then we would definitely be able to conclude that
$a_{F}(T)\ll_{F}\det(T)^{k/2}$ – for all $T>0$. Notice that the Dirichlet
series (5.2) defining $\mathcal{R}(F,F;s)$ would now converge in the region
$\mathrm{Re}(s)>k$. Then Theorem 1.2 implies that $F$ must be a cusp form,
since $\frac{k}{2}<k-\frac{n+1}{2}$. ∎

###  5.3. Some variants

In analogy with the usual Fourier expansion, one can ask the following
question based on the (extended) Petersson pairing of the Fourier-Jacobi
coefficients.

######  Question 5.5.

Let $F\in M^{n}_{k}$ be such that its Fourier Jacobi coefficients $\phi_{F,m}$
(with $m\in\Lambda_{r}^{+}$, $0\leq r\leq n-1$) satisfy

(5.10)
$\displaystyle|\langle\phi_{F,m},\phi_{F,m}\rangle_{reg}|\ll_{F}m^{\delta}$

for some $\delta<k-\frac{n-r+1}{2}$. Then $F$ is a cusp form.

Here $\langle\,,\,\rangle_{reg}$ denotes the regularized Petersson pairing of
$\phi_{F,m}$, which may not be positive definite. See [BD18a] for the case
$n=2,r=1$ – from there the regularization is known and also the bound
$\displaystyle\langle\phi_{F,m},\phi_{F,m}\rangle_{reg}\ll_{F}m^{k}$ – i.e.,
the Hecke bound holds! Whereas for cusp form $G$ of degree $2$, one expects
that $\displaystyle\langle\phi_{G,m},\phi_{G,m}\rangle\ll_{F}m^{k-1}$ for all
$m\geq 1$, cf. [Koh11].

###  5.4. Vector-valued Siegel modular forms

As we have remarked in the introduction, the lack of a good theta expansion of
Jacobi forms makes our approach not immediately useful. We however sketch
another possible approach based on restriction of domains – hopefully it can
be made to work in the future.

Let $F\in S_{\rho}(\Gamma_{n})$ where $\rho$ is a polynomial representation of
$\mathrm{GL}_{n}(\mathbf{C})$. Following Freitag [Fre79], we put
$h(z):=\rho(A)\cdot F(A^{t}z^{*}A)$, where $z=(z_{1},z_{2},\ldots,z_{n})$,
$z^{*}=\mathrm{diag}(z_{1},z_{2},\ldots,z_{n})$, and
$A\in\mathrm{GL}_{n}(\mathbf{R})$ chosen in such a way that $\Gamma_{A}$
($\subset\mathrm{SL}_{2}(\mathbf{R})^{n}$, a certain subgroup, see [Fre79] for
the precise definition) is commensurable with a Hilbert modular group over
some totally real field $K$ and moreover that none of the components of $h(z)$
is constant. Then it is known from [Fre79] that all of the components of $h$
is a Hilbert modular form with respect to some weights of $\rho$.

The growth condition (say Hecke’s bound) for $\rho$ reads as
$\left\lVert\rho(T^{-1/2})a_{F}(T)\right\rVert\ll 1$. This should imply a
similar growth condition on the Fourier coefficients of $h$, which should
imply that each component of $h$, and hence $h$ is cuspidal. Here only
parallel weights play a role. Since the set of $A$ described above is dense in
$\mathrm{GL}_{n}(\mathbf{R})$, one should be able to prove that $\Phi(F)=0$.
The main point is the transition of the growth condition.

## References

  * [BC08] Siegfried Böcherer and F.. Chiera  “On Dirichlet series and Petersson products for Siegel modular forms”  In _Annales de L’Institut Fourier_ 58.3, 2008, pp. 801–824
  * [BD14] Siegfried Böcherer and Soumya Das  “Characterization of Siegel cusp forms by the growth of their Fourier coefficients”  In _Math. Ann._ 359.1–2 Springer, 2014, pp. 169–188
  * [BD16] Siegfried Böcherer and Soumya Das  “Cuspidality and the growth of Fourier coefficients: small weights”  In _Mathematische Zeitschrift_ 283 Springer, 2016, pp. 539–553
  * [BD18] Siegfried Böcherer and Soumya Das  “Cuspidality and the growth of Fourier coefficients of modular forms”  In _J. Reine Angew. Math_ 741 De Gruyter, 2018, pp. 161–178
  * [BD18a] Siegfried Böcherer and Soumya Das  “Petesson norms of not necessarily cuspidal Jacobi modular forms and applications”  In _Adv. Math._ 336 Elsevier, 2018, pp. 335–376
  * [BD22] Siegfried Böcherer and Soumya Das  “On fundamental Fourier coefficients of Siegel modular forms”  In _Journal of the Institute of Mathematics of Jussieu_ 21.6 Cambridge University Press, 2022, pp. 2001–2041
  * [BK19] S. Böcherer and W. Kohnen  “On the Fourier coefficients of Siegel Modular Forms”  In _Nagoya Math. J._ 234, 2019, pp. 1–16
  * [CN62] K. Chandrasekharan and R. Narasimhan  “Functional equations with multiple gamma factors and the average order of arithmetical functions”  In _Ann. of Math.(2)_ 76, 1962, pp. 93–136
  * [Das24] Soumya Das  “$L^{\infty}$-sizes of the spaces Siegel cusp forms of degree $n$ via Poincaré series”  In _arXiv preprint arXiv:2406.19335v1_ , 2024
  * [Dut97] Shamita Dutta  “On the Rankin-Selberg method for functions not of rapid decay on congruence subgroups”  In _Journal of Number Theory_ 62 Elsevier, 1997, pp. 115–126
  * [Fre79] Eberhard Freitag  “Ein Verschwindungssatz für Automorphe Formen zur Siegelschen Modulgruppe”  In _Math. Z._ 165, 1979, pp. 11–18
  * [Fre83] Eberhard Freitag  “Siegelsche modulfunktionen”  Springer-Verlag, 1983
  * [Kal84] V.. Kalinin  “Analytic properties of the convolution of Siegel modular forms of genus $n$”  In _Math. USSR Sb._ 48.1, 1984, pp. 193–200
  * [Kit86] Y. Kitaoka  “Lectures on Siegel modular forms and representation by quadratic forms”  Springer-Verlag, Berlin for T.I.F.R., 1986
  * [KM14] W. Kohnen and Y. Martin  “A characterization of degree $2$ cusp forms by the growth of their Fourier coefficients”  In _Forum Math._ 26.5 De Gruyter, 2014, pp. 1323–1331
  * [Koh10] W. Kohnen  “On certain generalized modular forms”  In _Funct. Approx. Comment. Math._ 43.5, 2010, pp. 23–29
  * [Koh11] W. Kohnen  “On the growth of Petersson norms of Fourier-Jacobi coefficients of Siegel cusp forms”  In _Bull. London Math. Soc._ 43, 2011, pp. 717–720
  * [Lin14] B. Linowitz  “Characterizing adelic Hilbert modular cusp forms by coefficient size”  In _Kyushu J. Math_ 68.1 Springer, 2014, pp. 105–111
  * [Miy89] T. Miyake  “Modular Forms”  Springer, Berlin, Heidelberg, 1989
  * [Miz15] Yoshinori Mizuno  “On characterisation of Siegel cusp forms by the Hecke bound”  In _Mathematika_ 61 Springer, 2015, pp. 89–100
  * [Zag81] Don Zagier  “The Rankin-Selberg method for automorphic functions which are not of rapid decay”  In _J. Fac. Sci. Univ. Tokyo Sect. IA Math._ 28, 1981, pp. 415–437
  * [Zie89] C. Ziegler  “Jacobi Forms of higher degree”  In _Abh. Math. Sem. Univ. Hamburg_ 59 Springer, 1989, pp. 191–224


# A theory for the stabilization of polar crystal surfaces by a liquid
environment

Stephen J. Cox  Yusuf Hamied Department of Chemistry, University of Cambridge,
Lensfield Road, Cambridge CB2 1EW, United Kingdom<EMAIL_ADDRESS>
###### Abstract

Polar crystal surfaces play an important role in the functionality of many
materials, and have been studied extensively over many decades. In this
article, a theoretical framework is presented that extends existing theories
by placing the surrounding solution environment on an equal footing with the
crystal itself; this is advantageous, e.g., when considering processes such as
crystal growth from solution. By considering the polar crystal as a stack of
parallel plate capacitors immersed in a solution environment, the equilibrium
adsorbed surface charge density is derived by minimizing the free energy of
the system. In analogy to the well-known diverging surface energy of a polar
crystal surface at zero temperature, for a crystal in solution it is shown
that the “polar catastrophe” manifests as a diverging free energy cost to
perturb the system from equilibrium. Going further than existing theories, the
present formulation predicts that fluctuations in the adsorbed surface charge
density become increasingly suppressed with increasing crystal thickness. We
also show how, in the slab geometry often employed in both theoretical and
computational studies of interfaces, an electric displacement field emerges as
an electrostatic boundary condition, the origins of which are rooted in the
slab geometry itself, rather than the use of periodic boundary conditions.
This aspect of the work provides a firmer theoretical basis for the recent
observation that standard “slab corrections” fail to correctly describe, even
qualitatively, polar crystal surfaces in solution.

##  I Introduction

Crystal morphology is a determining factor in the mechanical, rheological and
catalytic properties of ionic crystals Dandekar, Kuvadia, and Doherty (2013).
As such, there is considerable effort to establish simple and inexpensive
routes to control which facets of a crystal are exposed at its surfaces. In
this article, we will concern ourselves exclusively with polar crystal faces,
which we describe in more detail in Sec. II.1. Examples of polar crystal
surfaces include halite (111), wurtzite (0001), and zincblende (111). Polar
surfaces have been studied extensively from a solid state physics perspective
Tasker (1979); Noguera (2000); Goniakowski, Finocchi, and Noguera (2008),
owing to their use in microelectronic devices and the now well-established
fabrication techniques of thin metal oxide films. Polar surfaces are also
important for catalysis. For example, the 111 facets of MgO exhibit ultrahigh
catalytic activity Zhu  _et al._ (2006) for Claisen-Schmidt condensation of
benzaldehyde and acetophone, and the (0001) face of ZnO is active in methanol
Cheng, Akhter, and Kung (1983) decomposition.

What makes polar surfaces both challenging and interesting is that basic
electrostatic arguments show that their surface energy diverges with
increasing crystal thickness along a polar crystallographic direction Nosker,
Mark, and Levine (1970); Tasker (1979); Noguera (2000); this is the so-called
“polar catastrophe.” The fact that polar crystal surfaces are observed thus
implies a “polarity compensation mechanism” that overcomes the underlying
instability. Known polarity compensation mechanisms are Noguera (2000);
Goniakowski, Finocchi, and Noguera (2008): (i) charge transfer via
nonstoichiometric reconstruction, which amounts to an effective transfer of
ions from one side of the crystal to the other; (ii) electronic
reconstruction, where interfacial charges are modified by partial filling of
electronic interface states; and (iii) adsorption of foreign atoms or ions
that provide an appropriate amount of compensating interfacial charge. The
theoretical rationalization of the polar catastrophe from energetic
considerations of a crystal surface in contact with a low-density vapor phase
lends itself most naturally to describing mechanisms (i) and (ii).

Polarity compensation through adsorption of ions from a surrounding liquid
environment can be considered an extreme example of mechanism (iii), and
provides a potentially simple and cost-effective route to controlled exposure
of polar crystal facets. For example, using the molten salt synthesis route,
Xu _et al._ Xu  _et al._ (2009) have synthesized particles of ZnO, MgO and
Co3O4 that exclusively expose polar facets, and subsequent studies from Susman
_et al._ suggest that K+ and Cl- ions potentially stabilize the 111 facets of
MgO Susman  _et al._ (2018), while a complex interplay between the molten salt
and resulting crystal morphology is evident for NiO Susman  _et al._ (2020).
To understand the processes that underlie polarity compensation by ion
adsorption from a surrounding liquid, it will be useful to have a theory that
treats the environment on an equal footing to the crystal. The purpose of this
article is to develop such a theoretical framework.

A second motivation is to reiterate the central message of Tasker’s seminal
article Tasker (1979) on the theory of polar surfaces: to demonstrate that the
polar catastrophe has a true physical origin and is not the result of lattice
summation techniques. Such a (re)clarification is necessary as the use of
molecular simulations—and consequently, lattice summation techniques—has
exploded since Tasker’s original 1979 publication, and the technical nature of
treating long-range electrostatic interactions under periodic boundary
conditions can obfuscate physical intuition. Similar to Tasker, then, the
theory presented here will not rely on periodic boundary conditions.
Nonetheless, we will see implications for simulations that do employ periodic
boundary conditions, and we make two observations that are useful from a
practical viewpoint when modelling, with a slab geometry, the surfaces of
polar crystals in contact with a conducting medium such as an electrolyte:

  1. 1.

tinfoil Ewald approaches, briefly described in Sec. II.2, are appropriate if
we are interested in genuinely thin polar crystals.

  2. 2.

An electric displacement field emerges as a boundary condition.

The second point clarifies the role of the electric displacement field imposed
in previous simulation studies Sayer, Zhang, and Sprik (2017); Sayer, Sprik,
and Zhang (2019); Sayer and Cox (2020). Specifically, the current work
emphasizes that it is conceptually incorrect to associate this electric
displacement field with the removal of spurious interactions between periodic
replicas of the simulation cell, and places the observation made in Ref. Sayer
and Cox, 2020—that standard “slab corrections” lead to qualitatively incorrect
results—on a firmer theoretical footing.

The rest of the article is organized as follows. In Sec. II.1 we derive an
expression for the equilibrium adsorbed surface charge density at halite (111)
in contact with an electrolyte solution, as a function of crystal thickness.
We also show how the polar catastrophe manifests as a diverging free energy
cost for fluctuations away from equilibrium as the crystal thickness
increases. In Sec. II.2, we show that the theoretical predictions are
consistent with results from molecular dynamics simulations that employ a
tinfoil Ewald approach to treat electrostatic interactions. In Sec. III, we
discuss the emergence of the electric displacement field in the slab geometry,
and also outline a justification for the use of tinfoil Ewald approaches to
model genuinely thin polar crystals. We summarize our findings in Sec. IV,
before providing a brief summary of the simulation methods in Sec. V. We also
provide an appendix detailing a more formal justification for the use of
tinfoil Ewald approaches.

##  II The polar catastrophe from free energy minimization

###  II.1 Derivation

A schematic the system we are interested is shown in Fig. 1a. Here, a crystal
exposing polar facets is immersed in an electrolyte environment. We assume
that the crystal assumes a bulk truncated structure, i.e., nonstoichiometric
reconstruction [mechanism (i)] has not taken place. Moreover, we assume an
ionic model such that electronic reconstruction [mechanism (ii)] is not
possible. With these approximations, the exposed crystal faces have surface
charge density $\pm\sigma_{0}$. We now suppose that ions from solution adsorb
to the polar crystal faces with surface charge density $\mp\sigma^{(n)}$. The
reason for the “$(n)$” superscript will become clear shortly. To investigate
these surfaces, we imagine taking a cut of the system far from the edges of
the crystal, as indicated by the dotted box in Fig. 1a, and ignore any edge
effects. This amounts to assuming that the planes of charge span the entire
plane orthogonal to the surface normal. This “slab geometry,” which is shown
in greater detail in Fig. 1b, is similar to the textbook parallel plate
capacitor model Zangwill (2013), a fact we draw upon heavily in what follows.
Despite relying on the parallel plate capacitor model extensively, we should
acknowledge that neglecting edge effects in the slab geometry is not an
innocent omission; in Sec. III we will see that doing so ultimately gives rise
to a boundary condition that is absent in the system of interest (Fig. 1a).

Figure 1: Schematic of a polar crystal in solution. (a) We imagine that a
crystal exposing polar surfaces with charge density $\pm\sigma_{0}$ is
immersed in an electrolyte environment. We suppose that anions from solution
adsorb to the positive crystal surface, giving rise to adsorbed surface charge
density $-\sigma^{(n)}$. Similarly, at the negative crystal surface the
adsorbed surface charge density is $+\sigma^{(n)}$. At some macroscopic
distance away from the crystal, there is a boundary between the electrolyte
and its vapor phase. To investigate the properties of these polar surfaces, we
imagine taking a cut of the system far from the edges of the crystal, as
indicated by the dashed box, and ignore edge effects; this is the “slab
geometry.” (b) In the slab geometry, the polar crystal is modeled as $n+1$
planes with alternating charge density $\pm\sigma_{0}$, separated by a
distance $R$. The specific case of $n=5$ is shown. We will also consider cases
where the central planes are separated by a distance $R_{\rm c}$. The planes
of adsorbed surface charge density (dotted lines) are separated from the
positive and negative termini by $\ell_{+}$ and $\ell_{-}$, respectively. In
regions marked with filled circles, the electric field is
$E=4\pi\sigma^{(n)}$, while in those marked with empty squares
$E=4\pi(\sigma^{(n)}-\sigma_{0})$.

Also given in Fig. 1b is a more detailed description of the structure of the
polar crystal. For simplicity, we consider halite (111)—relevant to, e.g.,
MgO, NiO, and NaCl—in which $n+1$ crystal planes are separated by a distance
$R$, where $n$ is an odd integer, though results can be generalized to other
geometries. Later, we will also perform “simulation experiments” in which we
change the separation between the central planes in the crystal to $R_{\rm
c}\neq R$. In the slab geometry, the underlying average microscopic charge
density $\rho_{\rm c}(z)$ only varies in the direction normal to the surface,
which we indicate by $z$. The total adsorbed surface charge density is then
determined by

$\sigma^{(n)}=\int_{{\rm int}-}\mathrm{d}z\,\rho_{\rm c}(z)=-\int_{{\rm
int}+}\mathrm{d}z\,\rho_{\rm c}(z),$ (1)

where “int$-$” and “int$+$” indicate that the integrations are performed over
interfacial regions corresponding to the negatively and positively charged
crystal surfaces, respectively.111Equation 1 implicitly assumes that an equal
number of crystal planes with $\sigma_{0}$ and $-\sigma_{0}$ are included in
the domain of integration. While Eq. 1 provides a direct means to obtain
$\sigma^{(n)}$ from average microscopic properties, deciding where to locate
the planes of adsorbed surface charge density is less straightforward.
Nonetheless, it seems reasonable to place the plane with charge density
$\sigma^{(n)}$ a distance $\ell_{-}$ from the crystal’s negatively charged
surface, and that with $-\sigma^{(n)}$ a distance $\ell_{+}$ from the
crystal’s positively charged surface, where $\ell_{-}$ and $\ell_{+}$ both
have a clear microscopic interpretation, inasmuch as they can be related to
features of $\rho_{\rm c}(z)$. In what follows, we will in fact find an
optimal _apparent_ length scale entering the continuum model that has no such
clear microscopic interpretation.

Let us denote the free energy per unit area of the system as
$f(\sigma^{(n)})$. Despite discussing the system in terms of uniform planes of
charge, we have already acknowledged that the underlying system comprises
atomic and molecular entities. It is natural, then, to partition $f$ into
“capacitive” and “noncapacitive” contributions:

$f(\sigma^{(n)})=u_{\rm cap}(\sigma^{(n)})+f_{\rm nc}(\sigma^{(n)}).$ (2)

The capacitive contribution, $u_{\rm cap}$, captures the energy stored in the
electric fields assuming that the system comprises uniform planes of charge,
as shown in Fig. 1b. The noncapacitive term, $f_{\rm nc}$, captures
contributions from everything else, e.g., non-electrostatic interactions,
electrostatic interactions not captured by the simple capacitor model, solvent
effects, and any entropic contributions.

To find $u_{\rm cap}$, we view the crystal in the slab geometry as $(n+1)/2$
parallel plate capacitors, between which the electric field is
$E=4\pi(\sigma^{(n)}-\sigma_{0})$, along with $(n-1)/2$ capacitors between
which the electric field is $E=4\pi\sigma^{(n)}$. For both sets of capacitors,
the separation between plates is $R$. (For the moment, we consider $R_{\rm
c}=R$.) In addition, the electric field between an adsorbed plane and the
surface of the crystal is also $E=4\pi\sigma^{(n)}$. Recalling that the energy
density of an electric field is $|E|^{2}/8\pi$, it follows that

$\displaystyle u_{\rm cap}$
$\displaystyle=2\pi(nR+2\ell)\big{(}\sigma^{(n)}\big{)}^{2}$
$\displaystyle-2\pi(n+1)R\sigma_{0}\sigma^{(n)}+\pi(n+1)R\sigma_{0}^{2},$ (3)

where $2\ell=l_{+}+l_{-}$. (Throughout this formulation, we work in a unit
system in which $4\pi\epsilon_{0}=1$, where $\epsilon_{0}$ is the permittivity
of free space.) At equilibrium, $\sigma^{(n)}$ assumes a value that minimizes
$f$:

$\sigma^{(n,\rm eq)}=\frac{(n+1)R\sigma_{0}-f^{\prime}_{\rm
nc}/2\pi}{2nR+4\ell},$ (4)

where the prime indicates a partial derivative with respect to $\sigma^{(n)}$.
Equation 4 is a central result of this article. While noncapacitive
contributions are irrelevant as $n\to\infty$, their effects become
increasingly important as the thickness of the crystal decreases.

Although formally exact, analysis of Eq. 4 is complicated by the presence of
the noncapacitive contributions. To make exploratory progress, we simply
postulate that

$f_{\rm nc}(\sigma^{(n)})=4\pi a_{\rm nc}\big{(}\sigma^{(n)}\big{)}^{2}.$ (5)

The equilibrium adsorbed surface charge density then reads

$\sigma^{(n,\rm eq)}=\frac{(n+1)R\sigma_{0}}{2nR+4(\ell+a_{\rm nc})}.$ (6)

Within the phenomenological model specified by Eq. 5, one can view the effects
of noncapacitive contributions as modifying the apparent length scale
associated with ion adsorption, from $\ell$ to $\ell+a_{\rm nc}$, in the
continuum representation of the system. A similar expression for
$\sigma^{(n,\rm eq)}$ was derived previously by Hu Hu (2021) for the slab
geometry under periodic boundary conditions, within the framework of symmetry-
preserving mean field theory Pan, Yi, and Hu (2019); Hu (2014). Within that
framework, an effective length scale is rationalized by considering the slowly
varying components of $\rho_{\rm c}$, whose structural features need not
coincide with those of the full average microscopic charge density. We will
postpone further discussion of the effective length scale $\ell+a_{\rm nc}$ to
Sec. II.2 and Appendix A. Importantly, the physical interpretation of Eq. 6 is
somewhat different from that of Ref. Hu, 2021, where
$\sigma^{(n)}=\sigma^{(n,\rm eq)}$ was attributed to the presence of a
spurious electric field arising from the use of tinfoil Ewald sums. In
contrast, in this study Eq. 6 has been derived for a nonperiodic slab
geometry; the driving force for ion adsorption is a significant reduction in
the capacitive energy stored in the system upon ion adsorption.

To understand the origin of the polar catastrophe from this free energy
perspective, we consider the reversible work required to change the adsorbed
surface charge density by an amount $\delta\sigma$ from its equilibrium value,

$f(\sigma^{(n,\rm eq)}+\delta\sigma)=f(\sigma^{(n,\rm
eq)})+\frac{1}{2}f^{\prime\prime}(\sigma^{n,\rm
eq})(\delta\sigma)^{2}+\mathcal{O}\big{(}(\delta\sigma)^{3}\big{)}.$

Neglecting terms higher than second order in $\delta\sigma$, and combining
with Eqs. 2, II.1 and 5 gives

$f(\sigma^{(n,\rm eq)}+\delta\sigma)-f(\sigma^{(n,\rm
eq)})=2\pi\big{(}nR+2(\ell+a_{\rm nc})\big{)}(\delta\sigma)^{2}.$ (7)

Equation 7 is another key result of this study; it is analogous to the
familiar diverging surface energy of an unreconstructed polar crystal in
contact with vacuum Tasker (1979); Noguera (2000). We see that the reversible
work required to change the adsorbed surface charge density from
$\sigma^{(n,\rm eq)}$ by any finite amount diverges linearly with the crystal
thickness $nR$.

Despite the clear analogy to the well-known diverging surface energy of polar
crystal surfaces in vacuum, Eq. 7 contains greater information. This can be
seen more clearly by considering the behavior of the fluctuations of the total
charge $Q$ within a probe area $A^{(n)}$ at the surface of the crystal. At
equilibrium, the total adsorbed charge within such an area is $Q^{(n,\rm
eq)}=A^{(n)}\sigma^{(n,\rm eq)}$. Eq. 7 states that the area required to
observe $\delta Q=Q-Q^{(n,\rm eq)}\approx e$ with appreciable probability,
i.e., $\beta A^{(n)}[f(\sigma^{(n,\rm eq)}+\delta\sigma)-f(\sigma^{(n,\rm
eq)})]\approx 1$, is,

$A^{(n)}\approx 2\pi\beta e^{2}\big{(}nR+2(\ell+a_{\rm nc})\big{)},$ (8)

where $e$ is the unit of elementary charge, and $\beta=1/k_{\rm B}T$ ($k_{\rm
B}$ is Boltzmann’s constant, and $T$ is the temperature). As the thickness $n$
of the crystal increases, ever increasing probe areas are required to observe
appreciable fluctuations in the adsorbed surface charge. The theory presented
here not only reveals the strong constraint placed on the average adsorbed
surface charge density at polar crystal surfaces, but also that fluctuations
at the interface are strongly suppressed.

###  II.2 Validation of the theory with molecular simulations

Equation 6 provides an expression for the equilibrium adsorbed surface charge
density at polar halite (111) surfaces in contact with an electrolyte
solution, assuming an ionic model for the crystal. As such, Eq. 6 lends itself
naturally to validation with molecular simulations. But two comments are in
order. The first concerns the treatment of long-range electrostatic
interactions under periodic boundary conditions, for which lattice summation
techniques are often employed. For Ewald summation, which underlies many of
the most popular techniques used in molecular simulations, the expression for
the Coulomb energy for a set of charges $\\{q_{i}\\}$ with positions
$\\{\mathbf{r}_{i}\\}$ reads:

$\displaystyle U_{\rm C}$
$\displaystyle=\frac{1}{2}\sum_{\rm\mathbf{b}}\sideset{}{{}^{\prime}}{\sum}_{i,j}q_{i}q_{j}\phi_{\rm
SR}(\mathbf{r}_{i}-\mathbf{r}_{j}-\mathbf{b})$
$\displaystyle+\frac{1}{2\Omega}\sum_{\mathbf{k}\neq\mathbf{0}}|\tilde{\rho}(\mathbf{k})|^{2}\tilde{\phi}_{\rm
LR}(\mathbf{k})-\frac{1}{2}\sum_{i}q_{i}^{2}\phi_{\rm LR}(\mathbf{0})$
$\displaystyle+\frac{1}{2\Omega}\mathbf{M}\cdot\mathbf{J}\cdot\mathbf{M}.$ (9)

In Eq. 9, $\mathbf{b}$ is a lattice vector, and the prime indicates that the
term $i=j$ is omitted when $\mathbf{b}=\mathbf{0}$. The volume of the
simulation cell is $\Omega$, and $\tilde{\rho}(\mathbf{k})$ is the Fourier
transform of the charge density
$\rho(\mathbf{r})=\sum_{i}q_{i}\delta(\mathbf{r}-\mathbf{r}_{i})$. The
potentials $\phi_{\rm SR}(\mathbf{r})$ and $\phi_{\rm LR}(\mathbf{r})$ are
defined by a splitting of the Coulomb potential into short- and long-ranged
contributions, $1/|\mathbf{r}|=\phi_{\rm SR}(\mathbf{r})+\phi_{\rm
LR}(\mathbf{r})$, and $\tilde{\phi}_{\rm LR}$ is the Fourier transform of the
long-ranged part. Most important for the current discussion is the final
“surface term,” which depends upon the total dipole moment of the simulation
cell, $\mathbf{M}$, and the depolarization tensor, $\mathbf{J}$. This surface
term is determined by the specified summation order of the lattice sum, and
for a spherical summation order it is equal to
$2\pi|\mathbf{M}|^{2}/(2\epsilon^{\prime}+1)$, where $\epsilon^{\prime}$ is
the dielectric constant of a surrounding medium “at infinity.” Tinfoil Ewald
approaches take this surrounding medium to be a perfect conductor,
$\epsilon^{\prime}=\infty$, and amounts to ignoring the surface term. The
above is not intended as a detailed discussion of Ewald summation, for which
there is an extensive literature de Leeuw, Perram, and Smith (1980a, b); Smith
(1981); De Leeuw, Perram, and Smith (1986); Neumann (1983); Neumann and
Steinhauser (1983a, b, 1984); Redlack and Grindlay (1975); Kantorovich and
Tupitsyn (1999); Smith (2008). The reader is referred to Ref. Ballenegger,
2014 for a particularly clear discussion on the topic.

For systems in the slab geometry under periodic boundary conditions employing
tinfoil Ewald, where the slab is surrounded on either side by vacuum, it has
long been recognized that if a net dipole $M_{z}$ exists along $z$, then there
is a finite electric field in the vacuum region. It is therefore common to
employ “slab corrections” to mitigate possible spurious interactions between
periodic images. In classical simulations of liquid-solid interfaces, the
standard choice is the method of Yeh and Berkowitz Yeh and Berkowitz (1999),
which includes a surface term $2\pi M_{z}^{2}/\Omega$. A similar expression
exists in the surface science community (which typically uses periodic density
functional theory calculations)Neugebauer and Scheffler (1992); Bengtsson
(1999), where it is known as the “dipole correction.” In Ref. Sayer and Cox,
2020, however, it was demonstrated that typical slab correction schemes impose
vanishing adsorbed surface charge density; more generally, it was found that
an electric displacement field directly determines the adsorbed surface
charge, irrespective of the slab’s thickness. Such a scenario is clearly
incompatible with $\sigma^{(n,\rm eq)}$ predicted by Eq. 6. In contrast,
simulations performed with tinfoil Ewald approaches were found to give
adsorbed surface charge densities that increase with crystal thickness. We
will therefore compare the predictions of the theory presented above to
simulations that use a tinfoil Ewald approach; a theoretical justification for
doing so will be given in Sec. III, where it will also become clear why
standard slab corrections fail.

Figure 2: Theoretical predictions and molecular simulations paint a consistent
picture for ion adsorption at polar surfaces. (a) Symbols show $\sigma^{(n,\rm
eq)}$ obtained from 10 ns simulations of a halite slab exposing its (111) and
(11-1) surfaces to aqueous electrolyte solution. The dot-dashed line indicates
a best fit of the theoretical prediction (Eq. 6) to the simulation data, with
$\ell=R$ and $a_{\rm nc}=-0.568$ Å, where $a_{\rm nc}$ is the only fitting
parameter. A model that ignores noncapacitive contributions ($a_{\rm nc}=0$)
gives a poor description of the simulation data. $A\approx 15.95\times 13.82$
Å2 is the cross-sectional area of the simulation cell. Error bars indicate 95
% confidence intervals, estimated by splitting trajectories into five samples.
(b) Average number density profiles $\rho_{\rm n}(z)$ for Na+ ions (solid
orange line) and water hydrogen atoms (dashed blue line) above halite (11-1).
(c) Average charge density $\rho_{\rm c}(z)$ above halite (11-1). In both (c)
and (d), the region occupied by the crystal is indicated by the shaded gray
region. The vertical purple dotted line indicates the separation $R$ from the
crystal’s surface, which aligns well with the first peaks in $\rho_{\rm n}(z)$
and $\rho_{\rm c}(z)$. The vertical blue dot-dashed line indicates the
separation $R+a_{\rm nc}$, which lacks a clear microscopic interpretation.

The second comment concerns the parameters that determine $\sigma^{(n,\rm
eq)}$. Here, we will consider crystals where the positions of the ions have
been clamped, and whose charges do not fluctuate; $\sigma_{0}$ and $R$ are
then straightforward input parameters. The extent to which agreement between
theory and simulation is observed is an indicator of how faithfully the theory
presented in this article describes the ionic model of polar surfaces.
Moreover, we will consider the specific example of NaCl (111) in contact with
an aqueous NaCl solution. We then anticipate that ions in direct contact with
the crystal will form a plane separated from the crystal by a distance
$\ell\approx R$. The only remaining unknown parameter is $a_{\rm nc}$, the
length scale associated with noncapacitive contributions.

In a first validation step, we compare $\sigma^{(n,\rm eq)}$ obtained from
molecular simulations of a concentrated aqueous NaCl solution in contact with
halite (111), where $a_{\rm nc}$ is simply treated as a free fitting
parameter. As seen in Fig. 2a, Eq. 6 with $\ell=R=1.628$ Å and $a_{\rm
nc}=-0.568$ Å captures the trend seen in the molecular simulations very well.
To determine the significance of noncapacitive contributions, it is
instructive to compare the effective length scale $\ell+a_{\rm nc}=1.06$ Å in
the continuum model to real length scales in the system. To that end, Fig. 2b
shows the number density profile of Na+ ions at halite (11-1) with $n=5$. As
expected, the plane separated from halite (11-1) by $R$ coincides with the
plane of Na+ directly adsorbed to the surface. In contrast, the plane at
$\ell+a_{\rm nc}$ from (11-1) is located in a region that cannot be readily
associated with ion adsorption; it does not appear that $a_{\rm nc}$ simply
captures effects due to thermal fluctuations.

The catch-all nature of $f_{\rm nc}$, and the _ad hoc_ assertion of its
quadratic form (Eq. 5), makes direct physical interpretation of the value of
$a_{\rm nc}$ challenging. Recent work suggests that dielectric boundaries for
systems with water as solvent are closely associated with water’s hydrogen
density Cox, Mandadapu, and Geissler (2021); Cox and Geissler (2022). Also
plotted in Fig. 2b, therefore, is the number density profile of water’s
hydrogen atoms. While some hydrogen density is found marginally closer to
halite (11-1) than that of the Na+ ions, it is difficult to draw any firm
conclusions on a possible relationship between $a_{\rm nc}$ and the solvent.
More likely is that noncapacitive contributions account for the rather drastic
approximation of collapsing the entire adsorbed surface charge density into a
single plane. As can be seen in Fig. 2c, $\rho_{\rm c}(z)$ exhibits pronounced
structure beyond the first peak at $\ell=R$. In Sec. III, we will see that the
effective length scale $\ell+a_{\rm nc}$ in the continuum model is that which
equalizes the electrostatic potential in the solution on either side of the
crystal, a condition that is automatically satisfied by molecular simulations
of a crystal slab immersed in an electrolyte solution, in which tinfoil Ewald
approaches are used (see Appendix A).

Despite the lack of a simple physical interpretation, the importance of
$f_{\rm nc}$ can be readily seen by plotting Eq. 6 with $a_{\rm nc}=0$, as
shown in Fig. 2a. Clearly, a model that only considers capacitive
contributions in which the system is approximated as planes of uniform charge
density poorly describes the simulation data, with predicted adsorbed surface
charge densities much lower than that observed.

A possible cause for concern is that simply treating $a_{\rm nc}$ as free
fitting parameter amounts to nothing more than a convenient fix. Challenges in
interpreting the value of $a_{\rm nc}$ notwithstanding, it is reasonable to
assume that noncapacitive contributions largely arise from effects that are
localized to the solution and interface regions. While changes to the
structure of the crystal in regions far from the interface will affect $u_{\rm
cap}$, we do not anticipate significant impact on the coefficient $a_{\rm
nc}$. In a second validation step, we therefore perform a “simulation
experiment” in which the separation between the crystal’s central planes is
varied (i.e., $R_{\rm c}\neq R$). The manner in which $\sigma^{(n,\rm eq)}$ is
affected by $R_{\rm c}\neq R$ depends on whether the central planes bound a
region with $E=4\pi\sigma^{(n)}$ [$(n+1)/2$ is even] or
$E=4\pi(\sigma^{(n)}-\sigma_{0})$ [$(n+1)/2$ is odd]. Following the same
approach leading to Eq. 6 gives

$\displaystyle\sigma^{(n,\rm eq)}=$ $\displaystyle\frac{\left[(n-1)R+2R_{\rm
c}\right]\sigma_{0}}{\left[(2n-2)R+2R_{\rm c}+4(\ell+a_{\rm nc})\right]}$
[$(n+1)/2$ is odd], (10a) $\displaystyle\sigma^{(n,\rm eq)}=$
$\displaystyle\frac{(n+1)R\sigma_{0}}{\left[(2n-2)R+2R_{\rm c}+4(\ell+a_{\rm
nc})\right]}$ [$(n+1)/2$ is even]. (10b)

Results obtained from simulations with $R_{\rm c}=3R/4$, $R/2$ and $R/4$ are
presented in Fig. 3, along with the theoretical predictions described by Eqs.
10a and 10b, in which no further fitting has been performed. That is, the same
value $a_{\rm nc}=-0.568$ Å has been used. The good agreement between the
theoretical predictions and the simulation data lends support to the notion
that fitting $a_{\rm nc}$ captures genuine noncapacitive effects that are
localized to the solution and interface regions.

Figure 3: “Simulation experiments” further validate the theory. Symbols show
$\sigma^{(n,\rm eq)}$ obtained from simulations with $R_{\rm c}\neq R$, as
indicated in the legend. Dashed and dotted lines indicate theoretical
predictions described by Eqs. 10a and 10b, respectively, with $a_{\rm
nc}=-0.568$ Å (no further fitting). The good agreement between simulation and
theory suggests that $a_{\rm nc}$ is insensitive to structural changes away
from the interface. The blue dot-dashed line is the same as Fig. 2a, i.e.,
$R_{\rm c}=R$.

Results presented so far demonstrate consistency between the theory presented
in Sec. II.1 for the equilibrium adsorbed surface charge density at halite
(111), and molecular simulations that use a tinfoil Ewald approach. In Sec.
III, we will discuss why it is appropriate to compare $\sigma^{(n,\rm eq)}$
for the system of interest (Fig. 1a) to simulations that employ the slab
geometry under periodic boundary conditions. In doing so, we will also shed
light on the role of the electric displacement field in the slab geometry
Zhang and Sprik (2016a); Zhang (2018); Sayer, Sprik, and Zhang (2019); Sayer
and Cox (2019); Zhang  _et al._ (2020); Sayer and Cox (2020).

##  III The emergence of the electric displacement field in a slab geometry

Recall the procedure employed in Sec. II.1 to investigate the polar surfaces
of a crystal in solution. First, a particle of finite size, with polar
surfaces exposed, was immersed in an electrolyte, and ions from solution were
imagined to adsorb to its surfaces. Then, we constructed the slab geometry by
taking a cut of the system far from the edges of the crystal, and considered
the free energy per unit area. Why did we adopt this procedure, rather than
starting immediately with the slab geometry? Consider again the system of
interest, which, to avoid notational clutter, is shown again in Fig. 4a. As
the particle has finite size, the ions adsorbed to the positive and negative
crystal faces originate from the same pool of ions in the surrounding
electrolyte solution. The remaining electrolyte solution is electroneutral. As
the electrolyte solution is conducting, the total electric field in its
interior vanishes, $\mathbf{E}_{\rm e}=\mathbf{0}$. Now consider the vapor
region. We will assume that the vapor pressure is sufficiently low that this
region can be approximated as vacuum. In the absence of an external field, the
total electric field there too vanishes, $\mathbf{E}_{\rm v}=\mathbf{0}$. It
then follows that the surface charge density at the electrolyte-vapor
interface is zero, i.e., $\sigma_{\rm
ev}=\hat{\mathbf{r}}\cdot{(\mathbf{E}_{\rm v}-\mathbf{E}_{\rm e})}=0$, where
$\hat{\mathbf{r}}$ is the outward normal to the electrolyte surface. This
scenario is possible with the procedure employed in Sec. II.1, as the
electrolyte solution remains uncharged. Moreover, there is a clear physical
pathway by which the ions have adsorbed from solution to the surface.

Imagine that we instead start directly with the slab geometry. As the crystal
spans the entire plane orthogonal to the surface normal, there are now two
distinct regions of electrolyte. Upon ion adsorption, electroneutrality must
be preserved in each region separately (the ions cannot pass through the
crystal). As the electrolyte is a conductor, the remaining charge in solution
must be located at the electrolyte-vapor boundary, $|\sigma_{\rm
ev}|=\sigma^{(n)}$, as indicated in Fig. 4b. As before, $\mathbf{E}_{\rm
e}=\mathbf{0}$, but the electrolyte is now polarized, with uniform
polarization $\hat{\mathbf{z}}\cdot\mathbf{P}_{\rm e}=\sigma^{(n)}$. The
electric displacement field in the electrolyte then satisfies

$\hat{\mathbf{z}}\cdot\mathbf{D}_{\rm e}=\hat{\mathbf{z}}\cdot(\mathbf{E}_{\rm
e}+4\pi\mathbf{P}_{\rm e})=4\pi\sigma^{(n)}.$ (11)

In the vapor, we assume that atomic density is sufficiently low that
$\mathbf{P}_{\rm v}=\mathbf{0}$. Matching the electric displacement field
across the boundaries leads us to conclude that

$\hat{\mathbf{z}}\cdot\mathbf{E}_{\rm v}=\hat{\mathbf{z}}\cdot(\mathbf{D}_{\rm
v}-4\pi\mathbf{P}_{\rm v})=4\pi\sigma^{(n)}.$ (12)

This result for $\mathbf{E}_{\rm v}$ in the slab geometry is most severe.
Integrating the electrostatic energy density $|\mathbf{E}_{\rm v}|^{2}/8\pi$
over the macroscopic volume occupied by the vapor phase indicates a
prohibitively expensive energy cost for $\sigma^{(n)}\neq 0$. Another route to
arrive at the conclusion $\sigma^{(n)}=0$ _when starting in the slab geometry_
is that, as the electric field external to the parallel plate capacitors
vanishes, there is no driving force for ion adsorption.

Figure 4: Emergence of the electric displacement field in a slab geometry. (a)
The adsorbed surface charge densities $\sigma^{(n)}$ and $-\sigma^{(n)}$
comprise ions originating from the same pool of electrolyte, which remains
electroneutral. Both the surface charge density at the electrolyte-vapor
interface, $\sigma_{\rm ev}$, and the electric field in vacuum,
$\mathbf{E}_{\rm v}$, are zero. As the electrolyte is conducting, its interior
is equipotential, $V_{\rm a}=V_{\rm b}$. In contrast, the slab geometry
creates two distinct regions of electrolyte either side of the crystal slab,
as shown in (b). Unlike Fig. 1b, we have now included the vapor phase, and
only the outermost crystal planes are shown explicitly. As electroneutrality
must be preserved in each region separately, $\sigma^{(n)}=\sigma^{(n,\rm
eq)}$ implies a polarized electrolyte with a finite displacement field,
$\hat{\mathbf{z}}\cdot\mathbf{D}_{\rm e}=4\pi\sigma^{(n,\rm eq)}$. While this
satisfies $V_{\rm a}=V_{\rm b}$, the electric field in vapor is now finite,
$\hat{\mathbf{z}}\cdot\mathbf{E}_{\rm v}=4\pi\sigma^{(n,\rm eq)}$. If
$\sigma^{(n)}=\hat{\mathbf{z}}\cdot\mathbf{E}_{\rm v}=0$, then $V_{\rm
a}-V_{\rm b}=2\pi(n+1)R\sigma_{0}$. Panel (c) depicts the system under a
_periodic_ slab geometry, without the vapor region. While for
$\sigma^{(n)}=\sigma^{(n,\rm eq)}$ it is still the case that
$\hat{\mathbf{z}}\cdot\mathbf{D}_{\rm e}=4\pi\sigma^{(n,\rm eq)}$, the
electric field in the electrolyte vanishes, $\mathbf{E}_{\rm e}=\mathbf{0}$,
and $\mathbf{D}_{\rm e}\cdot\mathbf{E}_{\rm e}=0$. Moreover, $V_{\rm a}=V_{\rm
b}$.

The slab geometry, despite widespread use in both theoretical and computer
simulation studies, does not correspond to a physical reality. Consider the
electrostatic potential in the regions of the electrolyte indicated by the
dotted and dot-dashed circles in Fig. 4a, which we denote $V_{\rm a}$ and
$V_{\rm b}$, respectively. As the electrolyte is conducting, it follows
immediately that $V_{\rm a}=V_{\rm b}$. In the slab geometry, the crystal
provides an insulating layer between the two distinct regions of electrolyte.
For $\sigma^{(n)}=0$, the electrostatic potential difference across the
crystal is,

$V_{\rm a}-V_{\rm b}=2\pi(n+1)R\sigma_{0}.$ (13)

Even for $n=1$, this electrostatic potential difference is substantial:
$V_{\rm a}-V_{\rm b}\approx 20$ V. In contrast, for
$\sigma^{(n)}=\sigma^{(n,\rm eq)}$ (Eq. 6),

$V_{\rm a}-V_{\rm b}=0,$ (14)

where, as discussed in Appendix A, we have treated the adsorbed planes as if
they are separated from the crystal’s surfaces by the effective distance
$\ell+a_{\rm nc}$. The slab geometry fundamentally breaks the system: for
$\sigma^{(n)}=0$, $\mathbf{E}_{\rm v}=\mathbf{0}$ is accurately captured, but
there exists an erroneous potential difference across the crystal; for
$\sigma^{(n)}=\sigma^{(n,\rm eq)}$, the potential difference across the
crystal is correct, but at the cost of introducing a finite field in the vapor
phase, $\hat{\mathbf{z}}\cdot\mathbf{E}_{\rm v}=4\pi\sigma^{(n,\rm eq)}$. This
is a consequence of the finite electric displacement field in the electrolyte,
$\hat{\mathbf{z}}\cdot\mathbf{D}_{\rm e}=4\pi\sigma^{(n,\rm eq)}$ (Eq. 11),
that emerges in the slab geometry.

What are the implications for molecular simulations? Let us consider cases
where periodic boundary conditions are used in all three Cartesian directions,
restricting ourselves to cases where periodic replicas of the crystal slab
sandwich the electrolyte solution such that there is no vapor region, as shown
in Fig. 4c. Similar to the nonperiodic slab geometry, Fig. 4b, the electrolyte
is polarized, $\hat{\mathbf{z}}\cdot\mathbf{P}_{\rm e}=\sigma^{(n)}$, and
there is a finite electric displacement field in the electrolyte,
$\hat{\mathbf{z}}\cdot\mathbf{D}_{\rm e}=4\pi\sigma^{(n)}$. For simulations
that employ tinfoil Ewald approaches to compute electrostatic interactions,
there is no electrostatic potential difference across the simulation cell. As
$\mathbf{E}_{\rm e}=\mathbf{0}$, it then follows that $V_{\rm a}-V_{\rm b}=0$.
(In Appendix A, it is shown explicitly that $\sigma^{(n,\rm eq)}$ is the
adsorbed surface charge density that enforces $\mathbf{E}_{\rm e}=\mathbf{0}$
when using tinfoil Ewald approaches.) By not including a vapor region between
periodic images, the electrostatic energy is well behaved; as $\mathbf{D}_{\rm
e}\cdot\mathbf{E}_{\rm e}=0$, we can extend the electrolyte region between
periodic replicas without incurring any cost in electrostatic energy.
Simulations using tinfoil Ewald sums thus appear to give a faithful
representation of polar surfaces in solution—in the sense of taking a cut of
the system of interest as shown in Fig. 1a, albeit with a finite polarization
and electric displacement field in the electrolyte, the effects of which,
however, seem largely benign.

Implicit in the preceding discussion is that the polarization is treated in an
itinerant fashion—the ions are included in its definition—which has been shown
Caillol (1994); Cox and Sprik (2019) to satisfy key statistical mechanical
properties for electrolyte solutions, such as the Stillinger-Lovett sum rules
Stillinger Jr and Lovett (1968a, b). Any discussion on the multivaluedness of
polarization under periodic boundary conditions has also been neglected, along
with implications of whether the crystal or the electrolyte straddles the
boundaries of the simulation cell. The reader is referred to Refs. Zhang and
Sprik, 2016a; Sprik, 2018; Sayer, Sprik, and Zhang, 2019 for a discussion of
these issues in the context of molecular dynamics simulations.

As already mentioned, Ref. Sayer and Cox, 2020 demonstrated that imposing an
electric displacement field directly determines the adsorbed surface charge
density. Theoretical discussion in that work, however, was brief, and what
little there was relied heavily upon the framework provided by the “finite
field approach” of Zhang and Sprik Zhang and Sprik (2016b, a); Zhang  _et al._
(2020). The picture to emerge from the finite field approach is that of a set
of “virtual electrodes” that either controls the electrostatic potential
difference across the simulation cell (“constant $E$”) or controls the
electric displacement field at the cell boundaries (“constant $D$”). As the
finite field approach was formulated in the context of periodic boundary
conditions, it can be challenging to disentangle its deeper significance from
issues concerning conditional convergence of lattice sums and unwanted
interactions between periodic images. Moreover, the picture of virtual
electrodes can give the impression that one is forcing the system adopt a
particular adsorbed surface charge density with an arbitrary field. In this
article, the status of the electric displacement field as a control variable
has been clarified: it should be viewed as matching an electrostatic boundary
condition that emerges in the slab geometry—even in the absence of periodic
boundary conditions—rather than an _ad hoc_ external electric field that is
applied to the system. In Ref. Sayer and Cox, 2020, the surfaces of
macroscopically thick crystals in contact with an electrolyte solution were
effectively simulated by imposing $\hat{\mathbf{z}}\cdot\mathbf{D}_{\rm
e}=4\pi\sigma^{(\infty,\rm eq)}$, using the finite field approach. Similar
control of the adsorbed surface charge with an electric displacement field has
also been observed to work away from equilibrium Dufils, Sprik, and Salanne
(2021).

To end this discussion, it is worth considering standard slab correction
schemes again Yeh and Berkowitz (1999); Neugebauer and Scheffler (1992);
Bengtsson (1999), which typically introduce a vacuum region between periodic
replicas along $z$, and remove the electric field in the vacuum region by
enforcing $\hat{\mathbf{z}}\cdot\mathbf{D}_{\rm v}=0$. While such a line of
attack seems reasonable from the perspective of removing interactions between
periodic images, it presupposes that the true physical scenario corresponds to
an isolated slab in contact with the liquid (e.g., Fig. 4b). Slabs that
completely span the plane orthogonal to the surface normal, however, do not
exist; this fact cannot be ignored when modeling the surfaces of polar
crystals.

##  IV Conclusions

A theoretical framework to describe polarity compensation arising from the
solution environment has been presented in which the free energy of the system
has been separated into “capacitive” and “noncapacitive” contributions. For
crystals that are thin along a polar crystallographic direction, noncapacitive
contributions were found to play a significant role. As the thickness of the
crystal increases, the capacitive contribution dominates, and the adsorbed
surface charge density at equilibrium agrees with known expressions for the
required charge imbalance at polar surfaces Tasker (1979); Noguera (2000).

The theory presented here only attempts to capture the net adsorbed charge at
each polar crystal surface, and says nothing about the complex structural
details of the solution near the interface. A major source of the noncapactive
contributions likely arises from the drastic approximation of assuming that
all adsorbed charge lies in a single plane. In addition, the theory assumes an
ionic model for the crystal. While the severity of such an approximation will
depend on the system under investigation Noguera (2000), it may nonetheless
provide a useful starting point for more sophisticated theoretical approaches.

Like previous theoretical treatments for polar crystal surfaces in contact
with vacuum Tasker (1979); Noguera (2000), the presented theory relies heavily
on modeling the crystal as a stack of parallel plate capacitors. It was shown
that in such a slab geometry, an electric displacement field arises naturally
in the solution. This electric displacement field is a consequence of the
artificial change of topology imposed on the system, in which the crystal
blocks the passage of ions, resulting in a uniformly polarized electrolyte.
Our analysis reveals that this is a consequence of the slab geometry itself,
rather than the lattice summation techniques that are typically used to treat
electrostatic interactions under periodic boundary conditions. We argue that
tinfoil Ewald approaches, which impose zero electropotential difference in the
electrolyte either side of the crystal, are appropriate when attempting to
model the polar surfaces of crystals in solution, if the crystal is genuinely
thin. In contrast, if one is interested in crystals that are macroscopic in
extent, the electric displacement field can be chosen to enforce
$\sigma^{(\infty,\rm eq)}$, as demonstrated in Ref. Sayer and Cox, 2020.

Finally, our theoretical prediction (Eq. 7) for the fluctuations in the
adsorbed surface charge density is intriguing, and suggests that dynamics at
the interface may change significantly as the thickness of the crystal
increases. Static energy calculations already reveal that the behavior of thin
films in contact with vacuum can be complex, with structural relaxations that
we have not considered in this work permitting uncompensated polarity to a
certain extent Goniakowski, Noguera, and Giordano (2007). The morphology of
crystals is generally determined by the relative growth rates along different
crystallographic directions Dandekar, Kuvadia, and Doherty (2013). Considering
the potential interplay between structural relaxations, and the thickness
dependence of both the equilibrium adsorbed surface charge density and its
fluctuations, one can easily imagine that the growth mechanisms of a polar
surface from either a melt or supersaturated solution will be a complex
affair. In addition, we have implicitly assumed $\sqrt{\mathcal{A}}\gg nR$
throughout, where $\mathcal{A}$ is the surface area of the exposed polar
facets in the crystal (Fig. 1a), whereas the polar catastrophe has nontrivial
dependence on crystal morphology Nosker, Mark, and Levine (1970); Noguera and
Goniakowski (2013). Future work will focus on disentangling these different
contributions to the growth mechanisms of polar crystal surfaces.

##  V Methods

The simulation methodology closely resembles that of Ref. Sayer and Cox, 2020,
and indeed, for the $R_{\rm c}=R$ system, we reused trajectories for $n=3$,
$5$, $7$, $11$ and $23$. All simulations used the SPC/E water model Berendsen,
Grigera, and Straatsma (1987), whose geometry was constrained using the RATTLE
algorithm Andersen (1983) and the Joung-Cheatham NaCl force field Joung and
Cheatham III (2008). Dynamics were propagated using the velocity Verlet
algorithm with a time step of 2 fs. The temperature was maintained at 298 K
with a Nosé-Hoover chain Shinoda, Shiga, and Mikami (2004); Tuckerman  _et
al._ (2006), with a 0.2 ps damping constant. The particle-particle particle-
mesh Ewald method was used to account for long-ranged interactions Hockney and
Eastwood (1988), with parameters chosen such that the root mean square error
in the forces were a factor $10^{5}$ smaller than the force between two unit
charges separated by a distance of 1 ÅKolafa and Perram (1992). A cutoff of 10
Å was used for non-electrostatic interactions.

The electrolyte comprised 600 water molecules and 20 NaCl ion pairs. The
crystal consisted of alternating layers of Na+ and Cl- ions, separated by
$R=1.628$ Å, and each layer comprised 16 ions. The lateral dimensions of the
simulation cell were $L_{x}=15.952$ Å and $L_{y}=13.815$ Å along $x$ and $y$,
respectively. In the slab geometry with $n=3$ and $R_{\rm c}=R$, the length of
the simulation cell along $z$ was $L=94.841$ Å, and $L$ was increased with $n$
accordingly e.g. for $n=5$, $L$ was increased by $2R$. Each simulation was 10
ns long post equilibration. The LAMMPS simulation package was used throughout
Plimpton (1995).

###### Acknowledgements.

I am grateful to Rob Jack, Michiel Sprik and Tom Sayer for many enlightening
discussions on this topic, and Rick Remsing for reading a draft of the
manuscript. This work was supported by a Royal Society University Research
Fellowship (URF\R1\211144). For the purpose of open access, I have applied a
Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript
version arising.

## Author Declarations

### Conflict of Interest

The author has no conflicts to disclose.

## Data Availability

The data that support the findings of this study are openly available in the
data repository at https://doi.org/10.17863/CAM.83994.

##  Appendix A Derivation of $\sigma^{(n,\rm eq)}$ under tinfoil Ewald
periodic boundary conditions

In the main article, $\sigma^{(n,\rm eq)}$ was derived for a nonperiodic slab
geometry (Eqs. 4 and 6). In this appendix, an expression for $\sigma^{(n,\rm
eq)}$ is derived by considering the average molecular charge distribution
under periodic boundary conditions. The derivation strongly hints that
noncapacitive contributions arise from approximating $\sigma^{(n,\rm eq)}$ as
being confined to a single plane close to the crystal’s surface (Fig. 1). The
derivation largely resembles that of Ref. Hu, 2021, but differs in that we do
not consider a smeared average molecular charge distribution. Moreover, we
take the present derivation as evidence that tinfoil Ewald sums provide a
reliable estimate of $\sigma^{(n,\rm eq)}$, rather than ion adsorption
resulting from a spurious treatment of electrostatic interactions.

When the average molecular charge density $\rho_{\rm c}$ only varies along
$z$, the average electrostatic potential under tinfoil Ewald periodic boundary
conditions is Wirnsberger  _et al._ (2016); Pan, Yi, and Hu (2017); Cox and
Geissler (2018)

$V(z)=4\pi\int_{\rm cell}\\!\mathrm{d}z^{\prime}\,\rho_{\rm
c}(z^{\prime})J(z-z^{\prime}),$ (15)

with

$J(z)=\mathrm{const.}+\frac{z^{2}}{2L}-\frac{|z|}{2},$ (16)

where $L$ is the length of the simulation cell along $z$. For simplicity, let
us first consider the case $n=1$. In the simulations discussed in this
article, the crystal ions do not move; approximating these as uniformly
charged planes we have,

$\rho_{\rm c}(z)=\rho_{\rm c}^{(\rm
soln)}(z)+\sigma_{0}\big{[}\delta(z-R/2)-\delta(z+R/2)\big{]},$ (17)

where $\rho_{\rm c}^{\rm(soln)}$ is the average molecular charge distribution
arising from the electrolyte solution (both solvent and ions). Let us assume
that the crystal is centered at $z=0$. In the region $R/2<z\leq L/2$ we then
have

$V(z)=4\pi\int_{\rm cell}\\!\mathrm{d}z^{\prime}\,\rho^{\rm(soln)}_{\rm
c}(z^{\prime})J(z-z^{\prime})+4\pi\sigma_{0}\bigg{[}\frac{-Rz}{L}+\frac{R}{2}\bigg{]},$
(18)

and an average electric field

$\hat{\mathbf{z}}\cdot\mathbf{E}(z)=-4\pi\int_{\rm
cell}\\!\mathrm{d}z^{\prime}\,\rho^{\rm(soln)}_{\rm
c}(z^{\prime})\frac{\mathrm{d}J(z-z^{\prime})}{\mathrm{d}z}+4\pi\sigma_{0}\frac{R}{L}.$
(19)

At equilibrium, the electric field in the electrolyte vanishes, thus

$\int_{\rm cell}\\!\mathrm{d}z^{\prime}\,\rho^{\rm(soln,eq)}_{\rm
c}(z^{\prime})\frac{\mathrm{d}J(z-z^{\prime})}{\mathrm{d}z}=\sigma_{0}\frac{R}{L},$
(20)

where $z$ is understood to represent a point in the bulk electrolyte. The
equilibrium electrolyte charge distribution satisfies Eq. 20, no matter how
complex its behavior. For a general $\rho_{\rm c}^{(\rm soln,eq)}$, however,
Eq. 20 is not straightforward to analyze.

To simplify matters, we introduce a model for $\rho_{\rm c}^{(\rm soln)}$
whereby the adsorbed surface charge at each surface is confined to a single
plane:

$\rho_{\rm c}^{(\rm
soln)}(z)\approx\sigma^{(n=1)}\big{[}\delta\big{(}z+(R+2\bar{\ell}\big{)}/2)-\delta\big{(}z-(R+2\bar{\ell})/2\big{)}\big{]}.$
(21)

The relationship between $\bar{\ell}$ and the true average molecular charge
distribution is left unspecified, and it should not necessarily be considered
a physically meaningful length scale. For this model, the average
electrostatic potential for $(R+2\bar{\ell})/2<z\leq L/2$ is

$\displaystyle V(z)=$
$\displaystyle-4\pi\sigma^{(n=1)}\bigg{[}\frac{-(R+2\bar{\ell})z}{L}+\frac{R+2\bar{\ell}}{2}\bigg{]}$
$\displaystyle+4\pi\sigma_{0}\bigg{[}\frac{-Rz}{L}+\frac{R}{2}\bigg{]},$ (22)

while the electric field is

$\hat{\mathbf{z}}\cdot\mathbf{E}_{\rm
e}=-4\pi\sigma^{(n=1)}\frac{R+2\bar{\ell}}{L}+4\pi\sigma_{0}\frac{R}{L}.$ (23)

Again enforcing the equilibrium condition,
$\hat{\mathbf{z}}\cdot\mathbf{E}_{\rm e}=0$, we find

$\sigma^{(n=1,\rm eq)}=\frac{\sigma_{0}}{1+2\bar{\ell}/R},$ (24)

which agrees with Eq. 6, provided that $\bar{\ell}=\ell+a_{\rm nc}$. Comparing
Eqs. 6, 20 and 24 suggests that noncapacitive contributions largely account
for the approximation of treating the entire adsorbed surface charge density
as if it were confined to a single plane.

Following the same procedure for general $n$, we find

$\displaystyle\hat{\mathbf{z}}\cdot\mathbf{E}_{\rm
e}=-4\pi\sigma^{(n)}\frac{(nR+2\bar{\ell})}{L}$
$\displaystyle-4\pi\sigma_{0}\frac{R}{L}(-1)^{(n+1)/2}\sum_{j=1}^{(n+1)/2}(2j-1)(-1)^{j+1}$
$\displaystyle=-4\pi\sigma^{(n)}\frac{(nR+2\bar{\ell})}{L}+2\pi\sigma_{0}\frac{(n+1)R}{L}(-1)^{(n+1)/2}i^{n+1}$
$\displaystyle=-4\pi\sigma^{(n)}\frac{(nR+2\bar{\ell})}{L}+2\pi\sigma_{0}\frac{(n+1)R}{L}.$
(25)

To arrive at this result, we recall that $n+1$ is an even integer, such that
$i^{n+1}=(-1)^{(n+1)/2}$. At equilibrium,

$\sigma^{(n,\rm eq)}=\frac{(n+1)\sigma_{0}}{2nR+4\bar{\ell}}.$ (26)

Again, we find agreement with Eq. 6, provided that $\bar{\ell}=\ell+a_{\rm
nc}$. As we derived $\sigma^{(n,\rm eq)}$ by enforcing
$\hat{\mathbf{z}}\cdot\mathbf{E}_{\rm e}=0$, and as $V(-L/2)=V(L/2)$ (Eqs. 15
and 16), it immediately follows that $V_{\rm a}=V_{\rm b}$ in the bulk
electrolyte either side of the slab (Fig. 4, see also Fig. 4 of Ref Sayer and
Cox, 2020). Note that $\sigma^{(n,\rm eq)}$ is independent of $L$. In the
limit $L\to\infty$, any value of $\sigma^{(n)}$ satisfies
$\hat{\mathbf{z}}\cdot\mathbf{E}_{\rm e}=0$, which reflects the fact that the
electric field external to a set of infinite parallel plate capacitors
vanishes. For $\sigma^{(n)}\neq\sigma^{(n,\rm eq)}$, however, $V_{\rm a}\neq
V_{\rm b}$, and the system becomes increasingly unstable as $n$ increases (Eq.
7).

## References

  * Dandekar, Kuvadia, and Doherty (2013) P. Dandekar, Z. B. Kuvadia, and M. F. Doherty, Annu. Rev. Mater. Res. 43, 359 (2013).
  * Tasker (1979) P. Tasker, J. Phys. C: Solid State Phys. 12, 4977 (1979).
  * Noguera (2000) C. Noguera, J. Phys.: Condens. Matter 12, R367 (2000).
  * Goniakowski, Finocchi, and Noguera (2008) J. Goniakowski, F. Finocchi, and C. Noguera, Rep. Prog. Phys. 71, 016501 (2008).
  * Zhu  _et al._ (2006) K. Zhu, J. Hu, C. Kuebel, and R. Richards, Angew. Chem. Intl. Ed. 45, 7277 (2006).
  * Cheng, Akhter, and Kung (1983) W. Cheng, S. Akhter, and H. Kung, J. Catal. 82, 341 (1983).
  * Nosker, Mark, and Levine (1970) R. Nosker, P. Mark, and J. Levine, Surf. Sci. 19, 291 (1970).
  * Xu  _et al._ (2009) T. Xu, X. Zhou, Z. Jiang, Q. Kuang, Z. Xie, and L. Zheng, Cryst. Growth Des. 9, 192 (2009).
  * Susman  _et al._ (2018) M. D. Susman, H. N. Pham, A. K. Datye, S. Chinta, and J. D. Rimer, Chem. Mater. 30, 2641 (2018).
  * Susman  _et al._ (2020) M. D. Susman, H. N. Pham, X. Zhao, D. H. West, S. Chinta, P. Bollini, A. K. Datye, and J. D. Rimer, Angew. Chem. Intl. Ed. 59, 15119 (2020).
  * Sayer, Zhang, and Sprik (2017) T. Sayer, C. Zhang, and M. Sprik, J. Chem. Phys. 147, 104702 (2017).
  * Sayer, Sprik, and Zhang (2019) T. Sayer, M. Sprik, and C. Zhang, J. Chem. Phys. 150, 041716 (2019).
  * Sayer and Cox (2020) T. Sayer and S. J. Cox, J. Chem. Phys. 153, 164709 (2020).
  * Zangwill (2013) A. Zangwill,  _Modern Electrodynamics_ (Cambridge University Press, Cambridge, United Kingdom, 2013).
  * Note (1) Equation 1 implicitly assumes that an equal number of crystal planes with $\sigma_{0}$ and $-\sigma_{0}$ are included in the domain of integration.
  * Hu (2021) Z. Hu, “Comment on ‘Macroscopic surface charges from microscopic simulations’ [J. Chem. Phys. 153, 164709 (2020)],”  (2021), arXiv:2106.12291 [physics.chem-ph].
  * Pan, Yi, and Hu (2019) C. Pan, S. Yi, and Z. Hu, Phys. Chem. Chem. Phys. 21, 14858 (2019).
  * Hu (2014) Z. Hu, Chem. Commun. 50, 14397 (2014).
  * de Leeuw, Perram, and Smith (1980a) S. W. de Leeuw, J. W. Perram, and E. R. Smith, Proc. R. Soc. London A 373, 27 (1980a).
  * de Leeuw, Perram, and Smith (1980b) S. W. de Leeuw, J. W. Perram, and E. R. Smith, Proc. R. Soc. London A 373, 57 (1980b).
  * Smith (1981) E. R. Smith, Proc. R. Soc. London A 375, 475 (1981).
  * De Leeuw, Perram, and Smith (1986) S. De Leeuw, J. W. Perram, and E. Smith, Annu. Rev. Phys. Chem. 37, 245 (1986).
  * Neumann (1983) M. Neumann, Mol. Phys. 50, 841 (1983).
  * Neumann and Steinhauser (1983a) M. Neumann and O. Steinhauser, Chem. Phys. Lett. 95, 417 (1983a).
  * Neumann and Steinhauser (1983b) M. Neumann and O. Steinhauser, Chem. Phys. Lett. 102, 508 (1983b).
  * Neumann and Steinhauser (1984) M. Neumann and O. Steinhauser, Chem. Phys. Lett. 106, 563 (1984).
  * Redlack and Grindlay (1975) A. Redlack and J. Grindlay, J. Phys. Chem. Solids 36, 73 (1975).
  * Kantorovich and Tupitsyn (1999) L. Kantorovich and I. Tupitsyn, J. Phys.: Condens. Matter 11, 6159 (1999).
  * Smith (2008) E. R. Smith, J. Chem. Phys. 128, 174104 (2008).
  * Ballenegger (2014) V. Ballenegger, J. Chem. Phys. 140, 161102 (2014).
  * Yeh and Berkowitz (1999) I.-C. Yeh and M. L. Berkowitz, J. Chem. Phys. 111, 3155 (1999).
  * Neugebauer and Scheffler (1992) J. Neugebauer and M. Scheffler, Phys. Rev. B 46, 16067 (1992).
  * Bengtsson (1999) L. Bengtsson, Phys. Rev. B 59, 12301 (1999).
  * Cox, Mandadapu, and Geissler (2021) S. J. Cox, K. K. Mandadapu, and P. L. Geissler, J. Chem. Phys. 154, 244502 (2021).
  * Cox and Geissler (2022) S. J. Cox and P. L. Geissler, “Dielectric response of thin water films: A thermodynamic perspective,” (2022), DOI:10.1039/d2sc01243j.
  * Zhang and Sprik (2016a) C. Zhang and M. Sprik, Phys. Rev. B 94, 245309 (2016a).
  * Zhang (2018) C. Zhang, J. Chem. Phys. 149, 031103 (2018).
  * Sayer and Cox (2019) T. Sayer and S. J. Cox, Phys. Chem. Chem. Phys. 21, 14546 (2019).
  * Zhang  _et al._ (2020) C. Zhang, T. Sayer, J. Hutter, and M. Sprik, J. Phys.: Energy 2, 032005 (2020).
  * Caillol (1994) J.-M. Caillol, J. Chem. Phys 101, 6080 (1994).
  * Cox and Sprik (2019) S. J. Cox and M. Sprik, J. Chem. Phys. 151, 064506 (2019).
  * Stillinger Jr and Lovett (1968a) F. H. Stillinger Jr and R. Lovett, J. Chem. Phys. 48, 3858 (1968a).
  * Stillinger Jr and Lovett (1968b) F. H. Stillinger Jr and R. Lovett, J. Chem. Phys. 49, 1991 (1968b).
  * Sprik (2018) M. Sprik, Mol. Phys. 116, 3114 (2018).
  * Zhang and Sprik (2016b) C. Zhang and M. Sprik, Phys. Rev. B 93, 144201 (2016b).
  * Dufils, Sprik, and Salanne (2021) T. Dufils, M. Sprik, and M. Salanne, J. Phys. Chem. Lett. 12, 4357 (2021).
  * Goniakowski, Noguera, and Giordano (2007) J. Goniakowski, C. Noguera, and L. Giordano, Phys. Rev. Lett. 98, 205701 (2007).
  * Noguera and Goniakowski (2013) C. Noguera and J. Goniakowski, Chem. Rev. 113, 4073 (2013).
  * Berendsen, Grigera, and Straatsma (1987) H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987).
  * Andersen (1983) H. C. Andersen, J. Comput. Phys. 52, 24 (1983).
  * Joung and Cheatham III (2008) I. S. Joung and T. E. Cheatham III, J. Phys. Chem. B 112, 9020 (2008).
  * Shinoda, Shiga, and Mikami (2004) W. Shinoda, M. Shiga, and M. Mikami, Phys. Rev. B 69, 134103 (2004).
  * Tuckerman  _et al._ (2006) M. E. Tuckerman, J. Alejandre, R. López-Rendón, A. L. Jochim, and G. J. Martyna, J. Phys. A 39, 5629 (2006).
  * Hockney and Eastwood (1988) R. W. Hockney and J. W. Eastwood,  _Computer simulation using particles_ (CRC Press, 1988).
  * Kolafa and Perram (1992) J. Kolafa and J. W. Perram, Mol. Sim. 9, 351 (1992).
  * Plimpton (1995) S. Plimpton, J. Comput. Phys. 117, 1 (1995).
  * Wirnsberger  _et al._ (2016) P. Wirnsberger, D. Fijan, A. Šarić, M. Neumann, C. Dellago, and D. Frenkel, J. Chem. Phys. 144, 224102 (2016).
  * Pan, Yi, and Hu (2017) C. Pan, S. Yi, and Z. Hu, Phys. Chem. Chem. Phys. 19, 4861 (2017).
  * Cox and Geissler (2018) S. J. Cox and P. L. Geissler, J. Chem. Phys. 148, 222823 (2018).