# The link between hyperuniformity, Coulomb energy, and Wasserstein distance
to Lebesgue for two-dimensional point processes

Martin Huesmann, Thomas Leblé Institute for Mathematical Stochastics,
University of Münster<EMAIL_ADDRESS>Université de Paris-
Cité, CNRS, MAP5 UMR 8145, F-75006 Paris, France<EMAIL_ADDRESS>
(August 28, 2024)

###### Abstract

We investigate the interplay between three possible properties of stationary
point processes: i) Finite Coulomb energy with short-scale regularization, ii)
Finite $2$-Wasserstein transportation distance to the Lebesgue measure and
iii) Hyperuniformity. In dimension $2$, we prove that i) implies ii), which is
known to imply iii), and we provide simple counter-examples to both converse
implications. However, we prove that ii) implies i) for processes with a
uniformly bounded density of points, and that i) - finiteness of the
regularized Coulomb energy - is equivalent to a certain property of
quantitative hyperuniformity that is just slightly stronger than
hyperuniformity itself.

Our proof relies on the classical link between $H^{-1}$-norm and
$2$-Wasserstein distance between measures, on the screening construction of
[SS15] for Coulomb gases (of which we present an adaptation to $2$-Wasserstein
space which might be of independent interest), and on the recent necessary and
sufficient conditions given in [SWY23] for the existence of stationary
“electric” fields compatible with a given stationary point process.

##  1 Introduction

###  1.1 Setting and main results

Let $\mathsf{d}\geq 1$ be the ambient dimension and let
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ be a
stationary point process, i.e. a random locally finite configuration of points
in $\mathbb{R}^{\mathsf{d}}$ whose distribution is invariant under
translations.

(HU)

For $r>0$, we let $\mathrm{B}_{r}$ be the closed Euclidean ball of radius $r$
centered at the origin and if $\mathbf{X}$ is a point configuration we denote
by $|\mathbf{X}\cap\mathrm{B}_{r}|$ the number of its points in
$\mathrm{B}_{r}$. We say that a point process
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ is
_hyperuniform_ when (see Section 1.2 for more context):

$\sigma(r):=\frac{\mathrm{Var}[|\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}\cap\mathrm{B}_{r}|]}{|\mathrm{B}_{r}|}\longrightarrow
0\text{ as }r\to\infty.$ (1.1)

($\texttt{Wass}_{2}$)

We say that $\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$
is at finite $2$-Wasserstein distance of the Lebesgue measure when there
exists a stationary random coupling between the points of
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ and the
Lebesgue measure, whose quadratic cost per unit volume is bounded on average
(see Section 1.3).

(Coul)

We say that a stationary point process
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ has _finite
regularized Coulomb energy_ when there exists a stationary _electric field_
$\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ compatible
with $\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ (a
random vector field which could be the gradient of the Coulomb potential
generated by the random points of
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$) that
satisfies:
$\mathbb{E}\left[|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\eta}(0)|^{2}\right]<+\infty$,
where $\eta\in(0,1)$ is any short-distance truncation lengthscale. We refer to
Section 1.4 for a precise definition.

There has been recent interest in studying the link between those three
notions when the dimension $\mathsf{d}=2$. In [DFHL24], it is proven that
$\texttt{Wass}_{2}$ implies HU, confirming statements in the statistical
physics literature. In [LRY24], a partial converse is stated as follows: if a
point process is hyperuniform (HU) and if its two-point correlation function
satisfies _“a mild integrability condition”_ then it is at finite distance to
Lebesgue ($\texttt{Wass}_{2}$). We discuss those results in Section 1.5.

We argue here that the integrability condition of [LRY24] essentially
corresponds to the property Coul of having finite regularized Coulomb energy,
and that Coul implies indeed $\texttt{Wass}_{2}$. Moreover, we introduce a
slight reinforcement of the notion of hyperuniformity, which we call
$\texttt{HU}_{\star}$ (see Definition 1.1), and we show that
$\texttt{HU}_{\star}$ implies Coul. In summary, for $\mathsf{d}=2$:

$\texttt{HU}_{\star}\implies\texttt{Coul}\implies\texttt{Wass}_{2}\implies\texttt{HU}.$

Moreover, we prove that the first implication is an equivalence but that the
other two are not.

#### From finite Coulomb energy to finite $2$-Wasserstein distance (…to
hyperuniformity)

Our first result is the implication $\texttt{Coul}\implies\texttt{Wass}_{2}$,
which is in fact true for all $\mathsf{d}\geq 2$.

######  Theorem 1 ($\texttt{Coul}\implies\texttt{Wass}_{2}$).

In any dimension $\mathsf{d}\geq 2$, Coul implies $\texttt{Wass}_{2}$. More
precisely, there exists a constant $C$ depending only on the dimension such
that, for all stationary point processes
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$:

$\mathsf{Wass}^{2}_{2}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}},\mathrm{Leb})\leq
C\left(\mathsf{Coul}_{1}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})+1\right),$
(1.2)

thus effectively bounding the $2$-Wasserstein distance to Lebesgue (defined in
Section 1.3) by the regularized Coulomb energy $\mathsf{Coul}_{1}$ (see
Definition 1.4).

The Wasserstein-Coulomb inequality (1.2) can be thought of as the infinite-
volume equivalent of classical inequalities bounding the $2$-Wasserstein
distance between measures by the Sobolev $H^{-1}$-norm of their difference
(see [Led19, Sec. 2] and the references therein). Combined with the main
result of [DFHL24] this yields (in dimension $2$):

$\texttt{Coul}\implies\texttt{Wass}_{2}\implies\texttt{HU}.$ (1.3)

We prove Theorem 1 in Section 2.

It is natural to ask whether the converse implications hold in (1.3). Thanks
to simple counter-examples, we show that the answer is negative:

######  Theorem 2 ($\texttt{Wass}_{2}\centernot\implies\texttt{Coul}$,
$\texttt{HU}\centernot\implies\texttt{Wass}_{2}$).

In dimension $2$:

  1. 1.

There exists a stationary point process which has finite $2$-Wasserstein
distance to Lebesgue but infinite regularized Coulomb energy. Hence
$\texttt{Wass}_{2}$ does not imply Coul.

  2. 2.

There exists a stationary point process which is hyperuniform yet has infinite
$1$-Wasserstein distance to Lebesgue. In particular, its $2$-Wasserstein
distance is also infinite, and thus HU does not imply $\texttt{Wass}_{2}$.

The constructions leading to Theorem 2 are presented in Section 5.

#### Finite $2$-Wasserstein distance plus density bound implies finite Coulomb
energy

The fact that Coul implies $\texttt{Wass}_{2}$ but not vice-versa is related
to the slight difference in behavior between the Coulomb energy and the
$2$-Wasserstein distance with respect to an abnormal concentration of points
at a given place. In short, the cost of such situations in terms of Coulomb
energy is larger than the cost in $2$-Wasserstein distance by a logarithmic
factor, which gives us enough freedom to construct counter-examples. In a
related direction, it is observed in [DFHL24] that finiteness of the
$(2+\varepsilon)$-Wasserstein distance for some $\varepsilon>0$ is enough to
imply $\texttt{HU}_{\star}$ and thus, by our results, $\mathsf{Coul}$.

However, if we assume that the points of
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ _never_
gather too much, in the sense of a uniform, almost sure bound on the local
density of points for
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$, then we
get a converse statement.

######  Theorem 3 ($\texttt{Wass}_{2}$ plus density bound implies Coul).

Let $\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ be a
two-dimensional stationary point process satisfying $\texttt{Wass}_{2}$.
Assume that there exists $\bar{\rho}\geq 1$ such that ($\mathrm{B}_{1}(x)$
being the unit disk around $x$):

$\sup_{x\in\mathbb{R}^{\mathsf{d}}}\left|\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}\cap\mathrm{B}_{1}(x)\right|\leq\bar{\rho}\text{
almost surely.}$ (1.4)

Then Coul is satisfied and in fact, for some constant $C$ not depending on
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$, we have:

$\mathsf{Coul}_{1}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})\leq
C\bar{\rho}\left(\mathsf{Wass}^{2}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}},\mathrm{Leb})+1\right)$
(1.5)

Clearly, $B_{1}(x)$ can be replaced in (1.4) by $B_{r}(x)$ for any fixed $r$.
We prove Theorem 3 in Section 3. It can be thought of as the infinite-volume
version of inequalities relating the $H^{-1}$-norm of the difference of two
measures to their $2$-Wasserstein distance under certain density assumptions
as in [Loe06, Prop. 3.1].

#### From a slight reinforcement of hyperuniformity to finite Coulomb energy

In dimension $2$, we have $\mathsf{Coul}\implies\texttt{HU}$ (via
$\texttt{Wass}_{2}$, see (1.3)) but some point processes that are only very
weakly hyperuniform (in the sense that their rescaled number variance
$\sigma(r)$ goes to $0$ very slowly as $r\to\infty$) do not have finite
Coulomb energy, see Theorem 2. In order to formulate a converse statement, we
introduce the following property, which is slightly stronger than
hyperuniformity itself.

######  Definition 1.1 ($\texttt{HU}_{\star}$).

We say that a stationary point process satisfies $\texttt{HU}_{\star}$ when
the rescaled number variance $\sigma(r)$ (see (1.1)) not only goes to $0$ as
$r\to\infty$ but does so _in such a way that_ :

$\sum_{n\geq 0}\sigma(2^{n})<+\infty.$ (1.6)

Note that if we know a priori that $r\mapsto\sigma(r)$ is non-increasing, then
convergence of the series implies hyperuniformity in the usual sense, however
we cannot rule out non-monotonous behaviors for $\sigma$ hence we prefer to
define $\texttt{HU}_{\star}$ as HU _plus_ said convergence.

We call condition (1.6) a “slight reinforcement of hyperuniformity” because
any point process such that $\sigma(r)\to 0$ as some negative power of $r$, or
even as $(\log r)^{-c}$ with $c>1$, satisfies $\texttt{HU}_{\star}$. However,
it is a fact that $\texttt{HU}_{\star}$ _is strictly stronger than_ HU since
there exists point processes such that $\lim_{r\to\infty}\sigma(r)=0$ but
$\liminf_{r\to\infty}\sigma(r)\log r>0$, see [DFHL24, Thm. 3] \- amusingly,
such a behavior is never discussed in the statistical physics literature about
hyperuniformity, presumably because those decays are so slow that they are
deemed “unphysical”.

In dimension $2$, we find that making the ever-so-slightly stronger assumption
$\texttt{HU}_{\star}$ instead of HU is enough to deduce Coul \- in fact
$\texttt{HU}_{\star}$ is _equivalent_ to finiteness of the regularized Coulomb
energy.

######  Theorem 4 ($\texttt{HU}_{\star}\iff\texttt{Coul}$).

Let $\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ be a
two-dimensional stationary point process.

  1. 1.

If $\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$
satisfies $\texttt{HU}_{\star}$, then it has finite regularized Coulomb energy
i.e. satisfies Coul. More precisely, there exists a constant $C$ (not
depending on
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$) such that:

$\mathsf{Coul}_{1}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})\leq
C\left(\sum_{n\geq 0}\sigma(2^{n})+1\right).$ (1.7)

  2. 2.

Conversely, Coul implies $\texttt{HU}_{\star}$.

The proof of Theorem 4 relies on a recent result by Sodin-Wennman-Yakir and
goes by showing that $\texttt{HU}_{\star}\iff\texttt{SC}\iff\texttt{Coul}$,
where SC is a “spectral condition” introduced in [SWY23], see Section 4.

######  Remark 1.2.

Our proof of the implication $\texttt{Coul}\implies\texttt{HU}_{\star}$ uses
results which are not explicitly written down in a quantitative way in
[SWY23]. However, an inspection of [SWY23] together with our Lemmas 4.4 and
4.7 suggests that one could write the implication
$\texttt{Coul}\implies\texttt{HU}_{\star}$ as an effective inequality of the
form:

$\sum_{n\geq 0}\sigma(2^{n})\leq
C\left(\mathsf{Coul}_{1}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})+1\right).$

#### Applications

We now mention some consequences of our results (in dimension $2$).

  * •

Since the Poisson point process is not hyperuniform, it must have infinite
regularized Coulomb energy, and thus also infinite renormalized Coulomb energy
in the sense of Sandier-Serfaty (see Section 1.4), which answers a question
raised in [LS17, Sec. 3.2]. This could also be deduced directly from the
results of [SWY23].

  * •

Following the phenomenological classification of Torquato [Tor18], all
hyperuniform point processes of class-I (i.e. $\sigma(r)=\mathcal{O}(r^{-1})$)
, class-II (i.e. $\sigma(r)=\mathcal{O}(\log r\times r^{-1})$) and class-III
(i.e. $\sigma(r)=\mathcal{O}(r^{-\alpha})$ with $\alpha\in(0,1)$) satisfy
$\texttt{HU}_{\star}$ and thus have both finite regularized Coulomb energy and
finite $2$-Wasserstein distance to the Lebesgue measure. In particular, they
can all be written as the perturbation of a stationary lattice by stationary
displacements with finite second moment - in total agreement with the phrasing
of [LRY24]: _“most planar hyperuniform point processes are $L^{2}$-perturbed
lattices”_. Of course, by the triangular inequality, this also means that they
are at finite $2$-Wasserstein distance of each other.

  * •

Minimizers of the $2d$ Coulomb gas free energy functional introduced in [LS17]
all have finite renormalized Coulomb energy (by definition), hence they have
finite regularized Coulomb energy (which is a less stringent condition), and
are thus hyperuniform - in fact, they even satisfy $\texttt{HU}_{\star}$.

Note, however, that those point processes describe the limiting behavior of
objects named “empirical fields” in the statistical mechanics literature,
which are defined as “averages” (over translations) of the finite $2d$ Coulomb
gas. It is not at all obvious whether they coincide or not with the actual
limit points (in the space of point processes) of the finite random point
configuration corresponding to the Gibbs measure of a $2d$ Coulomb gas
_without spatial averaging_ (the mere existence of those limit points is hard
to prove see [AS21, Cor. 1.1]). Hyperuniformity in this case has been recently
proven [Leb21] \- interestingly, the upper bound on the rescaled number
variance given there is not enough to ensure $\texttt{HU}_{\star}$ and
finiteness of the Coulomb energy for those limit points is yet an open
question (although a positive answer is very plausible).

  * •

Hardcore point processes, having an almost sure uniform lower bound on the
distance between two points (one can e.g. obtain such systems by considering
Gibbs point processes with hardcore interactions), satisfy a density bound of
the form (1.4), and our Theorem 3 applies to them.

###  1.2 Hyperuniformity

The notion of _hyperuniform_ or _super-homogeneous_ distributions of points
originated in the theoretical chemistry and statistical physics literature
twenty years ago [GJL02, TS03]. We refer to the more recent survey by S.
Torquato [Tor18] for many examples of hyperuniform processes and for an
overview of the considerable literature. A simple way to think of the
Definition (1.1) is to say that for a random point configuration to be
hyperuniform means that its number variance (in large balls) is asymptotically
much smaller than that of a Poisson point process.
One reason that makes hyperuniform systems interesting is that they can
exhibit a form of “order within disorder”. Indeed, the smallness of the number
variance expressed by (1.1) is a form of order which is not necessarily coming
from any long-range persistence of correlations within the system. For
instance:

  1. 1.

The Poisson point process: is not hyperuniform, has absolutely no
correlations, it represents “disorder”.

  2. 2.

The stationary lattice: is “as hyperuniform as it gets”, possesses strong
correlations, it represents “order”.

  3. 3.

The Ginibre ensemble (a two-dimensional point process arising from non-
Hermitian random matrix theory [Gin65]): is “as hyperuniform” as a lattice and
yet has a (truncated) $2$-point correlation function that decays faster than
any exponential, it represents “order within disorder”.

###  1.3 Wasserstein distances for stationary point processes

We start by a brief reminder on some notions of “finite-volume” optimal
transport that will be useful for us. This is by now a classical topic for
which there exists numerous expository texts, e.g. [San15, Vil09, FG23]. We
then recall the definition (and some properties) of Wasserstein distances in
_infinite-volume_ for stationary point processes as recently developed in
[EHJM23].

#### Wasserstein distance in finite volume

Let $\mathbf{m}_{0},\mathbf{m}_{1}$ be two measures on
$\mathbb{R}^{\mathsf{d}}$ with the same finite mass. We denote by
$\mathsf{cpl}(\mathbf{m}_{0},\mathbf{m}_{1})$ the set of all _couplings_
between $\mathbf{m}_{0}$ and $\mathbf{m}_{1}$, i.e. the set of all measures on
$\mathbb{R}^{\mathsf{d}}\times\mathbb{R}^{\mathsf{d}}$ having
$\mathbf{m}_{0},\mathbf{m}_{1}$ as marginals. We recall that for
$p\in[1,+\infty)$, the “usual” $p$-Wasserstein distance between
$\mathbf{m}_{0}$ and $\mathbf{m}_{1}$ is defined as:

$\mathrm{W}_{p}(\mathbf{m}_{0},\mathbf{m}_{1}):=\left(\inf_{\Pi\in\mathsf{cpl}(\mathbf{m}_{0},\mathbf{m}_{1})}\int_{\mathbb{R}^{\mathsf{d}}\times\mathbb{R}^{\mathsf{d}}}|x-y|^{p}\mathtt{d}\Pi(x,y)\right)^{1/p}.$
(1.8)

We will mostly deal here with the $2$-Wasserstein distance, associated to the
“quadratic cost” $|x-y|^{2}$.

#### The Benamou-Brenier formulation of optimal transport

Let $\Lambda\subset\mathbb{R}^{\mathsf{d}}$ be an hypercube (or in fact any
convex set, or the whole $\mathbb{R}^{\mathsf{d}}$) and let $\mathbf{m}_{0}$
and $\mathbf{m}_{1}$ be two probability measures on $\Lambda$ possessing
bounded densities $\mathrm{m}_{0},\mathrm{m}_{1}$ with respect to the Lebesgue
measure. We call “an admissible pair” any family
$(\mathbf{m}_{t},v_{t})_{t\in[0,1]}$ (indexed by the “time” parameter $t$),
where $t\mapsto\mathbf{m}_{t}$ is an absolutely continuous family of
probability measures on $\Lambda$ interpolating between $\mathbf{m}_{0}$ and
$\mathbf{m}_{1}$, with bounded densities $t\mapsto\mathrm{m}_{t}$, and
$t\mapsto v_{t}$ is a measurable family of vector fields on
$\mathbb{R}^{\mathsf{d}}$ (the _velocity fields_), such that:

  1. 1.

The pair $(\mathbf{m}_{t},v_{t})$ solves the _continuity equation_

$\frac{\mathtt{d}}{\mathtt{d}t}\mathrm{m}_{t}+\mathrm{div}(\mathrm{m}_{t}v_{t})=0$
(1.9)

in the sense of distributions, i.e. for all smooth compactly supported test
function $f:(0,1)\times\mathbb{R}^{\mathsf{d}}\to\mathbb{R}$ we have:

$\int_{0}^{1}\left(\int_{\mathbb{R}^{\mathsf{d}}}\partial_{t}f(t,x)\mathrm{m}_{t}(x)\mathtt{d}x\right)\mathtt{d}t+\int_{0}^{1}\left(\int_{\mathbb{R}^{\mathsf{d}}}\nabla_{x}f(t,x)\cdot
v_{t}(x)\mathrm{m}_{t}(x)\mathtt{d}x\right)\mathtt{d}t=0.$

  2. 2.

The Neumann condition $\mathrm{m}_{t}v_{t}\cdot\vec{n}=0$ is satisfied on
$\partial\Lambda$ for all $t\in[0,1]$.

  3. 3.

The density $\mathrm{m}_{t}$ satisfies (see [Vil09, Corollary 19.5])

$\sup_{t\in[0,1]}\|\mathrm{m}_{t}\|_{L^{\infty}}\leq\max\left(\|\mathrm{m}_{0}\|_{\infty},\|\mathrm{m}_{1}\|_{\infty}\right).$
(1.10)

The _action_ of such a pair is then defined by:

$\mathcal{A}\left((\mathbf{m}_{t},v_{t})_{t}\right):=\int_{0}^{1}\left(\int_{\Lambda}|v_{t}|^{2}\mathtt{d}\mathbf{m}_{t}\right)\mathtt{d}t.$
(1.11)

The celebrated Benamou-Brenier formula of [BB00] (see [San15, Sec. 5.3] for a
presentation) asserts that the Wasserstein distance between $\mathbf{m}_{0}$
and $\mathbf{m}_{1}$ can be achieved as the infimum of the action over
admissible pairs:

$\mathrm{W}^{2}_{2}(\mathbf{m}_{0},\mathbf{m}_{1})=\inf_{(\mathbf{m}_{t},v_{t})_{t}\text{
admissible}}\mathcal{A}\left((\mathbf{m}_{t},v_{t})_{t}\right).$ (1.12)

Moreover, the continuity equation and Benamou-Brenier formula localise to
convex subsets of $\Lambda$ and $m_{t}v_{t}$ admit internal (and external)
traces on $\partial B_{r}$ (for $r>0$) such that we can safely consider
$m_{t}v_{t}$ on $\partial B_{r}$, see [GHO21, Sec. 2.3].

#### Optimal transport and Wasserstein distances for stationary point
processes

In order to define a notion of Wasserstein distance between invariant point
processes or, more generally, invariant random measures, we introduce suitable
cost functions, resembling a Wasserstein distance themselves:

  * •

For two Radon measures $\mathrm{M}_{0},\mathrm{M}_{1}$ we let
$\mathsf{cpl}(\mathrm{M}_{0},\mathrm{M}_{1})$ be the set of all Radon measures
on $\mathbb{R}^{\mathsf{d}}\times\mathbb{R}^{\mathsf{d}}$ with marginals
$\mathrm{M}_{0}$ and $\mathrm{M}_{1}$. Then we define the cost function
$\mathsf{w}_{p}$ as a “Wasserstein cost per unit volume” between $M_{0}$ and
$M_{1}$:

$\mathsf{w}_{p}(\mathrm{M}_{0},\mathrm{M}_{1}):=\left(\inf_{q\in\mathsf{cpl}(\mathrm{M}_{0},\mathrm{M}_{1})}\limsup_{n\to\infty}\frac{1}{n^{\mathsf{d}}}\int_{\Lambda_{n}\times\mathbb{R}^{2}}|x-y|^{p}\mathtt{d}q(x,y)\right)^{{\frac{1}{p}}}.$
(1.13)

  * •

For laws $\mathsf{P}_{0},\mathsf{P}_{1}$ of shift-invariant random measures,
i.e. stationary probability measures on the space of Radon measures, we define
(as in [EHJM23]) their $p$-Wasserstein distance by

$\mathsf{Wass}_{p}(\mathsf{P}_{0},\mathsf{P}_{1}):=\left(\inf_{Q\in\mathsf{Cpl}_{s}(\mathsf{P}_{0},\mathsf{P}_{1})}\int\mathsf{w}^{p}_{p}(\mathrm{M}_{0},\mathrm{M}_{1})\
Q(d\mathrm{M}_{0},d\mathrm{M}_{1})\right)^{\frac{1}{p}},$ (1.14)

where $\mathsf{Cpl}_{s}(\mathsf{P}_{0},\mathsf{P}_{1})$ denotes the set of all
couplings between $\mathsf{P}_{0}$ and $\mathsf{P}_{1}$ which are stationary
under the diagonal action of $\mathbb{R}^{\mathsf{d}}$, i.e.
$Q=Q\circ(\theta_{z},\theta_{z})^{-1}$ with
$(\theta_{z}\mathrm{M})(A)=\mathrm{M}(A+z)$ for all Radon measures
$\mathrm{M}$, Borel sets $A\subset\mathbb{R}^{\mathsf{d}}$ and
$z\in\mathbb{R}^{\mathsf{d}}$.

It is shown in [EHJM23], that $\mathsf{Wass}_{p}$ defines a geodesic metric
between invariant random measures of equal intensity.

We are particularly interested in the distance between (the law of a) random
point configuration
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$, or of its
“spread-out” version
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\eta}$ (see
(1.24)) to either the Lebesgue measure $\mathrm{Leb}$ (seen as a deterministic
Radon measure) or to the “stationary lattice”
$\mathsf{P}^{\mathbb{Z}^{\mathsf{d}}}$ defined as (the law of)
$\sum_{x\in\mathbb{Z}^{\mathsf{d}}}\delta_{x+\tau}$ (with a
$[0,1]^{\mathsf{d}}$-uniform random variable $\tau$). We will often make a
slight abuse of notation and write
$\mathsf{Wass}_{p}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}},\mathrm{Leb})$
or
$\mathsf{Wass}_{p}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}},\mathsf{P}^{\mathbb{Z}^{\mathsf{d}}})$
although $\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ is
a random measure, $\mathrm{Leb}$ is a deterministic measure, and
$\mathsf{P}^{\mathbb{Z}^{\mathsf{d}}}$ _the law of_ a random measure. The
definition of
$\mathsf{Wass}_{p}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}},\mathrm{Leb})$
reduces to:

$\mathsf{Wass}_{p}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}},\mathrm{Leb})=\mathbb{E}\left[\mathsf{w}^{p}_{p}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}},\mathrm{Leb})\right]^{1/p},$
(1.15)

with $\mathsf{w}_{p}$ the $p$-cost “per unit volume” between a realization of
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ and the
Lebesgue measure as in (1.13). Let us make two additional comments:

  * •

The perturbed lattice and the Lebesgue measure are at finite $p$-Wasserstein
distance of each other for all $p$ (in fact for $p=+\infty$) hence, by the
triangular inequality, it is equivalent to be at finite distance of either.

  * •

Thanks to the alternative representation of [EHJM23, Prop. 2.11.] of the
Wasserstein distance between point processes as the cost of optimal
_matchings_ , the condition
$\mathsf{Wass}_{p}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}},\mathsf{P}^{\mathbb{Z}^{\mathsf{d}}})<+\infty$
is equivalent to the fact that
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ can be
written as a _perturbed lattice_ i.e. as the point process
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}=\sum_{x\in\mathbb{Z}^{\mathsf{d}}}\delta_{x+v_{x}}$
where the displacements $(v_{x})_{x\in\mathbb{Z}^{\mathsf{d}}}$ are random
perturbations whose law is invariant under shifts in $\mathbb{Z}^{\mathsf{d}}$
and such that _the $p$-th moment of $v_{0}$ is finite._

##### Summary.

We end this section with a quick summary of the various Wasserstein distances
and the corresponding notation:

  * •

$\mathrm{W}_{p}$ denotes the usual $p$-Wasserstein distance between two finite
measures on $\mathbb{R}^{n}$ (see (1.8)).

  * •

$\mathsf{w}_{p}$ is the $p$-Wasserstein distance _per unit volume_ between two
point configurations, or more generally infinite Radon measures (see (1.13)).

  * •

$\mathsf{Wass}_{p}$ is the $p$-Wasserstein distance between (the law of) two
random point configurations (or two random Radom measures), defined in (1.14).

###  1.4 Infinite-volume Coulomb energy

We let $\mathsf{g}$ be (a certain multiple of) the Coulomb kernel in dimension
$\mathsf{d}$, namely:

$\mathsf{g}(x):=\begin{cases}-\log|x|&(\mathsf{d}=2)\\\
|x|^{-(\mathsf{d}-2)}&(\mathsf{d}\geq 3)\end{cases}.$

#### From finite to infinite systems

We recall the basic setup for (one-component) Coulomb systems in finite
volume. Let $N\geq 1$ and let $\mathrm{X}_{N}=(x_{1},\dots,x_{N})$ be a
$N$-tuple of points in the box
$\Lambda_{N}:=\left[-\frac{N}{2},\frac{N}{2}\right]^{\mathsf{d}}$. The Coulomb
energy of $\mathrm{X}_{N}$ can be defined as:

$\mathrm{H}_{N}(\mathrm{X}_{N}):=\frac{1}{2}\iint_{(x,y)\in\Lambda_{N}\times\Lambda_{N},x\neq
y}\mathsf{g}(x-y)\mathtt{d}\left(\mathbf{X}_{N}-\mathrm{Leb}\right)(x)\mathtt{d}\left(\mathbf{X}_{N}-\mathrm{Leb}\right)(y),$
(1.16)

with $\mathbf{X}_{N}:=\sum_{i=1}^{N}\delta_{x_{i}}$ a finite point
configuration. The usual physical interpretation of (1.16) is that the $N$
point charges $x_{1},\dots,x_{N}$, all living in $\Lambda_{N}$, interact both
with each other and with a “neutralizing background” given by the Lebesgue
measure on the box $\Lambda_{N}$. The energy $\mathrm{H}_{N}$ is bounded below
by a constant times $N$, but it is not bounded above due to the singularity of
$\mathsf{g}$.

The study of systems with pairwise Coulomb interactions is an old topic.
However, due to the long-range nature of the Coulomb kernel, properly defining
the Coulomb energy (per unit volume) for an infinite point configuration, or
for stationary point processes is challenging.

#### Infinite systems

The notion of a “renormalized energy” was introduced by Sandier-Serfaty
[SS12], and applied to two-dimensional Coulomb systems in [SS15], with various
extensions (see in particular [RS16] for higher dimensional Coulomb gases and
[PS17] for an generalization to Riesz gases). It provides a/the framework to
work with infinite Coulomb systems.

Before presenting this formalism, let us briefly evoke what a “naive” approach
could look like. In view of the finite interaction energy (1.16) and of
standard practice in statistical mechanics (for _short-range_ interactions!)
one could (try to) define the Coulomb energy per unit volume of a point
configuration $\mathbf{X}$ as:

$\lim_{N\to\infty}\frac{1}{|\Lambda_{N}|}\
\frac{1}{2}\iint_{(x,y)\in\Lambda_{N}\times\Lambda_{N},x\neq
y}\mathsf{g}(x-y)\mathtt{d}\left(\mathbf{X}-\mathrm{Leb}\right)(x)\mathtt{d}\left(\mathbf{X}-\mathrm{Leb}\right)(y),$
(1.17)

or perhaps more accurately as (note the change in the domain of integration):

$\lim_{N\to\infty}\frac{1}{|\Lambda_{N}|}\
\frac{1}{2}\iint_{(x,y)\in\Lambda_{N}\times\mathbb{R}^{\mathsf{d}},x\neq
y}\mathsf{g}(x-y)\mathtt{d}\left(\mathbf{X}-\mathrm{Leb}\right)(x)\mathtt{d}\left(\mathbf{X}-\mathrm{Leb}\right)(y).$
(1.18)

The second expression is in fact closer to the truth, as it takes into account
the effect of the whole configuration on the points in $\Lambda_{N}$, but
precisely because of this, and the long-range nature of the Coulomb kernel,
the finiteness of this expression for fixed $N$ is very much unclear and might
only be obtained through a careful control of the overall repartition of
charges. Ultimately, it is simply not possible in general to rigorously derive
“infinite-volume” energies of the type (1.18) from the finite $N$ energy
(1.16).

The general spirit of the works by Serfaty et al. can be phrased as an
“electric approach”. Instead of working directly with the pairwise interaction
energy $\mathsf{g}(x-y)$ between points as in (1.17), (1.18), one puts the
emphasis on “electric fields”, which are intrinsically global objects whose
energy is nonetheless easy to localize. Near each point of the configuration,
those fields have a singularity which needs to be taken care of, hence the
adjective “renormalized”.

In the rest of Section 1.4, we recall the main lines of the “electric
approach”, our goal being to define the “regularized Coulomb energy”
$\mathsf{Coul}_{\eta}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})$
($\eta\in(0,1)$) for a stationary point process
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$, which is
presented in Definition 1.4.

#### Electric approach in finite volume

The starting point of the “electric” approach to _infinite_ Coulomb systems is
the following consideration valid for _finite_ systems. Let us introduce the
electric potential $\Phi$ and the _true electric field_ $\mathrm{E}$ generated
by $\mathrm{X}_{N}$ as:

$\Phi:x\mapsto\int_{y\in\Lambda_{N}}\mathsf{g}(x-y)\mathtt{d}\left(\mathbf{X}_{N}-\mathrm{Leb}\right)(y),\quad\mathrm{E}(x):=\nabla\Phi=\int_{y\in\Lambda_{N}}\nabla\mathsf{g}(x-y)\mathtt{d}\left(\mathbf{X}_{N}-\mathrm{Leb}\right)(y).$
(1.19)

Note that $\Phi$ is a scalar field, while $\mathrm{E}$ is a (gradient) vector
field. The “true electric field” satisfies the distributional identity:

$-\mathrm{div}(\mathrm{E})=c_{\mathsf{d}}\left(\mathbf{X}_{N}-\mathrm{Leb}\right),\text{
with }\mathrm{E}(x)\to 0\text{ as }|x|\to\infty$ (1.20)

for some constant $c_{\mathsf{d}}$ depending on $\mathsf{d}$ (the value is
$2\pi$ for $\mathsf{d}=2$). Returning to (1.16), integrating by parts, and
neglecting the condition $\\{x\neq y\\}$ in the integral, one would get:

$\mathrm{H}_{N}(\mathrm{X}_{N})\approx-\frac{1}{2c_{\mathsf{d}}}\int_{\mathbb{R}^{\mathsf{d}}}\Phi\
\mathrm{div}(\mathrm{E})\approx\frac{1}{2c_{\mathsf{d}}}\int_{\mathbb{R}^{\mathsf{d}}}|\mathrm{E}|^{2}.$
(1.21)

The heuristic identity (1.21) does not make sense as such because $\mathrm{E}$
fails to be in $L^{2}$ due to the point charges self-interactions (which are
avoided in (1.16) by imposing the condition $\\{x\neq y\\}$). However, (1.21)
points to the fact (which is old and common knowledge in physics) that the
Coulomb interaction energy can be phrased in terms of the integral of an
“electric” energy density $|\mathrm{E}|^{2}$, where $\mathrm{E}$ is the right
“electric field” satisfying (1.20).

#### Electric fields: global, local, screened

We refer to vector fields in
$\cap_{p\in[1,2)}L^{p}_{loc}(\mathbb{R}^{\mathsf{d}},\mathbb{R}^{\mathsf{d}})$
as (global) “electric fields”. If $\mathbf{X}$ is a point configuration in
$\mathbb{R}^{\mathsf{d}}$, and $\mathrm{E}$ is an electric field, we say that
$\mathrm{E}$ is compatible with $\mathbf{X}$ when the following identity
holds:

$-\mathrm{div}(\mathrm{E})=c_{\mathsf{d}}\left(\mathbf{X}-\mathrm{Leb}\right).$
(1.22)

We also work with _local_ electric fields: if $\Lambda$ is a square in
$\mathbb{R}^{\mathsf{d}}$ and $\mathbf{X}$ a point configuration in $\Lambda$,
we say that
$\mathrm{E}\in\cap_{p\in[1,2)}L^{p}_{loc}(\Lambda,\mathbb{R}^{\mathsf{d}})$ is
a (local) electric field compatible with $\mathbf{X}$ when (1.22) holds in
$\Lambda$. To be clear, by (1.22) we mean that

$\int\varphi(x)\mathtt{d}\left(\mathbf{X}-\mathrm{Leb}\right)(x)=-\frac{1}{c_{\mathsf{d}}}\int\nabla\varphi\cdot\mathrm{E}$

for all $\varphi\in C^{\infty}_{c}(\mathbb{R}^{\mathsf{d}})$ (global electric
fields), or for all $\varphi\in C^{\infty}_{c}(\Lambda)$ (local electric
fields in $\Lambda$). We say that a local electric field $\mathrm{E}$ is
_screened_ when it has Neumann boundary conditions i.e. its normal component
satisfies

$\mathrm{E}\cdot\vec{n}=0\text{ along $\partial\Lambda$}.$ (1.23)

###### The importance of screened fields

Screened fields are of outmost importance because they can be glued together
in the following sense: if $\Lambda^{1},\Lambda^{2}$ are two squares with one
side in common, $\mathbf{X}^{1},\mathbf{X}^{2}$ are two point configurations
in $\Lambda^{1},\Lambda^{2}$ and $\mathrm{E}^{1},\mathrm{E}^{2}$ are two local
electric fields compatible with $\mathbf{X}^{1},\mathbf{X}^{2}$, then
defining:

$\mathbf{X}^{t}:=\mathbf{X}^{1}+\mathbf{X}^{2},\quad\mathrm{E}^{t}:=\mathrm{E}^{1}\mathsf{1}_{\Lambda_{1}}+\mathrm{E}^{2}\mathsf{1}_{\Lambda_{2}}$

yields a “total” point configuration and a “total” electric field, but _in
general, $\mathrm{E}^{t}$ is not compatible with $\mathbf{X}^{t}$_. However,
if $\mathrm{E}^{1},\mathrm{E}^{2}$ are both screened, or more generally if
_their normal components agree on the common side of
$\Lambda_{1},\Lambda_{2}$_, then $\mathrm{E}^{t}$ is compatible with
$\mathbf{X}^{t}$.

###### The role of abstract fields

An important remark is that we deal with “abstract” fields $\mathrm{E}$, not
necessarily of the form (cf. (1.19)):

$\mathrm{E}(x)\overset{?}{=}\int_{\mathbb{R}^{\mathsf{d}}}\nabla\mathsf{g}(x-y)\mathtt{d}\left(\mathbf{X}-\mathrm{Leb}\right)(y),$

which could be a candidate to solve (1.22), but whose convergence is unclear
(note that $\nabla\mathsf{g}(x)\propto\frac{1}{|x|^{\mathtt{d}-1}}$ is far
from being integrable at infinity). Considering the whole family of compatible
electric fields is a fundamental feature of the approach by Sandier-Serfaty.
The price to pay is that in general we do not know “who $\mathrm{E}$ is”,
besides that it solves (1.22).

#### Spreading out charges and regularizing point configurations.

To handle the singularity of the electric field near a point charge, one
introduces a short-scale regularization, or “truncation” of the field, and a
corresponding operation of “spreading out” the point charges. Let us start
with the latter, as it is simpler to describe. The usual way to proceed is as
follows (see [RS16, Sec. 2.1]):

  1. 1.

Take some radial non-negative “bump function” $\chi$ supported on the unit
ball, with mass $1$, and for $\eta>0$ let

$\chi_{\eta}:=\frac{1}{\eta^{\mathsf{d}}}\chi\left(\frac{\cdot}{\eta}\right).$

In general, the choice of the bump function $\chi$ is not really important and
does not affect the limiting objects or the values of the limiting quantities
as $\eta\to 0$. A common choice is to take $\chi=\mathsf{1}_{\mathrm{B}_{1}}$,
which would work for us here, but we may also take $\chi$ to be smooth. We fix
such a $\chi$ for the rest of the paper.

  2. 2.

If $x\in\mathbb{R}^{\mathsf{d}}$, define $\delta^{(\eta)}_{x}$ as a Dirac mass
“spread-out” according to $\chi_{\eta}$, namely as the measure with density
$\chi_{\eta}(\cdot-x)$. If we think of a Dirac mass $\delta_{x}$ as a point
charge placed at $x$, then $\delta^{(\eta)}_{x}$ is a “spread-out” charge at
$x$.

This operation goes back at least to Onsager [Ons39], see e.g. [LL01, Sec.
9.7] for a study of its properties in the context of Coulomb energies. In
short:

    1. (a)

By Newton’s theorem, the field generated by a spread-out charge coincides with
the field generated by the original charge outside the support of the “spread”
(here $\mathrm{B}_{\eta}$).

    2. (b)

Spreading out all charges slightly decreases the pairwise interaction energy
between them.

    3. (c)

After this operation, a single charge self-interaction becomes _finite_ \- but
it blows up proportionally to $\mathsf{g}(\eta)$ as $\eta\to 0$, as one could
expect.

  3. 3.

If $\mathbf{X}$ is a point configuration, replacing each Dirac mass
$\delta_{x}$ by its “spread-out” version $\delta^{(\eta)}_{x}$ for
$x\in\mathbf{X}$ yields a “regularized” version of $\mathbf{X}$, denoted by
$\mathbf{X}_{\eta}$.

Formally, we have defined two unambiguous operations:

$\delta_{x}\to\delta^{(\eta)}_{x}:=\chi_{\eta}\ast\delta_{x}\text{ (for
$x\in\mathbb{R}^{\mathsf{d}}$)},\quad\mathbf{X}\to\mathbf{X}_{\eta}:=\sum_{x\in\mathbf{X}}\delta^{(\eta)}_{x}.$
(1.24)

Let us emphasize that $\mathbf{X}\to\mathbf{X}_{\eta}$ is an operation at the
level of Radon measures on $\mathbb{R}^{\mathsf{d}}$. The image of a point
configuration $\mathbf{X}$ is no longer purely atomic, it is a locally finite
measure with locally bounded density with respect to the Lebesgue measure (the
density of one regularized Dirac is bounded by
$\|\chi\|_{\infty}\eta^{-{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathsf{d}}}$
but the density of $\mathbf{X}_{\eta}$ depends also on the number of points
per unit volume).

##### Truncation of the field.

We now study a corresponding operation at the level of electric fields. With
$\chi$ chosen as before and for $\eta>0$, let $\mathsf{f}_{\eta}$ be the mean-
zero solution to:

$-\Delta\mathsf{f}_{\eta}=c_{\mathsf{d}}\left(\delta^{(\eta)}_{0}-\delta_{0}\right),\quad\mathsf{f}_{\eta}\equiv
0\text{ outside }\mathrm{B}_{\eta}.$ (1.25)

If $\mathrm{E}$ is an electric field (global or local) compatible with
$\mathbf{X}$, we define the “truncated field” $\mathrm{E}_{\eta}$ as:

$\mathrm{E}_{\eta}:=\mathrm{E}+\sum_{x\in\mathbf{X}}\nabla\mathsf{f}_{\eta}(x-\cdot).$
(1.26)

######  Claim 1.3 (Properties of the truncated field).

It is straightforward that:

  1. 1.

If $\mathrm{E}$ is compatible with $\mathbf{X}$ in the sense of (1.22), then
$\mathrm{E}_{\eta}$ is compatible with $\mathbf{X}_{\eta}$ in the sense:

$-\mathrm{div}(\mathrm{E}_{\eta})=c_{\mathsf{d}}\left(\mathbf{X}_{\eta}-\mathrm{Leb}\right).$
(1.27)

  2. 2.

If $\mathrm{E}$ is a gradient vector field, then so is $\mathrm{E}_{\eta}$.

  3. 3.

(Local version.) If $\mathrm{E}$ is compatible with $\mathbf{X}$ in $\Lambda$
and $\min_{x\in\mathbf{X}}\mathrm{dist}(x,\partial\Omega)>\eta$, then (1.27)
holds in $\Lambda$. Moreover, if $\mathrm{E}$ is “screened”, then so is
$\mathrm{E}_{\eta}$.

#### Stationary random electric fields

We call a random electric field
$\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ stationary
when its law is invariant under translations
$\mathrm{E}\to\mathrm{E}(\cdot-t)$ for all $t\in\mathbb{R}^{\mathsf{d}}$. We
say that a random electric field
$\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ is
compatible with a point process
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ when the
random measure
$\frac{1}{c_{\mathsf{d}}}\mathrm{div}(\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})+\mathrm{Leb}$
and $\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ have
the same distribution. Note however that if
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ is
stationary, it does not imply that so is
$\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$.

#### Regularized Coulomb energy

Finally, we may define the notion that will be relevant for us, namely “finite
regularized Coulomb energy”.

######  Definition 1.4.

Let $\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ be a
stationary point process. We say that it has _finite regularized Coulomb
energy_ when there exists $\eta\in(0,1]$ and a _stationary_ random electric
field $\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$
compatible with
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ such that,
using the truncation/regularization presented in (1.26):

$\mathbb{E}\left[|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\eta}(0)|^{2}\right]<+\infty.$
(1.28)

Since stationary, compatible electric fields are not unique, we define:

$\mathsf{Coul}_{\eta}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}):=\inf_{\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}\text{
stat., comp. with
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$}}\mathbb{E}\left[|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\eta}(0)|^{2}\right].$
(1.29)

######  Lemma 1.5 (Two properties of the regularized energy).

The following holds:

  1. 1.

If
$\mathsf{Coul}_{\eta}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})$
is finite for some $\eta\in(0,1]$, then the number of points of
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ in
$\mathrm{B}_{r}$ has a finite second moment for all $r>0$ and more precisely:

$\mathbb{E}\left[\left(|\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}\cap\mathrm{B}_{r}|-|\mathrm{B}_{r}|\right)^{2}\right]\leq
C\left(\mathsf{Coul}_{\eta}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})+1\right)r^{2}.$
(1.30)

  2. 2.

If
$\mathsf{Coul}_{1}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})$
is finite, then so is
$\mathsf{Coul}_{\eta}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})$
for all $\eta\in(0,1)$ and we have:

$\mathsf{Coul}_{\eta}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})\leq
C_{\eta}\left(\mathsf{Coul}_{1}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})+1\right)$
(1.31)

###### Proof of Lemma 1.5.

The first item (control of the discrepancies by the Coulomb energy) was proven
in [LS17, Sec. 8.5]. Now, in view of (1.26), we can control the difference
$\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{1}(0)-\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\eta}(0)$
almost surely by:

$\left|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{1}(0)-\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\eta}(0)\right|\leq|\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}\cap\mathrm{B}_{2}|\times\left|\mathsf{f}_{\eta}-\mathsf{f}_{1}\right|_{\infty},$

and thus
$\mathbb{E}[|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\eta}(0)|^{2}]\leq
C\mathbb{E}[|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{1}(0)|^{2}]+C_{\eta}\mathbb{E}[|\mathbf{X}\cap\mathrm{B}_{2}|^{2}]$,
which yields (1.31) using (1.30) (for $\eta=1,r=2$). ∎

#### Renormalized energy

For completeness, we conclude by recalling the definition of the
“renormalized” energy of Sandier-Serfaty in this context. The renormalized
Coulomb energy of a stationary random electric field
$\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ is obtained
by taking the limit $\eta\to 0$ of
$\mathbb{E}\left[|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\eta}(0)|^{2}\right]$
in a renormalized fashion, namely by substracting a divergent part of order
$\mathsf{g}(\eta)$:

$\mathsf{Coul}(\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}):=\lim_{\eta\to
0}\left(\mathbb{E}\left[|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\eta}(0)|^{2}\right]-c_{\mathsf{d}}\mathsf{g}(\eta)\right),$

and the Coulomb energy of a point process is then again phrased as an infimum:

$\mathsf{Coul}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}):=\inf_{\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}\text{
stat., comp. with
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$}}\mathsf{Coul}(\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}).$

The fact that this is the “correct” notion of an infinite-volume Coulomb
energy is supported by a result of $\Gamma$-convergence connecting the finite-
volume energy $\mathrm{H}_{N}$ (as defined in (1.16)) to $\mathsf{Coul}$. We
refer e.g. to [LS17, Sec. 2.7 & 3.] for more details.

#### “Intrinsic” Coulomb energy

One frustrating aspect of the electric approach to defining infinite-volume
Coulomb energies is that it lacks any explicit formula (besides perhaps in the
case of lattices). Given a certain point process, computing - or simply
bounding above - its average Coulomb energy requires to construct a compatible
stationary electric field and then compute its average (renormalized) energy…
which sounds overwhelming.

In [Leb16], inspired by computations done in [BS13], the second author tried
to remedy this issue by introducing a so-called “intrinsic” version of the
Coulomb energy for stationary point processes, denoted by
$\mathbb{W}^{\mathrm{int}}$ in [Leb16] and which we will here denote by
$\mathsf{Coul}^{\mathrm{int}}$. The main advantage of
$\mathsf{Coul}^{\mathrm{int}}$ is that it is _completely explicit and depends
only on the two-point correlation function of the process_. In the “good
cases” (and focusing on dimension $2$):

  1. 1.

The intrinsic energy can be written as:

$\mathsf{Coul}^{\mathrm{int}}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}):=\int_{\mathbb{R}^{2}}-\log|v|(\rho_{2}-1)(v)\mathtt{d}v,$
(1.32)

where $\rho_{2}$ is the two-point correlation function (or measure) of the
process.

  2. 2.

The following “electric-intrinsic” inequality holds, which provides a
condition for finiteness of $\mathsf{Coul}$:

$\mathsf{Coul}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})\leq\mathsf{Coul}^{\mathrm{int}}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}).$
(1.33)

There are, however, two drawbacks in that regard:

  1. 1.

The nice expression (1.32) for $\mathsf{Coul}^{\mathrm{int}}$ is only valid
when $\rho_{2}-1$ has enough integrability at infinity.

  2. 2.

In general, in dimension $2$, the “electric-intrinsic” inequality (1.33) has
to be corrected and written as:

$\mathsf{Coul}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})\leq\mathsf{Coul}^{\mathrm{int}}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})+\mathcal{D}_{\mathrm{log}}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}),$
(1.34)

where $\mathcal{D}_{\mathrm{log}}$ was defined in [Leb16, Sec. 5] as:

$\mathcal{D}_{\mathrm{log}}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}):=C\limsup_{r\to\infty}\left(\sigma(r)\log
r\right)$ (1.35)

where $\sigma$ is the rescaled number variance as in (1.1) and $C$ is some
constant whose value is irrelevant.

It is clear that our condition $\texttt{HU}_{\star}$ (see (1.6)) is related to
the condition $\mathcal{D}_{\mathrm{log}}=0$. We will also explain below (see
Proposition 1.7) how the “intrinsic” approach to Coulomb energy connects (some
of) our results with (some of) the results of [LRY24].

###  1.5 Connection with the literature

#### Finite Wasserstein distance to Lebesgue implies hyperuniformity

Since finite Wasserstein distance to the Lebesgue measure or to a stationary
lattice are equivalent, let us adopt the latter point of view. It has been
known for a long time that:

  1. 1.

Stationary lattices are hyperuniform of class-I i.e.
$\sigma(r)=\mathcal{O}(r^{-1})$.

  2. 2.

Stationary lattices with i.i.d. perturbations admitting a _finite first
moment_ remain hyperuniform of class-I, as proven in [GS75].

The question of whether _dependent_ perturbations could break hyperuniformity
has been raised and studied in the statistical physics literature [Gab04],
with arguments and computations hinting at the fact that - in dimension $2$ \-
applying dependent perturbations with _finite second moment_ to a stationary
lattice should always preserve hyperuniformity. In the recent paper [DFHL24]
this assertion is proven, and examples are given showing that the moment
condition is sharp. In the language of Wasserstein distances, [DFHL24] shows
that:

  1. 1.

If
$\mathsf{Wass}_{2}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}},\mathsf{P}^{\mathbb{Z}^{\mathsf{d}}})<+\infty$
then $\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ is
hyperuniform i.e. the rescaled number variance $\sigma(r)\to 0$ as
$r\to\infty$.

  2. 2.

For all $p<2$, there exists
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ such that
$\mathsf{Wass}_{p}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}},\mathsf{P}^{\mathbb{Z}^{\mathsf{d}}})<+\infty$
and yet (!) $\sigma(r)=+\infty$ for all $r>0$.

It is remarkable that the dimension $2$ plays here a crucial role: for
$\mathsf{d}\geq 3$ it is possible to break hyperuniformity with arbitrarily
small dependent perturbations. Note also that for $\mathsf{d}=1$, the
situation is in a sense better, as finiteness of the $1$-Wasserstein distance
to the stationary lattice is enough to imply hyperuniformity.

#### A partial converse: the result of [LRY24]

In the recent paper [LRY24], the following result is proven (in dimension
$2$):

######  Proposition 1.6 (Thm. 2, (ii) in [LRY24]).

Let $\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ be a
stationary point process such that:

  1. 1.

$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ is
hyperuniform.

  2. 2.

The “correlation measure” $\mathcal{C}$ of
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ is well-
defined as a signed measure, has finite mass, and satisfies:

$\int_{v\in\mathbb{R}^{2},|v|\geq
1}\log|v||\mathcal{C}|(\mathtt{d}v)<+\infty.$ (1.36)

Then the process
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ can be
written as a stationary perturbation, with finite second moment, of a
stationary lattice (and thus $\texttt{Wass}_{2}$ is satisfied).

The “correlation measure” mentioned here is given by
$\mathcal{C}:=(\rho_{2}-\mathrm{Leb})+\delta_{0}$, where $\rho_{2}$ is the
“two-point (or pair) correlation measure” as above, hence the condition
(1.36), which is presented in [LRY24] as _“a mild logarithmic integrability
condition on the reduced pair correlation measure”_ , can be rephrased as:

$\int_{v\in\mathbb{R}^{2},|v|\geq 1}\log|v||\rho_{2}(v)-1|\mathtt{d}v<+\infty$

(at least when $\rho_{2}$ has a density with respect to the Lebesgue measure),
which is very much reminiscent from (1.32) and seems to be controlling some
kind of logarithmic (hence Coulomb) interaction energy of the underlying point
process. Indeed, we find:

######  Proposition 1.7 (Consequences of (1.36)).

Under the assumptions of Proposition 1.6:

  * •

$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ satisfies
$\texttt{HU}_{\star}$, i.e. it has finite regularized Coulomb energy, and thus
(by our Theorem 1) finite $2$-Wasserstein distance to Lebesgue/a lattice.

  * •

Moreover, we have
$\mathcal{D}_{\mathrm{log}}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})=0$
(see (1.35)).

###### Proof of Proposition 1.7.

It is known (see [Cos21, Sec. 1]) that the rescaled number variance in
$\mathrm{B}_{r}$ can be expressed for $r>0$ as:

$\sigma(r)=1+\int_{\mathbb{R}^{2}}j_{r}(v)\mathcal{C}(\mathtt{d}v),$ (1.37)

with
$j_{r}:=\frac{1}{|\mathrm{B}_{r}|}\mathsf{1}_{\mathrm{B}_{r}}\ast\mathsf{1}_{\mathrm{B}_{r}}$,
and that if $\mathcal{C}$ is well-defined as a signed measure, then
hyperuniformity is equivalent (see [Cos21, Sec. 2.2]) to the “sum rule”:

$\int_{\mathbb{R}^{2}}\mathcal{C}(\mathtt{d}v)=-1.$ (1.38)

The function $j_{r}$ is supported on $\mathrm{B}_{2r}$, with $j_{r}(0)=1$, and
is $\mathcal{O}(r^{-1})$-Lipschitz, thus it satisfies:

$j_{r}(v)=1+\mathcal{O}\left(\frac{|v|}{r}\right)\text{ for $|v|\leq
2r$},\quad j_{r}(v)=0\text{ for $|v|\geq 2r$}.$

If $\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ is
hyperuniform, we may use (1.37) and (1.38) to write:

$\sigma(r)=1+\int_{v\in\mathbb{R}^{2},|v|\leq
2r}\mathcal{C}(\mathtt{d}v)+\int_{v\in\mathbb{R}^{2},|v|\leq
2r}\mathcal{O}\left(\frac{|v|}{r}\right)|\mathcal{C}|(\mathtt{d}v)\\\
=-\int_{v\in\mathbb{R}^{2},|v|>2r}\mathcal{C}(\mathtt{d}v)+\int_{v\in\mathbb{R}^{2},|v|\leq
2r}\mathcal{O}\left(\frac{|v|}{r}\right)|\mathcal{C}|(\mathtt{d}v).$ (1.39)

Now, if the condition (1.36) of [LRY24] holds, we can write, on the one hand:

$\left|\int_{v\in\mathbb{R}^{2},|v|>2r}\mathcal{C}(\mathtt{d}v)\right|\leq\frac{1}{\log
r}\int_{v\in\mathbb{R}^{2},|v|>2r}\log|v||\mathcal{C}|(\mathtt{d}v)=o\left(\frac{1}{\log
r}\right)\text{ as }r\to\infty,$

and on the other hand:

$\int_{v\in\mathbb{R}^{2},|v|\leq
2r}\frac{|v|}{r}|\mathcal{C}|(\mathtt{d}v)\leq\frac{1}{\sqrt{r}}\int_{v\in\mathbb{R}^{2},|v|\leq\sqrt{r}}|\mathcal{C}|(\mathtt{d}v)+\frac{2}{\log(r)}\int_{v\in\mathbb{R}^{2},|v|>\sqrt{r}}\log|v||\mathcal{C}|(\mathtt{d}v)\\\
=o\left(\frac{1}{\log r}\right)\text{ as }r\to\infty,$

which gives us $\sigma(r)\log r\to 0$ as $r\to\infty$ and thus indeed
$\mathcal{D}_{\mathrm{log}}=0$ (see (1.35)). This proves the second item.

Moreover, going back to (1.39), taking $r=2^{n}$ and summing over $n\geq 1$,
we find:

$\sum_{n=1}^{+\infty}\sigma(2^{n})\leq\sum_{n=1}^{+\infty}\int_{v\in\mathbb{R}^{2},|v|>2^{n+1}}|\mathcal{C}|(\mathtt{d}v)+C\sum_{n=1}^{+\infty}\int_{v\in\mathbb{R}^{2},|v|\leq
2^{n+1}}\frac{|v|}{2^{n}}|\mathcal{C}|(\mathtt{d}v)\\\
\leq\int_{v\in\mathbb{R}^{2},|v|\geq
1}|\mathcal{C}|(\mathtt{d}v)\left(\sum_{n=1}^{\log|v|}1\right)+C\int_{v\in\mathbb{R}^{2}}|v||\mathcal{C}|(\mathtt{d}v)\left(\sum_{n=\log|v|}^{+\infty}2^{-n}\right)\\\
\leq\int_{v\in\mathbb{R}^{2},|v|\geq
1}\log|v||\mathcal{C}|(\mathtt{d}v)+C\int_{v\in\mathbb{R}^{2}}|\mathcal{C}|(\mathtt{d}v)<+\infty,$

because of the assumptions on $\mathcal{C}$ and the condition (1.36). Thus
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ satisfies
indeed $\texttt{HU}_{\star}$. ∎

######  Remark 1.8.

Since the condition (1.36) (together with hyperuniformity) is found to imply
$\texttt{HU}_{\star}$ and since on the other hand we prove that
$\texttt{HU}_{\star}$ implies Coul, and thus $\texttt{Wass}_{2}$, we do
recover the result of [LRY24, Thm. 2, (ii)]. Of course, here we need to first
pass from $\texttt{HU}_{\star}$ to Coul, and the proof of this implication
uses the analysis of [SWY23] which is a much longer route than the one taken
in [LRY24] (which is, in spirit, very close to the implication
$\texttt{Coul}\implies\texttt{Wass}_{2}$).

Because (as we just proved) the condition of [LRY24] also implies that
$\mathcal{D}_{\mathrm{log}}=0$, it seems to us that another way to deduce Coul
(and thus $\texttt{Wass}_{2}$) from it would be to return to the analysis of
[Leb16] and to find an “intrinsic” formulation for the _regularized_
logarithmic energy (which is in fact easier to handle than the renormalized
one, so the arguments of [Leb16] would readily apply). Without a doubt, since
the regularization corresponds to a short-scale truncation of the interaction,
one would be led to something like:

$\mathsf{Coul}^{\mathrm{int}}_{\eta}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}):=\int_{v\in\mathbb{R}^{2}}-\widetilde{\log}(v)\
(\rho_{2}-1)(v)\mathtt{d}v,$

with $\widetilde{\log}(v)=\log|v|+\mathsf{f}_{\eta}(v)$, which coincides with
$\log|v|$ outside $\mathrm{B}_{\eta}$ and is bounded on $\mathrm{B}_{\eta}$,
and presumably to an “electric-intrinsic” inequality of the form:

$\mathsf{Coul}_{\eta}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})\leq\mathsf{Coul}^{\mathrm{int}}_{\eta}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})+\mathcal{D}_{\mathrm{log}}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}).$

Then hyperuniformity together with (1.36) would imply that
$\mathcal{D}_{\mathrm{log}}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})=0$
and
$\mathsf{Coul}^{\mathrm{int}}_{\eta}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})<+\infty$,
and therefore that Coul, and thus $\texttt{Wass}_{2}$, are satisfied.

#### Wasserstein distance versus discrepancies in finite volume: results and
questions of [BDGZ24]

In the recent preprint [BDGZ24], the authors estimate the Wasserstein distance
in large square boxes between two-dimensional point processes and their “mean”
in terms of controls on the size of “discrepancies” i.e. the difference
between the number of points and its expectation in any square box. Among
their results, the one closest to our concerns is [BDGZ24, Thm.2, case $p=2$],
which we reformulate here in order to make a comparison easier. Let us
consider a stationary point process
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ in
$\mathbb{R}^{2}$, and assume for convenience that
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ has
intensity $1$, i.e.
$\mathbb{E}[|\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}\cap\mathrm{B}_{1}|]=|\mathrm{B}_{1}|$.
If $\Lambda$ is a square, we let
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\Lambda}$
be the restriction of the point process to the square $\Lambda$. For $r>0$ let
$\Lambda_{r}=\left[-\frac{r}{2},\frac{r}{2}\right]^{2}.$

######  Proposition 1.9 ([BDGZ24, Thm.2]).

For $r>0$, let $\bar{\sigma}(r)$ be (a non-increasing upper bound on) the
rescaled number variance of
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ in _the
square_ $\Lambda_{r}$. Then for all $N>1$, we have:

$\frac{1}{|\Lambda_{N}|}\mathbb{E}\left[\mathrm{W}_{2}^{2}\left(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\Lambda_{N}},\frac{|\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}\cap\Lambda_{N}|}{|\Lambda_{N}|}\mathrm{Leb}_{\Lambda_{N}}\right)\right]\leq
C\left(1+\int_{1}^{N}\frac{\sqrt{\bar{\sigma}(s)}}{s}\mathtt{d}s\right).$
(1.40)

Note that the quantity $\bar{\sigma}$ does not control the variance in disks
(as $\sigma$ does), but in squares. In general, the notion of hyperuniformity
can be shape-dependent (see [Cos21, Sec. 2.1] for a discussion), but although
the technique of [BDGZ24] does use the fact that $\Lambda_{N}$ is a square,
this does not seem like a major difference. For non-increasing $\bar{\sigma}$,
it is equivalent to have
$\int_{1}^{+\infty}\frac{\sqrt{\bar{\sigma}(s)}}{s}\mathtt{d}s<+\infty$ and
$\sum_{n\geq 0}\sqrt{\bar{\sigma}(2^{n})}<+\infty$. Hence the result of
[BDGZ24] implies:

$\text{If }\sum_{n\geq 0}\sqrt{\bar{\sigma}(2^{n})}<+\infty,\text{ then
}\limsup_{N\to\infty}\frac{1}{|\Lambda_{N}|}\mathbb{E}\left[\mathrm{W}_{2}^{2}\left(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\Lambda_{N}},\frac{|\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}\cap\Lambda_{N}|}{|\Lambda_{N}|}\mathrm{Leb}_{\Lambda_{N}}\right)\right]<+\infty.$
(1.41)

The assumption on the left-hand side of (1.41) looks very similar to
$\texttt{HU}_{\star}$, with one minor difference (computing the number
variance on squares rather than on disks) and a more significant one which is
the presence of a square root (since $\bar{\sigma}\to 0$, we have
$\sqrt{\bar{\sigma}}\gg\bar{\sigma}$). Of course, if $\bar{\sigma}$ behaves as
a power-law (which is by far the most common case among hyperuniform point
processes) then both $\sum_{n\geq 0}\sqrt{\bar{\sigma}(2^{n})}$ and
$\sum_{n\geq 0}\bar{\sigma}(2^{n})$ will be finite, but strictly speaking the
former assumption is a bit stronger. The conclusion in the right-hand side of
(1.41) is reminiscent of the definitions (1.13), (1.14) given in Section 1.3.
Indeed, arguing as for instance in [HS13, Section 4], see in particular
Corollary 4.5, the right hand side of (1.41) implies finiteness of
$\mathsf{Wass}_{2}^{2}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}},\mathrm{Leb})$.

Let us turn to questions raised in [BDGZ24]. The authors ask: _“Can we find
conditions on two-dimensional point processes for which the Wasserstein
distance is of the optimal order?”_ , to which our Theorems 1 and 4 answer
that (in the case of stationary point processes) $\texttt{HU}_{\star}$, or
equivalently Coul, are sufficient conditions for having $2$-Wasserstein
distances of optimal order.

Moreover, we show in Lemma 5.4 the existence (for any $p\geq 1$) of a
hyperuniform point process
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ such that
$\mathsf{Wass}_{p}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}},\mathrm{Leb})=\infty$,
which provides a negative answer to the possibility raised [BDGZ24, Sec. 4]
that _“it may still be true that a general hyperuniform point process attains
the optimal transport cost”_.

####  $L^{p}$ version of Theorem 1

As explained above, Theorem 1 can be thought of as an infinite volume
equivalent of the classical $\mathrm{W}_{2}-H^{-1}$ inequality. This
inequality has natural $L^{p}$ versions, the $\mathrm{W}_{p}-H^{-1,p}$
inequalities, see [Led19, Theorem 2]. Hence, one could expect that in our
setup
$\mathsf{Wass}_{p}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}},\mathrm{Leb})$
should be upper bounded by
$\mathbb{E}[|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{1}(0)|^{p}]+1$
for any electric field compatible with
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$. However,
there is no obvious interpretation of
$\mathbb{E}[|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{1}(0)|^{p}]$
so the case $p=2$ is the interesting and natural one for us - thus we
formulated our results for $p=2$ only.

It is also worth noting that the Riesz energy (when the Coulomb potential
$\mathsf{g}$ is replaced by a power law $\frac{1}{|\cdot|^{s}}$ for some $s$)
does not seem to have a natural connection to Wasserstein distances.

###  1.6 Notation, plan of the paper

We let $\Lambda_{R}$ be the square $[-\frac{1}{2}R,\frac{1}{2}R]^{\mathtt{d}}$
and write $\mathrm{Leb}_{A}$ for the restriction of the Lebesgue measure to
the set $A$.

We denote deterministic point configurations (resp. fields) with bold face
letters as $\mathbf{X}$ (resp. $\mathrm{E}$) and random point configurations
(resp. fields) with additional bullets as superscripts as
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ (resp.
$\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$).

We sometimes write $A\preceq B$ to express the fact that $A$ is smaller than
$B$ times some universal constant.

#### Plan of the paper.

  * •

In Section 2 we prove that finite regularized Coulomb energy implies finite
Wasserstein distance in any dimension (Theorem 1).

  * •

In Section 3 we prove that finite Wasserstein distance plus a uniform density
bound implies finite regularized Coulomb energy (Theorem 3). We introduce (see
Section 3.2) a new procedure of “screening in Wasserstein space”, strongly
inspired by the “screening” lemma of Sandier-Serfaty, which is a key tool for
the study of Coulomb gases.

  * •

In Section 4, we prove the equivalence between $\texttt{HU}_{\star}$ and Coul
(Theorem 4). To do so, we rely on the results of [SWY23], which we present,
and we show that both $\texttt{HU}_{\star}$ and Coul are equivalent to the
“spectral condition” put forward in [SWY23].

  * •

Finally, in Section 5, we construct simple counter-examples to the
implications $\texttt{Wass}_{2}\implies\texttt{Coul}$ and
$\texttt{HU}\implies\texttt{Wass}_{2}$ (Theorem 2).

##  2 From Coulomb energy to finite Wasserstein distance

###  2.1 The basic $W_{2}$-$H^{-1}$ inequality

The basic tool is the following classical lemma:

######  Lemma 2.1.

Let $\mathbf{m}$ be a probability measure on $\Lambda_{1}$ with bounded
density with respect to the Lebesgue measure and let $v$ be a vector field on
$\Lambda_{1}$ satisfying:

$-\mathrm{div}(v)=\mathbf{m}-\mathrm{Leb}_{\Lambda_{1}},\quad
v\cdot\vec{n}=0\text{ on }\partial\Lambda_{1}.$ (2.1)

Then there exists a coupling $\Pi$ between $\mathbf{m}$ and
$\mathrm{Leb}_{\Lambda_{1}}$ (as probability measures on $\Lambda_{1}$) such
that:

$\iint_{\Lambda_{1}\times\Lambda_{1}}|x-y|^{2}\mathtt{d}\Pi(x,y)\leq
C\int_{\Lambda_{1}}|v|^{2}.$ (2.2)

###### Proof of Lemma 2.1.

The idea is to find a simple interpolation $(\mathbf{m}_{t})_{t\in[0,1]}$
between $\mathrm{Leb}_{\Lambda_{1}}$ and $\mathbf{m}$ and to use $v$ to
construct an associated velocity field $(v_{t})_{t\in[0,1]}$ such that
$(\mathbf{m}_{t},v_{t})_{t\in[0,1]}$ solves the continuity equation (1.9), so
that the associated action gives an upper bound on the $2$-Wasserstein
distance between the measures. A usual choice consists in simply taking a
linear interpolation
$\mathbf{m}_{t}:=(1-t)\mathrm{Leb}_{\Lambda_{1}}+t\mathbf{m}$, however if the
density of $\mathbf{m}$ does not have a positive lower bound then an
integrability issue arises as $t\to 1$. We use instead the choice suggested by
[Led19, Proof of Thm. 2] and set:

$\mathbf{m}_{t}:=\left(1-\theta(t)\right)\mathrm{Leb}_{\Lambda_{1}}+\theta(t)\mathbf{m},\quad\theta(t):=1-(1-t)^{2}.$
(2.3)

It is clear that $\mathbf{m}_{0}=\mathrm{Leb}_{\Lambda_{1}}$,
$\mathbf{m}_{1}=\mathbf{m}$, moreover we can easily compute, using (2.1):

$\frac{\mathtt{d}}{\mathtt{d}t}\mathbf{m}_{t}=\theta^{\prime}(t)(\mathbf{m}-\mathrm{Leb}_{\Lambda_{1}})=-\mathrm{div}\left(\theta^{\prime}(t)v\right)$

and thus, setting $v_{t}:=\frac{\theta^{\prime}(t)v}{\mathrm{m}_{t}}$
($\mathrm{m}_{t}$ denoting the density of $\mathbf{m}_{t}$), we do ensure that
$(\mathbf{m}_{t},v_{t})_{t\in[0,1]}$ solves the continuity equation. Observe
that the density $\mathrm{m}_{t}$ is bounded below by $1-\theta(t)$ for all
$t\in[0,1]$ and compute the action:

$\int_{0}^{1}\int_{\mathbb{R}^{\mathsf{d}}}|v_{t}|^{2}\mathtt{d}\mathbf{m}_{t}\mathtt{d}t\leq\int_{0}^{1}\left(\int_{\Lambda_{1}}\frac{\left(\theta^{\prime}(t)\right)^{2}|v(x)|^{2}}{\mathrm{m}^{2}_{t}(x)}\mathrm{m}_{t}(x)\mathtt{d}x\right)\mathtt{d}t\leq\int_{\Lambda_{1}}|v(x)|^{2}\mathtt{d}x\int_{0}^{1}\frac{\left(\theta^{\prime}(t)\right)^{2}}{1-\theta(t)}\mathtt{d}t,$

and a direct computation shows that the last integral over $\mathtt{d}t$ is
finite. ∎

We deduce the following.

######  Corollary 2.2.

Let $n\geq 1$, let $\mathbf{X}$ be a $n^{\mathsf{d}}$-point configuration in
$\Lambda_{n}$ \- seen here as a purely atomic measure with mass
$n^{\mathsf{d}}$, and let $\eta\in(0,1)$. Assume that
$\min_{x\in\mathbf{X}}\mathrm{dist}(x,\partial\Lambda_{n})>\eta$, and let
$\mathbf{X}_{\eta}$ be the “spread-out” version of $\mathbf{X}$ at distance
$\eta$ as in (1.24). Let $\mathrm{E}$ be an electric field compatible with
$\mathbf{X}$ in $\Lambda_{n}$ and screened in the sense of (1.23), i.e.

$-\mathrm{div}(\mathrm{E})=c_{\mathsf{d}}\left(\mathbf{X}-\mathrm{Leb}\right)\text{
in $\Lambda_{n}$,}\quad\mathrm{E}\cdot\vec{n}=0\text{ on
}\partial\Lambda_{n}.$ (2.4)

Then there exists a coupling $\Pi(\cdot,\cdot;\mathbf{X})$ between
$\mathbf{X}$ and $\mathrm{Leb}_{\Lambda_{n}}$ (as measures of mass
$n^{\mathsf{d}}$ on $\Lambda_{n}$) such that:

$\frac{1}{n^{\mathsf{d}}}\int_{\Lambda_{n}\times\Lambda_{n}}|x-y|^{2}\mathtt{d}\Pi(x,y;\mathbf{X})\leq
C\left(\frac{1}{n^{\mathsf{d}}}\int_{\Lambda_{n}}|\mathrm{E}_{\eta}|^{2}+1\right).$
(2.5)

###### Proof of Corollary 2.2.

Since $\mathrm{E}$ is compatible with $\mathbf{X}$, the regularization
$\mathrm{E}_{\eta}$ is compatible with $\mathbf{X}_{\eta}$ (Claim 1.3) and can
play the role of $v$ in (2.1) (the multiplicative constant $c_{\mathsf{d}}$ is
irrelevant). Thanks to Lemma 2.1 and a simple scaling, this gives an upper
bound on the Wasserstein distance between $\mathbf{X}_{\eta}$ and
$\mathrm{Leb}$. Then the cost of transporting the spread out charges
$\mathbf{X}_{\eta}$ to their centers is of order
$n^{\mathsf{d}}\times\eta^{2}$. ∎

###  2.2 Screening of electric fields in large boxes

In order to apply Corollary 2.2 to stationary point processes, we would like
to be able to work with the restriction of configurations (and fields) to
“large boxes”. Observe that Corollary 2.2 requires a _screened_ electric
field, and (which is in fact a consequence of the mere existence of such a
screened field) that the number of points of $\mathbf{X}$ in $\Lambda_{n}$
agree perfectly with the volume of $\Lambda_{n}$ \- a fact that is absolutely
not guaranteed for random samples.

To fix this, the key tool is an approximation procedure in arbitrarily large
boxes, provided by the “screening” of Sandier-Serfaty, here in the general
version of [PS17, Prop. 6.1] (see also [RS16, Prop. 5.6])

######  Lemma 2.3.

There exists $R_{0}$, $\eta_{0}>0$, and $C$ (all depending only on
$\mathsf{d}$) such that the following holds.

Let $\varepsilon\in(0,\frac{1}{2})$ and let $R$ be such that
$R\geq\frac{R_{0}}{\varepsilon^{2}}$. Assume for simplicity that
$|\Lambda_{R}|$ is an integer. Let $\mathbf{X}$ be a point configuration in
$\Lambda_{R}$ and let $\mathrm{E}$ be an electric field compatible with
$\mathbf{X}$ in $\Lambda_{R}$. Under the following _screenability condition_ :

$\int_{\Lambda_{R}}|\mathrm{E}_{\eta_{0}}|^{2}\leq
R^{\mathsf{d}}\times\left(\varepsilon^{2\mathsf{d}+3}\frac{R}{R_{0}}\right),$
(2.6)

there exists a point configuration $\widetilde{\mathbf{X}}$ in $\Lambda_{R}$
and a local electric field $\widetilde{\mathrm{E}}$ such that:

  1. 1.

$\widetilde{\mathbf{X}}$ has exactly $|\Lambda_{R}|$ points in $\Lambda_{R}$.

  2. 2.

$\widetilde{\mathbf{X}}$ and $\mathbf{X}$ coincide on
$\Lambda_{R(1-\varepsilon)}$.

  3. 3.

All the points of $\widetilde{\mathbf{X}}$ are at distance at least $\eta_{0}$
from $\partial\Lambda_{R}$.

  4. 4.

$\widetilde{\mathrm{E}}$ is compatible with $\widetilde{\mathbf{X}}$ in
$\Lambda_{n}$ and is _screened_.

  5. 5.

We have:

$\int_{\Lambda_{R}}|\widetilde{\mathrm{E}}_{\eta_{0}}|^{2}\leq
C\int_{\Lambda_{R}}|\mathrm{E}_{\eta_{0}}|^{2}.$ (2.7)

######  Remark 2.4.

The screening condition as written in (2.6) is not optimal because we have
condensed two conditions (see [PS17, (6.2)]) into a single one. However, this
simplifies the exposition and is enough for our purposes. Also, the constant
$C$ in (2.7) is in fact of the form $1+C^{\prime}\varepsilon$, but we are not
trying here to derive optimal multiplicative constants for our inequalities.

Let us also record the following fact, which is not part of the original
screening statement but which follows directly from items 3., 4., 5. of Lemma
2.3 in view of Corollary 2.2.

######  Corollary 2.5.

With the assumptions and the notation of Lemma 2.3, we can ensure that there
exists a coupling $\Pi(\cdot,\cdot;\widetilde{\mathbf{X}})$ between
$\widetilde{\mathbf{X}}$ and $\mathrm{Leb}_{\Lambda_{R}}$ such that:

$\frac{1}{R^{\mathsf{d}}}\int_{\Lambda_{R}\times\Lambda_{R}}|x-y|^{2}\mathtt{d}\Pi(x,y;\widetilde{\mathbf{X}})\leq
C\left(\frac{1}{R^{\mathsf{d}}}\int_{\Lambda_{R}}|\mathrm{E}_{\eta_{0}}|^{2}+1\right).$
(2.8)

###  2.3 An approximate sequence of point processes

We may reduce the proof of Theorem 1 to the following proposition.

######  Proposition 2.6.

Let $\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ be a
stationary point process such that
$\mathsf{Coul}_{1}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})$
is finite. Then there exists a sequence
$(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{n})_{n\geq
1}$ of stationary point processes such that:

  1. 1.

$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{n}$
converges in distribution to
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ as
$n\to\infty$.

  2. 2.

$\mathsf{Wass}^{2}_{2}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{n},\mathrm{Leb})\leq
C\left(\mathsf{Coul}_{1}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})+1\right)(1+o_{n}(1))$,
for some constant $C$ depending only on the dimension.

Since $\mathsf{Wass}^{2}_{2}$ is lower semi-continuous on the space of
stationary point processes ([EHJM23, Prop. 2.10.]), this implies (1.2) and
proves Theorem 1.2.

###### Proof of Proposition 2.6.

Let $\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ be a
stationary point process.

##### Step 1. Applying the screening.

######  Lemma 2.7.

There exists some constant $c>0$ (depending on the dimension) such that the
following holds. For all $n\geq 1$, there exists an event $\Omega_{n}$ and a
random point configuration $\widetilde{\mathbf{X}}_{n}$, such that:

  1. 1.

$\lim_{n\to\infty}\mathbb{P}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}\in\Omega_{n})=1$

  2. 2.

Almost every configuration
$\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n}$
has exactly $n^{\mathsf{d}}$ points in $\Lambda_{n}$.

  3. 3.

$\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n}\cap\Lambda_{n-n^{1-c}}$
has the same law as
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}\cap\Lambda_{n-n^{1-c}}$
under the conditional probability $\mathbb{P}[\cdot|\Omega_{n}]$

  4. 4.

For almost every configuration
$\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n}$,
there exists a coupling
$\Pi(\cdot,\cdot;\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n})$
between
$\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n}$
and $\mathrm{Leb}_{\Lambda_{n}}$ such that:

$\frac{1}{n^{\mathsf{d}}}\mathbb{E}\left[\iint_{\Lambda_{n}\times\Lambda_{n}}|x-y|^{2}\mathtt{d}\Pi(x,y;\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n})\right]\leq
C\left(\mathsf{Coul}_{1}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})+1\right)(1+o_{n}(1)).$
(2.9)

###### Proof of Lemma 2.7.

Since we assume that
$\mathsf{Coul}_{1}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})$
is finite, there exists a stationary random field
$\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ compatible
with $\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ such
that:

$\mathsf{Coul}_{1}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})=\mathbb{E}[|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{1}(0)|^{2}]<+\infty.$

Choose $c:=\frac{1}{2\mathsf{d}+4}>0$ and for $n\geq 1$, let
$\varepsilon_{n}:=n^{-c}$. Let $\Omega_{n}$ be the event corresponding to the
screening condition (2.6) with $R=n$ and $\varepsilon=\varepsilon_{n}$ (note
that $n\geq\frac{R_{0}}{\varepsilon_{n}^{2}}$ for $n$ large enough so the
first condition of the screening lemma is satisfied). Then Markov’s inequality
gives:

$\mathbb{P}\left(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}\notin\Omega_{n}\right)\leq\frac{\mathbb{E}\left[\int_{\Lambda_{n}}|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\eta_{0}}|^{2}\right]}{n^{\mathsf{d}}\times\varepsilon_{n}^{2\mathsf{d}+3}\frac{n}{R_{0}}}=\frac{\mathbb{E}\left[\int_{\Lambda_{n}}|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\eta_{0}}|^{2}\right]}{n^{\mathsf{d}}}\times
o_{n}(1),$

because $\varepsilon_{n}$ has been chosen such that
$\varepsilon_{n}^{2\mathsf{d}+3}n\to+\infty$. Moreover, since
$\mathsf{Coul}_{1}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})$
is finite, we know that:

$\limsup_{n\to\infty}\frac{1}{n^{\mathsf{d}}}\mathbb{E}\left[\int_{\Lambda_{n}}|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\eta_{0}}|^{2}\right]<+\infty,$

(using stationarity) and we thus obtain
$\lim_{n\to+\infty}\mathbb{P}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}\notin\Omega_{n})=0$,
which proves the first item.

If $\mathbf{X}\in\Omega_{n}$, we can apply the screening construction of Lemma
2.3 in $\Lambda_{n}$, and we denote by $\widetilde{\mathbf{X}}$ the resulting
point configuration. Let
$\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n}$
be a random variable distributed as
$\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}$
under $\mathbb{P}[\cdot|\Omega_{n}]$.

By design, the screened configurations $\widetilde{\mathbf{X}}$ have exactly
$n^{\mathsf{d}}$ points in $\Lambda_{n}$ and coincide with the original
configuration in $\Lambda_{n-\varepsilon_{n}n}$. This proves the second and
third items.

Finally, using (2.8) and taking the expectation we get the existence of a
random coupling
$\Pi^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}(\cdot,\cdot;\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n})$
coupling
$\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n}$
and $\mathrm{Leb}_{\Lambda_{n}}$, such that:

$\frac{1}{n^{\mathsf{d}}}\mathbb{E}\left[\iint_{\Lambda_{n}\times\Lambda_{n}}|x-y|^{2}\mathtt{d}\Pi^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}(x,y;\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n})\right]\leq
C\left(\frac{1}{n^{\mathsf{d}}}\mathbb{E}\left[\int_{\Lambda_{n}}|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\eta_{0}}|^{2}\big{|}\Omega_{n}\right]+1\right)\\\
=C\left(\frac{1}{n^{\mathsf{d}}\mathbb{P}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}\in\Omega_{n})}\mathbb{E}\left[\int_{\Lambda_{n}}|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\eta_{0}}|^{2}\mathsf{1}_{\Omega_{n}}\right]+1\right)\leq
C\left(\frac{1}{n^{\mathsf{d}}}\mathbb{E}\left[\int_{\Lambda_{n}}|\mathrm{E}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}_{\eta_{0}}|^{2}\right]\left(1+o_{n}(1)\right)+1\right),$

which proves (2.9). ∎

##### Step 2. Forming a stationary approximation.

For $n\geq 1$, let
$\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n}$
be given by Lemma 2.7. We use it to form a stationary point process using a
standard procedure (see e.g. [Leb16, Sec. 5] for an extensive use in the same
context) as follows:

  1. 1.

We take a family of i.i.d. random variables
$\\{\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},(k)}_{n}\\}_{k\in
n\mathbb{Z}^{\mathsf{d}}}$ distributed as
$\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n}$.

  2. 2.

We pave $\mathbb{R}^{\mathsf{d}}$ by copies of $\Lambda_{n}$ centered on the
lattice points $k\in n\mathbb{Z}^{\mathsf{d}}$. We also let
$\tau^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n}$ be some random
variable, independent from the previous ones and uniformly distributed in
$\Lambda_{n}$.

  3. 3.

We define a global random point configuration
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n}$
by:

$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n}:=\bigcup_{k\in
n\mathbb{Z}^{\mathsf{d}}}\bigcup_{x\in\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},(k)}_{n}}\\{x+k\\},$
(2.10)

which amounts to pasting the random point configuration
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},(k)}_{n}$
onto the translated square $k+\Lambda_{n}$.

  4. 4.

We define
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{av}}_{n}$
as
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n}+\tau^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n}$,
which amounts to shifting
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}$
by an independent random vector uniformly distributed in $\Lambda_{n}$.

##### Step 3. Conclusion.

It remains to check that
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n}$
satisfies the properties stated in Proposition 2.6. By construction, it is
stationary. Moreover, by the first and third item of Lemma 2.7, we can
guarantee that it converges in distribution to
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}}$ as
$n\to\infty$.

For $k\in n\mathbb{Z}^{\mathsf{d}}$, let
$\Pi(\cdot,\cdot;\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},(k)}_{n})$
be the (random) coupling between
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},(k)}_{n}$
(seen as a (random) measure on $\Lambda_{n}$) and $\mathrm{Leb}_{\Lambda_{n}}$
given by Lemma 2.7. We can easily combine those couplings to form a (random)
coupling
$\Pi(\cdot,\cdot;\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n})$
between
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n}$
(seen as a (random) measure on $\mathbb{R}^{\mathsf{d}}$) and $\mathrm{Leb}$.
This is due to the fact that two couplings between two measures supported on
two squares that have a side in common will not interfere with each other -
this is somehow equivalent to “screened electric fields” which can be glued
together along a common side.

Finally, if $\tau\in\Lambda_{n}$ we can use
$\Pi(\cdot,\cdot;\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n})$
to define a coupling
$\Pi_{\tau}(\cdot,\cdot;\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n})$
between
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n}+\tau$
and $\mathrm{Leb}$, through a push forward
$\Pi(\cdot,\cdot;\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n})$
by the map $(x,y)\mapsto(x+\tau,y+\tau)$.

######  Claim 2.8.

For all $t\in\Lambda_{n}$, the following convergence holds almost surely and
in $L^{1}$:

$\lim_{R\to+\infty}\frac{1}{|\Lambda_{R}|}\iint_{\Lambda_{R}\times\mathbb{R}^{\mathsf{d}}}|x-y|^{2}\mathtt{d}\Pi_{t}(x,y;\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n})=\frac{1}{n^{\mathsf{d}}}\mathbb{E}\left[\iint_{\Lambda_{n}\times\Lambda_{n}}|x-y|^{2}\mathtt{d}\Pi(x,y;\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n})\right],$
(2.11)

where
$\Pi(\cdot,\cdot;\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n})$
is as given in Lemma 2.3.

###### Proof of Claim 2.8.

Let $\tau\in\Lambda_{n}$ be fixed and let $R\geq 1$. There exists two families
$\mathcal{R},\mathcal{R}^{\prime}$ of indices in $n\mathbb{Z}^{\mathsf{d}}$,
with $\mathcal{R}\subset\mathcal{R}^{\prime}$, both having cardinality
$\frac{R^{\mathsf{d}}}{n^{\mathsf{d}}}+\mathcal{O}(R^{\mathsf{d}-1})$ (with an
implicit constant that may depend on $n$), such that:

$\bigcup_{k\in\mathcal{R}}\left(\Lambda_{n}+k+t\right)\subset\Lambda_{R}\subset\bigcup_{k\in\mathcal{R}^{\prime}}\left(\Lambda_{n}+k+t\right),$
(2.12)

and these inclusions readily imply that:

$\sum_{k\in\mathcal{R}}\iint_{(\Lambda_{n}+k+t)\times\mathbb{R}^{\mathsf{d}}}|x-y|^{2}\mathtt{d}\Pi_{t}(x,y;\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n})\leq\iint_{\Lambda_{R}\times\mathbb{R}^{\mathsf{d}}}|x-y|^{2}\mathtt{d}\Pi_{t}(x,y;\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n})\\\
\leq\sum_{k\in\mathcal{R}^{\prime}}\iint_{(\Lambda_{n}+k+t)\times\mathbb{R}^{\mathsf{d}}}|x-y|^{2}\mathtt{d}\Pi_{t}(x,y;\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n}).$

By construction of the coupling
$\Pi_{t}(\cdot,\cdot;\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n})$
as a gluing of couplings in each translated copy of $\Lambda_{n}$, this
amounts to:

$\sum_{k\in\mathcal{R}}\iint_{\Lambda_{n}\times\Lambda_{n}}|x-y|^{2}\mathtt{d}\Pi(x,y;\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},(k)}_{n})\leq\iint_{\Lambda_{R}\times\mathbb{R}^{\mathsf{d}}}|x-y|^{2}\mathtt{d}\Pi_{t}(x,y;\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{glob}}_{n})\\\
\leq\sum_{k\in\mathcal{R}^{\prime}}\iint_{\Lambda_{n}\times\Lambda_{n}}|x-y|^{2}\mathtt{d}\Pi(x,y;\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},(k)}_{n}),$

and since the
$\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},(k)}_{n}$ are
i.i.d. of law
$\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n}$,
the law of large numbers implies that both:

$\frac{1}{|\mathcal{R}|}\sum_{k\in\mathcal{R}}\iint_{\Lambda_{n}\times\Lambda_{n}}|x-y|^{2}\mathtt{d}\Pi(x,y;\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},(k)}_{n}),\quad\frac{1}{|\mathcal{R}^{\prime}|}\sum_{k\in\mathcal{R}^{\prime}}\iint_{\Lambda_{n}\times\Lambda_{n}}|x-y|^{2}\mathtt{d}\Pi(x,y;\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},(k)}_{n})$

converge a.s. and in $L^{1}$ to
$\mathbb{E}\left[\iint_{\Lambda_{n}\times\Lambda_{n}}|x-y|^{2}\mathtt{d}\Pi(x,y;\widetilde{\mathbf{X}}^{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}_{n})\right]$
as $R\to\infty$. Since both families have cardinality equal to
$\frac{R^{\mathsf{d}}}{n^{\mathsf{d}}}+\mathcal{O}(R^{\mathsf{d}-1})$, this
proves the claim. ∎

Using (2.9) to estimate the right-hand side of (2.11) and taking the
expectation, we find:

$\mathbb{E}\left[\inf_{\Pi\in\mathsf{cpl}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}},\mathrm{av}}_{n},\mathrm{Leb})}\limsup_{R\to+\infty}\frac{1}{|\Lambda_{R}|}\iint_{\Lambda_{R}\times\mathbb{R}^{\mathsf{d}}}|x-y|^{2}\mathtt{d}\Pi(x,y)\right]\leq
C\left(\mathsf{Coul}_{1}(\mathbf{X}^{{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\bullet$}}}}}})+1\right)(1+o_{n}(1)),$
(2.13)

which, in view of the definitions (1.13) and (1.15) concludes the proof of
Proposition 2.6 and thus of Theorem 1. ∎

##  3 $\texttt{Wass}_{2}$ plus density bound implies Coul

In the previous section, we have proven
$\texttt{Coul}\implies\texttt{Wass}_{2}$. The proof of the converse
implication, under the (quite strong) bounded density assumption (1.4),
follows a symmetrical route and uses two ingredients:

  1. 1.

A reverse inequality bounding the $H^{-1}$ norm of the difference of two
measures in terms of their $2$-Wasserstein distance (cf. Lemma 3.1). This can
be found in [Loe06, Prop. 3.1] and we simply reformulate it here in a form
suited to our purposes.

  2. 2.

A new “screening lemma” (Lemma 3.3) in $2$-Wasserstein space, allowing the
construction of approximate point processes with good properties. This
combines the techniques of Sandier-Serfaty with the Benamou-Brenier point of
view on optimal transport.

For simplicity when proving Lemma 3.3, and because finiteness of
$2$-Wasserstein distance to Lebesgue is a less interesting property in
dimension $\mathsf{d}\geq 3$ (even the Poisson point process satisfies it), we
restrict ourselves to $\mathsf{d}=2$, but there is no obstruction to
proceeding similarly in any dimension.

###  3.1 A $H^{-1}$-$W_{2}$ inequality for measures with bounded densities

######  Lemma 3.1.

Let $\mathbf{m}$ be a probability measure on $\Lambda_{1}$ and let
$\bar{\rho}\geq 1$. Assume that $\mathbf{m}$ has a density $\mathrm{m}$ with
respect to the Lebesgue measure, with
$\|\mathrm{m}\|_{\infty}\leq\bar{\rho}<+\infty$. Then, there exists a vector
field $v$ satisfying:

$-\mathrm{div}(v)=\left(\mathbf{m}-\mathrm{Leb}\right)\text{ in
}\Lambda_{1},\quad v\cdot\vec{n}=0\text{ on }\partial\Lambda_{1},$

and such that:

$\int_{\Lambda_{1}}|v|^{2}\leq\bar{\rho}\times\mathrm{W}_{2}^{2}\left(\mathbf{m},\mathrm{Leb}_{\Lambda_{1}}\right).$
(3.1)

This is basically [Loe06, Prop. 3.1] and we sketch here the proof.

###### Proof.

Let $(\mathbf{m}_{t},v_{t})_{t\in[0,1]}$ be an optimal admissible pair in the
Benamou-Brenier formula (1.12) interpolating between
$\mathbf{m}_{0}:=\mathrm{Leb}$ and $\mathbf{m}_{1}:=\mathbf{m}$, and define
the vector field $v$ as $v:=\int_{0}^{1}\mathrm{m}_{t}v_{t}\mathtt{d}t$.
Integrating the continuity equation (1.9) shows that
$-\mathrm{div}(v)=\mathbf{m}_{1}-\mathbf{m}_{0}=\mathbf{m}-\mathrm{Leb}$ in
$\Lambda_{1}$, and since for all $t$ the vector field $v_{t}$ satisfies a
Neumann boundary condition, so does $v$. Moreover, we have by Cauchy-Schwarz:

$\int_{\Lambda_{1}}|v|^{2}\leq\int_{0}^{1}\left(\int_{\Lambda_{1}}\mathrm{m}^{2}_{t}(x)|v_{t}|^{2}(x)\mathtt{d}x\right)\mathtt{d}t\leq\sup_{t\in[0,1]}\|\mathrm{m}_{t}\|_{\infty}\mathrm{W}_{2}^{2}(\mathrm{Leb}_{\Lambda_{1}},\mathbf{m}).$

Using the uniform control (1.10) on the density of $\mathrm{m}_{t}$ concludes
the proof. ∎

######  Corollary 3.2.

Let $n\geq 1$, let $\mathbf{X}$ be a $n^{\mathsf{d}}$-point configuration in
$\Lambda_{n}$ and let $\eta\in(0,1)$. Assume that
$\eta<\min_{x\in\mathbf{X}}\mathrm{dist}(x,\partial\Lambda_{n})$ and that the
following “density bound” holds:

$\sup_{x\in\Lambda_{n}}\left|\mathbf{X}\cap\mathrm{B}_{1}(x)\right|\leq\bar{\rho}.$
(3.2)

Then there exists a screened electric field $\mathrm{E}$ compatible with
$\mathbf{X}$ in $\Lambda_{n}$ such that:

$\frac{1}{n^{\mathsf{d}}}\int_{\Lambda_{n}}|\mathrm{E}_{\eta}|^{2}\leq
C_{\eta}\bar{\rho}\left(\frac{1}{n^{\mathsf{d}}}\mathrm{W}_{2}^{2}(\mathbf{X},\mathrm{Leb}_{\Lambda_{n}})+1\right),$
(3.3)

with $C_{\eta}=C\|\chi\|_{\infty}\eta^{-\mathsf{d}}$ for some $C$ depending
only on $\mathsf{d}$.

###### Proof of Corollary 3.2.

Let $\mathbf{X}_{\eta}$ be the “spread-out” version of $\mathbf{X}$ at
distance $\eta$ as in (1.24). The $2$-Wasserstein distance between the
measures $\frac{1}{n^{\mathsf{d}}}\mathbf{X}$ and
$\frac{1}{n^{\mathsf{d}}}\mathbf{X}_{\eta}$ is clearly bounded by $\eta$, thus
by the triangular inequality we have:

$\mathrm{W}_{2}^{2}(\mathbf{X}_{\eta},\mathrm{Leb}_{\Lambda_{n}})\leq
C\left(\mathrm{W}_{2}^{2}(\mathbf{X},\mathrm{Leb}_{\Lambda_{n}})+n^{\mathsf{d}}\right).$
(3.4)

Regularizing $\mathbf{X}$ into $\mathbf{X}_{\eta}$ produces a measure with
bounded density with respect to the Lebesgue measure. Indeed, combining (3.2)
and the fact that spread-out Dirac masses $\delta^{(\eta)}_{x}$ have a bounded
density, we obtain:

$\|\frac{\mathtt{d}\mathbf{X}_{\eta}}{\mathtt{d}\mathrm{Leb}_{\Lambda_{n}}}\|_{\infty}\leq
C_{\eta}\bar{\rho}.$

We may now apply Lemma 3.1 to
$\frac{1}{n^{\mathsf{d}}}\widetilde{\mathbf{X}_{\eta}}$, where
$\widetilde{\mathbf{X}_{\eta}}$ is the push-forward of $\mathbf{X}_{\eta}$ by
$x\mapsto n^{-1}x$. By construction,
$\frac{1}{n^{\mathsf{d}}}\widetilde{\mathbf{X}_{\eta}}$ is a probability
measure on $\Lambda_{1}$ with density bounded by $C_{\eta}\bar{\rho}$ and such
that:

$\mathrm{W}_{2}^{2}\left(\frac{1}{n^{\mathsf{d}}}\widetilde{\mathbf{X}_{\eta}},\mathrm{Leb}_{\Lambda_{1}}\right)=\frac{1}{n^{\mathsf{d}+2}}\mathrm{W}_{2}^{2}\left(\widetilde{\mathbf{X}_{\eta}},\mathrm{Leb}_{\Lambda_{n}}\right).$

The vector field $v$ given by Lemma 3.1 can then be turned back into a vector
field $\mathrm{E}_{\eta}$ compatible with $\mathbf{X}_{\eta}$ by setting
$\mathrm{E}_{\eta}(x):=nv(n^{-1}x)$ and applying a simple scaling to (3.1)
yields:

$\int_{\Lambda_{n}}|\mathrm{E}_{\eta}|^{2}\leq
C_{\eta}\bar{\rho}\mathrm{W}_{2}^{2}(\mathbf{X}_{\eta},\mathrm{Leb}_{\Lambda_{n}}),$

which combined with (3.4) concludes the proof. ∎

Corollary 3.2 shows that for (finite) point configurations with a well-
controlled density of points, a bound on the Wasserstein distance (to
Lebesgue) can be turned into a bound on the energy of a compatible electric
field. This opens the way to proving an implication of the type
$\texttt{Wass}_{2}\implies\texttt{Coul}$, however our statement requires that
the number of points and the volume of the box coincide, which is not true in
general for a random sample of a large point configuration. We thus develop a
“screening procedure” which should correct the number of points without
changing the underlying configuration - or its distance to Lebesgue - too
much. Note that the construction of [EHJM23, Theorem 3.8] is somewhat similar
in spirit, but it is not suited to our purposes: it ensures that two slightly
modified point processes will have the same random number of points in large
boxes, however, these random numbers are in general different from the
Lebesgue mass of that box.

###  3.2 “Screening” in Wasserstein space

######  Lemma 3.3 (Wasserstein screening).

We can choose $\eta_{0}$ small enough, and $K,C$ large enough such that the
following holds.

Let $\varepsilon\in(0,\frac{1}{2})$ and let $R\geq 10$. Assume for simplicity
that $|\Lambda_{R}|$ is an integer. Let $\mathbf{X}$ be a point configuration
in $\Lambda_{R}$, and let $\Pi$ be a coupling between $\mathbf{X}$ and a
measure $\mathbf{m}_{0}$ on $\mathbb{R}^{2}$ with
$\mathbf{m}_{0}\leq\mathrm{Leb}$, such that the mass of $\mathbf{m}_{0}$ is
equal to $|\mathbf{X}|$. Assume that the points of $\mathbf{X}$ satisfy a
density bound of the form:

$\sup_{x\in\Lambda_{R}}\left|\mathbf{X}\cap\mathrm{B}_{1}(x)\right|\leq\bar{\rho}.$
(3.5)

Assume that the following “screenability conditions” hold:

  1. 1.

The parameters $R,\varepsilon$ and the bound (3.5) are such that
$\bar{\rho}\leq\frac{\varepsilon^{2}R}{K}$ (3.6)

  2. 2.

The number of points of $\mathbf{X}$ in $\Lambda_{R}$ is smaller than
$\frac{1}{K}\frac{\varepsilon^{7}R^{4}}{\bar{\rho}}$ (3.7)

  3. 3.

The quadratic cost of $\Pi$ is controlled by:

$\iint_{\Lambda_{R}\times\mathbb{R}^{2}}|x-y|^{2}\mathtt{d}\Pi(x,y)\leq\frac{1}{K}\frac{\varepsilon^{7}R^{4}}{\bar{\rho}}$
(3.8)

  4. 4.

The restriction of $\mathbf{m}_{0}$ to $\Lambda_{R(1-\varepsilon)}$ coincides
with the Lebesgue measure on $\Lambda_{R(1-\varepsilon)}$. (3.9)

Then there exists a point configuration $\widetilde{\mathbf{X}}$ in
$\Lambda_{R}$ such that:

  1. 1.

$\widetilde{\mathbf{X}}$ has exactly $|\Lambda_{R}|$ points in $\Lambda_{R}$.

  2. 2.

$\widetilde{\mathbf{X}}$ and $\mathbf{X}$ coincide on
$\Lambda_{R(1-\varepsilon)}$.

  3. 3.

All the points of $\widetilde{\mathbf{X}}$ are at distance at least $\eta_{0}$
from $\partial\Lambda_{R}$.

  4. 4.

The following density bound holds:

$\sup_{x\in\Lambda}\left|\widetilde{\mathbf{X}}\cap\mathrm{B}_{1}(x)\right|\leq
C\bar{\rho}.$ (3.10)

  5. 5.

We control the (usual) $2$-Wasserstein distance between
$\widetilde{\mathbf{X}}$ and $\mathrm{Leb}_{\Lambda_{R}}$ by:

$\mathrm{W}_{2}^{2}\left(\widetilde{\mathbf{X}},\mathrm{Leb}_{\Lambda_{R}}\right)\leq
C\left(1+\varepsilon\bar{\rho}\right)\times\iint_{\Lambda_{R}\times\mathbb{R}^{2}}|x-y|^{2}\mathtt{d}\Pi(x,y)+C|\mathbf{X}\cap\Lambda_{R}|+CR^{2}\left(\varepsilon+\varepsilon^{2}\bar{\rho}^{2}\right).$
(3.11)

###### Proof of Lemma 3.3.

We follow the strategy introduced by Sandier-Serfaty for their screening of
electric fields in the context of Coulomb gases, and transpose it into the
$2$-Wasserstein world using the following analogy:

Coulomb |  | Wasserstein
---|---|---
Electric field $\mathrm{E}$ | $\longleftrightarrow$ | Total flux $\int_{0}^{1}\mu_{t}\times\mathrm{v}_{t}\mathtt{d}t$
Electric energy $\int|\mathrm{E}|^{2}$ | $\longleftrightarrow$ | Total cost $\int_{0}^{1}\int\mu_{t}\times|\mathrm{v}_{t}|^{2}\mathtt{d}t$
Compatibility $-\mathrm{div}(\mathrm{E})=c_{\mathsf{d}}\left(\mathbf{X}-\mathrm{Leb}\right)$ | $\longleftrightarrow$ | Continuity equation $-\mathrm{div}(\mu_{t}\mathrm{v}_{t})=\frac{\mathtt{d}}{\mathtt{d}t}\mu_{t}$

One may quickly notice two non-trivial differences: in the Coulomb world, the
energy is the $L^{2}$-norm (squared) of the electric field, whereas on the
Wasserstein side we rather have:

$\text{(Cost) }\int_{0}^{1}\int\mu_{t}|\mathrm{v}_{t}|^{2}\mathtt{d}t\text{
vs.
}\int\left(\int_{0}^{1}\mu_{t}\mathrm{v}_{t}\mathtt{d}t\right)^{2}\leq\int_{0}^{1}\int|\mu_{t}|^{2}|\mathrm{v}_{t}|^{2}\mathtt{d}t\text{
($L^{2}$-norm squared of the total flux)},$

the integrands differing by one power of the density (hence the need to have
density bounds to control the cost in terms of the total flux). Moreover,
screening in the Coulomb world only involves finding _one_ compatible electric
field, whereas the Wasserstein picture, with the Benamou-Brenier point of
view, requires to construct a _family of interpolating measures_ and velocity
fields. Most of the technical ingredients can nevertheless be borrowed from
Coulomb screening as in [PS17, Prop. 6.1].

In short: “screening” consists in making ad hoc constructions near the
boundary of $\Lambda_{R}$ to pass from the existing data to the desired one.
Here we start with $\mathbf{m}_{0},\mathbf{X}$ and an optimal pair
interpolating between $\mathbf{m}_{0}$ and $\mathbf{X}$ (or more precisely a
regularized version thereof) in the sense of Benamou-Brenier - the coupling
$\Pi$ given in the assumption yields an upper bound on the Wasserstein
distance and thus on the action of that pair. We then use this data to produce
a _new_ pair, interpolating this time between $\mathrm{Leb}_{\Lambda_{R}}$ and
(a regularized version of) a new configuration $\widetilde{\mathbf{X}}$. The
new pair coincides with the old one within a large sub-square
$\Lambda_{r}\subset\Lambda_{R}$, but on $\Lambda_{R}\setminus\Lambda_{r}$ we
construct a new configuration (in fact a family of measures), by hand, in
order to end up with exactly the right number of points and to have a proper
interpolation as $t$ varies. This is done by solving some elliptic equations
and defining somewhat abstract “fields”, which will eventually produce the
velocity field needed for the Benamou-Brenier formula. The “action” of the
pair we produce is then evaluated thanks to “energy estimates” found in the
Coulomb screening lemmas.

##### Step 0. Regularizing the configuration and choosing an admissible pair

For $\eta_{0}$ small enough (we choose $\eta_{0}=\frac{1}{100}$ below) we let
$\mathbf{X}_{\eta_{0}}$ be the regularized point configuration as in (1.24).
Of course, there is a natural coupling between $\mathbf{X}$ and
$\mathbf{X}_{\eta_{0}}$ with quadratic cost bounded by
$|\mathbf{X}\cap\Lambda_{R}|\times\eta_{0}^{2}$, so by the triangular
inequality we have:

$\mathrm{W}_{2}^{2}(\mathbf{X}_{\eta_{0}},\mathbf{m}_{0})\preceq\mathrm{W}_{2}^{2}(\mathbf{X},\mathbf{m}_{0})+|\mathbf{X}\cap\Lambda_{R}|)\times\eta_{0}^{2}\preceq\iint_{\Lambda_{R}\times\mathbb{R}^{2}}|x-y|^{2}\mathtt{d}\Pi(x,y)+|\mathbf{X}\cap\Lambda_{R}|\times\eta_{0}^{2},$
(3.12)

the last inequality following from the fact that $\Pi$ is _a_ coupling between
$\mathbf{X}$ and $\mathbf{m}_{0}$, which provides an upper bound on the
Wasserstein distance. The density assumption (3.5) on the number of points of
$\mathbf{X}$ translates into a bound for the density of the measure
$\mathbf{X}_{\eta_{0}}$ with respect to the Lebesgue measure, namely:

$\|\frac{\mathtt{d}\mathbf{X}_{\eta_{0}}}{\mathtt{d}\mathrm{Leb}}\|_{\infty}\leq
C\bar{\rho}\eta_{0}^{-2}\leq C^{\prime}\bar{\rho}.$ (3.13)

We now fix an optimal family $(\mathbf{m}_{t},\mathrm{v}_{t})_{t\in[0,1]}$,
given by the Benamou-Brenier111The Wasserstein distance is a distance, hence,
symmetric. We always write the point configuration (or its spread out version)
first and then the reference measure. However, we will often use the Benamou-
Brenier formulation in “the other direction” i.e. with an interpolation
between the reference measure (at time $0$) and the point configuration (at
time $1$). formula (1.12), such that:

  1. 1.

The initial measure $\mathbf{m}_{0}$ is as in the assumption (in particular
$\mathbf{m}_{0}$ coincides with the Lebesgue measure on
$\Lambda_{R(1-\varepsilon)}$), and $\mathbf{m}_{1}=\mathbf{X}_{\eta_{0}}$.

  2. 2.

The associated cost satisfies:

$\int_{0}^{1}\left(\int_{\mathbb{R}^{2}}|\mathrm{v}_{t}|^{2}\mathtt{d}\mathbf{m}_{t}\right)\mathtt{d}t\leq\mathrm{W}_{2}^{2}(\mathbf{X}_{\eta_{0}},\mathbf{m}_{0}).$
(3.14)

  3. 3.

For all $t\in[0,1]$, $\mathbf{m}_{t}$ has a bounded density $\mathrm{m}_{t}$
with respect to the Lebesgue measure on $\mathbb{R}^{2}$, with:

$\sup_{t\in[0,1]}\|\mathrm{m}_{t}\|_{\infty}\leq C\bar{\rho},$ (3.15)

which is possible in view of (3.13) and (1.10).

##### Step 1. Finding a good boundary.

By a mean value argument, we find $r\in[R-\varepsilon
R,R-\frac{1}{2}\varepsilon R]$ such that:

$\int_{0}^{1}\left(\int_{\partial\Lambda_{r}}|\mathrm{v}_{t}|^{2}\mathrm{m}_{t}\right)\mathtt{d}t\preceq\frac{1}{\varepsilon
R}\int_{0}^{1}\left(\int_{\mathbb{R}^{2}}|\mathrm{v}_{t}|^{2}\mathrm{m}_{t}\right)\mathtt{d}t\preceq\frac{1}{\varepsilon
R}\mathrm{W}_{2}^{2}(\mathbf{X}_{\eta_{0}},\mathbf{m}_{0}),$ (3.16)

using (3.14).

We now define $\mathbf{m}^{\mathrm{int}}$ as the restriction to $\Lambda_{r}$
of the measure $\mathbf{X}_{\eta_{0}}$, which will not be modified. An
important technical remark is that since we have spread out the Dirac masses,
some might intersect $\partial\Lambda_{r}$ and we will need to take care of
them in the construction below.

For $t\in[0,1]$ we define the “integrated flux” at time $t$ as the following
vector field:

$\Phi_{t}:=\int_{0}^{t}\mathrm{m}_{s}\mathrm{v}_{s}\mathtt{d}s,$ (3.17)

and we will be particularly interested in the normal component of $\Phi_{t}$
along $\partial\Lambda_{r}$. Let us note already that:

  * •

(Compatibility.) We have $-\mathrm{div}\Phi_{1}=\mathbf{m}_{1}-\mathbf{m}_{0}$
in $\Lambda_{r}$.

  * •

(Flux along $\partial_{\Lambda_{r}}$.) We have, using an integration by parts
and the continuity equation:

$\int_{\partial\Lambda_{r}}\Phi_{1}\cdot\vec{n}=\int_{\partial\Lambda_{r}}\left(\Phi_{1}-\Phi_{0}\right)\cdot\vec{n}=\int_{0}^{1}\int_{\partial\Lambda_{r}}\frac{\mathtt{d}}{\mathtt{d}t}\Phi_{t}\cdot\vec{n}=\int_{0}^{1}\int_{\partial\Lambda_{r}}\mathrm{m}_{t}\mathrm{v}_{t}\cdot\vec{n}\\\
=\int_{0}^{1}\int_{\Lambda_{r}}-\mathrm{div}\left(\mathrm{m}_{t}\mathrm{v}_{t}\right)=\int_{0}^{1}\int_{\Lambda_{r}}\frac{\mathtt{d}}{\mathtt{d}t}\mathrm{m}_{t}=\mathbf{m}_{1}(\Lambda_{r})-\mathbf{m}_{0}(\Lambda_{r})=|\mathbf{X}_{\eta_{0}}\cap\Lambda_{r}|-|\Lambda_{r}|.$
(3.18)

  * •

(Energy along $\partial_{\Lambda_{r}}$.) We have, using (3.15):

$\int_{\partial\Lambda_{r}}|\Phi_{1}|^{2}\leq\int_{\partial\Lambda_{r}}\int_{0}^{1}\mathrm{m}_{t}^{2}|\mathrm{v}_{t}|^{2}\mathtt{d}s\leq
C\bar{\rho}\int_{0}^{1}\left(\int_{\partial\Lambda_{r}}|\mathrm{v}_{t}|^{2}\mathrm{m}_{t}\right)\mathtt{d}t\\\
\leq\frac{C\bar{\rho}}{\varepsilon
R}\mathrm{W}_{2}^{2}(\mathbf{X}_{\eta_{0}},\mathbf{m}_{0}).$ (3.19)

##### Step 2. Defining the “goal” measure

We now choose $\ell:=\varepsilon^{2}R$, so that $1\ll\ell\ll\varepsilon R$.

###### A first subdivision.

We begin by splitting $\Lambda_{R}\setminus\Lambda_{r}$ into a finite family
$\\{H_{i}\\}_{i\in I}$ of rectangles with sidelengths in
$[\frac{\ell}{2},2\ell]$ such that $m_{i}|H_{i}|\in\mathbb{N}$ for all $i$,
with the “corrected density at time 1” $m_{i}$ defined as:

$m_{i}:=1+\frac{1}{|H_{i}|}\left(\int_{\partial\Lambda_{r}\cap\partial
H_{i}}\Phi_{1}\cdot\vec{n}-n_{i}\right),\quad\text{ where
}n_{i}:=\sum_{x\in\mathbf{X}\cap\Lambda_{r}}\int_{H_{i}}\delta^{(\eta_{0})}_{x}.$
(3.20)

The fact that we can split $\Lambda_{R}\setminus\Lambda_{r}$ into a finite
family $\\{H_{i}\\}_{i\in I}$ of rectangles with sidelength in
$[\frac{\ell}{2},2\ell]$ is fairly elementary. The fact that we can moreover
guarantee $m_{i}|H_{i}|\in\mathbb{N}$ comes from a perturbative argument: we
justify below that $m_{i}$ is of order $1$, hence by “modifying the boundaries
of the $H_{i}$ a little bit”, the quantity $m_{i}|H_{i}|$ quickly reaches an
integer. This is a standard argument in Coulomb screening (see e.g. [PS17,
Proof of Prop. 6.1., Step 2.])

###### Control on $|m_{i}-1|$.

######  Claim 3.4.

Up to choosing $K$ large enough in the screenability condition (3.8), we have
$|m_{i}-1|\leq\frac{1}{100}$, in fact we can ensure that
$|m_{i}-1|\times\delta^{(\eta_{0})}_{0}\leq\frac{1}{100}$, with
$\delta^{(\eta_{0})}_{0}$ as in (1.24).

###### Proof of Claim: 3.4.

First observe that
$|\delta^{(\eta_{0})}_{0}|\leq\frac{1}{\eta_{0}^{2}}\|\chi\|_{\infty}.$ On the
one hand, using Cauchy-Schwarz’s inequality we find:

$\frac{1}{|H_{i}|}\left|\int_{\partial\Lambda_{r}\cap\partial
H_{i}}\Phi_{1}\cdot\vec{n}\right|\leq\frac{1}{\ell^{2}}\left(\int_{\partial\Lambda_{r}\cap\partial
H_{i}}1\right)^{\frac{1}{2}}\left(\int_{\partial\Lambda_{r}}|\Phi_{1}|^{2}\right)^{\frac{1}{2}}\leq\frac{1}{\ell^{2}}\ell^{\frac{1}{2}}\left(\int_{\partial\Lambda_{r}}\int_{0}^{1}|\mathrm{v}_{t}|^{2}\mathrm{m}^{2}_{t}\right)^{\frac{1}{2}},$

using the rough bound $\int_{\partial
H_{i}\cap\partial\Lambda_{r}}|\Phi_{1}|^{2}\leq\int_{\partial\Lambda_{r}}|\Phi_{1}|^{2}$
and the definition of $\Phi_{1}$. Next, using the uniform bound (3.15) on the
density $\mathrm{m}_{t}$ and (3.16), we obtain:

$\left(\int_{\partial\Lambda_{r}}\int_{0}^{1}|\mathrm{v}_{t}|^{2}\mathrm{m}^{2}_{t}\right)^{\frac{1}{2}}\leq\sup_{t\in[0,1]}\|\mathrm{m}_{t}\|^{\frac{1}{2}}_{L^{\infty}}\times\left(\int_{\partial\Lambda_{r}}\int_{0}^{1}|\mathrm{v}_{t}|^{2}\mathrm{m}_{t}\right)^{\frac{1}{2}}\leq
C\bar{\rho}^{\frac{1}{2}}\times\left(\frac{1}{\varepsilon
R}\mathrm{W}_{2}^{2}(\mathbf{X}_{\eta_{0}},\mathbf{m}_{0})\right)^{\frac{1}{2}},$

and thus, using (3.12), we find:

$\frac{1}{|H_{i}|}\left|\int_{\partial\Lambda_{r}\cap\partial
H_{i}}\Phi_{1}\cdot\vec{n}\right|\leq\frac{C}{\ell^{\frac{3}{2}}}\bar{\rho}^{\frac{1}{2}}\times\frac{1}{\varepsilon^{\frac{1}{2}}R^{\frac{1}{2}}}\times\left(\iint_{\Lambda_{R}\times\mathbb{R}^{2}}|x-y|^{2}\mathtt{d}\Pi(x,y)+|\mathbf{X}\cap\Lambda_{R}|\times\eta_{0}^{2}\right)^{\frac{1}{2}}.$

Since $\ell=\varepsilon^{2}R$, we see that the right-hand side can be made
arbitrarily small by choosing the constant $K$ in the screenability conditions
(2) and (3.8) large enough.

Controlling $n_{i}$ (as in (3.20)) is easier: by definition $n_{i}$ is bounded
by the number of points in a $1$-neighborhood of the boundary $\partial H_{i}$
of a rectangle whose sidelenghts are of order $\ell$, and thus $n_{i}\leq
C\bar{\rho}\ell$. In particular we have:

$\frac{n_{i}}{|H_{i}|}\leq
C\frac{\bar{\rho}\ell}{\ell^{2}}=\frac{C\bar{\rho}}{\varepsilon^{2}R},$

which can again be made arbitrarily small by choosing the constant $K$ in the
screenability condition (1) large enough. ∎

######  Remark 3.5.

It is helpful to realize that since $\ell\ll\varepsilon R$, our subdivision of
$\Lambda_{R}\setminus\Lambda_{r}$ into rectangles $\\{H_{i}\\}_{i\in I}$
consists in many ($\approx\frac{1}{\varepsilon}$) concentric “layers” of width
$\ell$. Among those layers, only the first one is actually concerned by the
definition of $n_{i},m_{i}$ as in (3.20), indeed for all rectangles $H_{i}$
besides those in this first layer, we have $n_{i}=0$ and $m_{i}=1$ (since
there is no sidelength of $H_{i}$ that touches $\partial\Lambda_{r}$).

###### A second subdivision.

Next, we subdivide each rectangle $H_{i}$ into a family of smaller rectangles
$\\{\mathcal{R}_{\alpha}\\}_{\alpha\in I_{i}}$ such that:

$\text{ for all $i\in I$, for all $\alpha\in I_{i}$,
}|\mathcal{R}_{\alpha}|=\frac{1}{m_{i}},\text{ the sidelengths of
$\mathcal{R}_{\alpha}$ are between $\frac{1}{3}$ and $3$.}$ (3.21)

This is again a standard construction, see e.g. [PS17, Lemma 6.3]. Since we
ensure that $m_{i}$ is very close to $1$, we can guarantee that the
sidelengths will also be close to $1$ (between $\frac{1}{3}$ and $3$ is a
conservative estimate, but the actual values do not matter anyway).

For each such rectangle, we keep track of the spread-out charges in
$\Lambda_{r}$ that might intersect $\partial\Lambda_{r}$ near
$\mathcal{R}_{\alpha}$ and write:

$n_{\alpha}:=\sum_{x\in\mathbf{X}\cap\Lambda_{r}}\int_{\mathcal{R}_{\alpha}}\delta^{(\eta_{0})}_{x}.$
(3.22)

######  Remark 3.6.

As in Remark 3.5, it is here helpful to keep in mind that if $H_{i}$ is not in
the “first layer” of rectangles neighboring $\partial\Lambda_{r}$, then
$m_{i}=1$ and $n_{\alpha}=0$ for all $\alpha\in I_{i}$. In fact, among each
rectangle $H_{i}$ in the “first layer”, only the $\mathcal{O}(\ell)$
rectangles $\mathcal{R}_{\alpha}$ which belong to “the first layer of the
first layer”, i.e. whose sides intersect $\partial\Lambda_{r}$, are such that
$n_{\alpha}$ may be $\neq 0$.

###### Defining the goal measure.

We define the “goal measure on the exterior” $\mathbf{m}^{\mathrm{ext}}$ as
being obtained by:

  1. 1.

Placing a spread-out Dirac mass at the center $p_{\alpha}$ of each rectangle
$\mathcal{R}_{\alpha}$ (for each $\alpha\in I_{i}$, $i\in I$).

  2. 2.

Completing (in some of the rectangles $\mathcal{R}_{\alpha}$ which touch
$\partial\Lambda_{r}$) the (spread-out) charges that are in $\Lambda_{r}$ and
which may intersect the boundary $\partial\Lambda_{r}$. Those should indeed
not be thrown away when re-defining the measure outside $\Lambda_{r}$.

We choose $\eta_{0}:=\frac{1}{100}$, so that for each $i\in I$ and each
$\alpha\in I_{i}$, we have
$\mathrm{dist}(p_{\alpha},\partial\mathcal{R}_{\alpha})\geq\frac{1}{3}\geq
10\eta_{0}$. We then let $\mathbf{m}^{\mathrm{ext}}$ be as follows:

$\mathbf{m}^{\mathrm{ext}}:=\sum_{i\in I}\sum_{\alpha\in
I_{i}}\delta^{(\eta_{0})}_{p_{\alpha}}+\sum_{i\in I}\sum_{\alpha\in
I_{i}}n_{\alpha},$ (3.23)

keeping in mind (see Remark 3.6) that by construction $n_{\alpha}=0$ besides
the “first layer of first layer” rectangles $\mathcal{R}_{\alpha}$ outside
$\Lambda_{r}$.

We also define the “total goal measure” $\widetilde{\mathbf{X}}_{\eta_{0}}$ as

$\widetilde{\mathbf{X}}_{\eta_{0}}:=\mathbf{m}^{\mathrm{int}}+\mathbf{m}^{\mathrm{ext}}.$
(3.24)

Note that by construction (and in particular because we made sure to complete
the spread-out charges across $\partial\Lambda_{r}$), the measure
$\widetilde{\mathbf{X}}_{\eta_{0}}$ is the regularization at scale $\eta_{0}$
of some underlying point configuration $\widetilde{\mathbf{X}}$, which
justifies our notation.

###### First properties of the total goal measure.

We may already check a certain number of things.

  1. 1.

By construction, the mass of $\mathbf{m}^{\mathrm{ext}}$ is given (using
(3.18) to compute $\int_{\partial\Lambda_{r}}\Phi_{1}\cdot\vec{n}$, and
(3.20)) by:

$|\mathbf{m}^{\mathrm{ext}}|=\sum_{i\in I}m_{i}|H_{i}|+\sum_{i\in
I}\sum_{\alpha\in I_{i}}n_{\alpha}=\sum_{i\in
I}|H_{i}|+\int_{\partial\Lambda_{r}}\Phi_{1}\cdot\vec{n}-\sum_{i\in
I}n_{i}+\sum_{i\in I}\sum_{\alpha\in I_{i}}n_{\alpha}\\\ =\sum_{i\in
I}|H_{i}|+\int_{\partial\Lambda_{r}}\Phi_{1}\cdot\vec{n}-\sum_{i\in
I}n_{i}+\sum_{i\in
I}n_{i}=|\Lambda_{R}\setminus\Lambda_{r}|-\left(|\mathbf{X}_{\eta_{0}}^{\mathrm{int}}\cap\Lambda_{r}|-|\Lambda_{r}|\right),$

and thus the total mass of $\widetilde{\mathbf{X}}_{\eta_{0}}$ is equal to:

$|\widetilde{\mathbf{X}}_{\eta_{0}}|=|\mathbf{m}^{\mathrm{int}}|+|\mathbf{m}^{\mathrm{ext}}|=|\mathbf{X}_{\eta_{0}}^{\mathrm{int}}\cap\Lambda_{r}|+\left(|\Lambda_{R}|-|\Lambda_{r}|-|\mathbf{X}_{\eta_{0}}^{\mathrm{int}}\cap\Lambda_{r}|+|\Lambda_{r}|\right)=|\Lambda_{R}|,$

which means that the underlying point configuration $\widetilde{\mathbf{X}}$
has exactly $|\Lambda_{R}|$ points and satisfies the first item of Lemma 3.3.

  2. 2.

Since we have chosen our “good boundary” $\Lambda_{r}$ with $r\geq
R-\varepsilon R$, and have not touched the configuration inside $\Lambda_{r}$,
the second item of Lemma 3.3 is satisfied.

  3. 3.

By construction, the points $p_{\alpha}$ placed in each $\mathcal{R}_{\alpha}$
are all at distance at least $\frac{1}{3}$ from the boundary, and the other
points are much further away, so the third item is satisfied.

  4. 4.

The density of points within $\Lambda_{r}$ has not been changed, and the
density of points outside $\Lambda_{r}$ is governed by $m_{i}$, which is close
to $1$. Hence the density bound (3.10) holds.

The rest of the proof is devoted to checking that $\widetilde{\mathbf{X}}$
satisfies the last item in our statement, namely the $2$-Wasserstein bound
(3.11). To do so, we will produce a curve
$(\widetilde{\mathbf{m}}_{t},\widetilde{\mathrm{v}}_{t})_{t\in[0,1]}$
satisfying the continuity equation, with a continuous family of measures
$(\widetilde{\mathbf{m}}_{t})_{t\in[0,1]}$ on $\Lambda_{R}$ such that:
$\widetilde{\mathbf{m}}_{0}=\mathrm{Leb}_{\Lambda_{R}}$,
$\widetilde{\mathbf{m}}_{1}=\widetilde{\mathbf{X}}_{\eta_{0}}$ and such that
the associated quadratic cost is controlled. This will give us a bound on the
$2$-Wasserstein distance between $\mathrm{Leb}_{\Lambda_{R}}$ and
$\widetilde{\mathbf{X}}_{\eta_{0}}$ \- which can easily be converted into a a
bound on the distance between $\mathrm{Leb}_{\Lambda_{R}}$ and
$\widetilde{\mathbf{X}}$ by the triangular inequality.

##### Step 3. Defining the goal field

We now define a vector field (which will eventually correspond to the value at
time $1$ of the velocity field $\widetilde{\mathrm{v}}$), meant to be
“compatible” with $\widetilde{\mathbf{X}}_{\eta_{0}}$. We follow the plan used
for Coulomb screening and will import the relevant estimates to control the
“energy” of the field. We start by constructing an “exterior” field
$\mathrm{E}^{\mathrm{ext}}$ corresponding to the “exterior” measure
$\mathbf{m}^{\mathrm{ext}}$ on $\Lambda_{R}\setminus\Lambda_{r}$, and which is
defined as a sum of various fields related to the rectangles $H_{i}$ and
$\mathcal{R}_{\alpha}$ in our subdivisions.

  1. 1.

Let us start with the rectangles $H_{i}$ which are _not_ in the “first layer”,
i.e. whose sides do not intersect $\partial\Lambda_{r}$. These are easy to
treat: since for them the “corrected” density $m_{i}$ is equal to $1$ there is
no need to “correct the flux” (see below) and we simply define, for each
$\alpha\in I_{i}$, a field $\mathrm{E}^{(\alpha,p)}$ corresponding to placing
a charge, namely we solve:

$-\mathrm{div}(\mathrm{E}^{(\alpha,p)})=\left(\delta^{(\eta_{0})}_{p_{\alpha}}-1\right)\text{
in }\mathcal{R}_{\alpha},\quad\mathrm{E}^{(\alpha,p)}\cdot\vec{n}=0\text{ on
}\partial\mathcal{R}_{\alpha},$ (3.25)

which is possible because $|\mathcal{R}_{\alpha}|=\frac{1}{m_{i}}=1$ by
construction. A simple estimate for “placing a charge without flux creation”
(see [PS17, Lemma 6.5]) gives us:

$\int_{\mathcal{R}_{\alpha}}\left|\mathrm{E}^{(\alpha,p)}\right|^{2}\leq C,$
(3.26)

with a constant $C$ depending only on our choice of $\eta_{0}$.

  2. 2.

Let us now focus on rectangles $H_{i}$ which _are_ in the “first layer”.

    1. (a)

Our first task will be to “complete charges” i.e. making sure that the parts
of spread-out charges coming from $\Lambda_{r}$ and intersecting
$\partial\Lambda_{r}$ along $\partial H_{i}$ are not forgotten. For each
$\alpha\in I_{i}$, we solve:

$-\mathrm{div}(\mathrm{E}^{(\alpha,c)})=n_{\alpha}\text{ in
}\mathcal{R}_{\alpha},\quad\mathrm{E}^{(\alpha,c)}\cdot\vec{n}=\frac{-n_{\alpha}}{\int_{\partial\mathcal{R}_{\alpha}\cap\partial\Lambda_{r}}1}\text{
on }\partial\mathcal{R}_{\alpha}$ (3.27)

with $n_{\alpha}$ as in (3.22) (note that we choose the boundary condition
appropriately to make this solvable). The energy estimate for “completing a
charge near the boundary” is found in [PS17, Lemma 6.6] and gives:

$\int_{\mathcal{R}_{\alpha}}\left|\mathrm{E}^{(\alpha,c)}\right|^{2}\leq
Cn_{\alpha}^{2}\leq C\bar{\rho}^{2},$

since $n_{\alpha}$ is controlled by the density bound for the points inside
$\Lambda_{r}$.

    2. (b)

The second task is to solve (3.25) again but with the correct density ($m_{i}$
instead of $1$), namely:

$-\mathrm{div}(\mathrm{E}^{(\alpha,p)})=\left(\delta^{(\eta_{0})}_{p_{\alpha}}-m_{i}\right)\text{
in }\mathcal{R}_{\alpha},\quad\mathrm{E}^{(\alpha,p)}\cdot\vec{n}=0\text{ on
}\partial\mathcal{R}_{\alpha},$ (3.28)

which is possible because $|\mathcal{R}_{\alpha}|=\frac{1}{m_{i}}$. The energy
estimate (3.26) still holds because $m_{i}$ is close to $1$ (see again [PS17,
Lemma 6.5]).

    3. (c)

Finally, we need to “correct the flux” on $H_{i}$, which means here solving:

$-\mathrm{div}(\mathrm{E}^{(i,f)})=\left(m_{i}-1\right)\text{ in
}H_{i},\quad\mathrm{E}^{(i,f)}\cdot\vec{n}=\begin{cases}\Phi_{1}\cdot\vec{n}+\sum_{\alpha\in
I_{i}}\frac{n_{\alpha}}{\int_{\partial\mathcal{R}_{\alpha}\cap\partial\Lambda_{r}}1}\mathsf{1}_{\partial\mathcal{R}_{\alpha}}&\text{
on }\partial H_{i}\cap\partial\Lambda_{r},\\\
\mathrm{E}^{(i,f)}\cdot\vec{n}=0&\text{ on }\partial
H_{i}\setminus\partial\Lambda_{r}\end{cases}.$ (3.29)

Note that this is compatible with our definition (3.20) of $m_{i}$. The energy
of $\mathrm{E}^{(i,f)}$ can be estimated using [PS17, Lemma 6.4], we obtain:

$\int_{H_{i}}\left|\mathrm{E}^{(i,f)}\right|^{2}\leq C\ell\left(\int_{\partial
H_{i}\cap\partial\Lambda_{r}}|\Phi_{1}|^{2}+\sum_{\alpha\in
I_{i}}n^{2}_{\alpha}\right)\leq C\ell\int_{\partial
H_{i}\cap\partial\Lambda_{r}}|\Phi_{1}|^{2}+C\ell^{2}\bar{\rho}^{2},$ (3.30)

using again our density assumption to bound each $n_{\alpha}$ by
$C\bar{\rho}$, and the fact that there are $\mathcal{O}(\ell)$ indices
$\alpha$ in $I_{i}$ such $n_{\alpha}\neq 0$.

  3. 3.

We then define the “exterior” field $\mathrm{E}^{\mathrm{ext}}$ as:

$\mathrm{E}^{\mathrm{ext}}:=\sum_{i\in I}\mathrm{E}^{(i,f)}+\sum_{i\in
I}\sum_{\alpha\in I_{i}}\mathrm{E}^{(\alpha,p)}+\sum_{i\in I}\sum_{\alpha\in
I_{i}}\mathrm{E}^{(\alpha,c)},$ (3.31)

knowing that $\mathrm{E}^{(i,f)}=0$ if $H_{i}$ is not in the “first layer”,
and that $\mathrm{E}^{(\alpha,c)}=0$ if $\alpha$ is not in “the first layer of
the first layer”.

  4. 4.

We finally define the “total goal field” (or field at time $1$)
$\mathrm{E}_{1}$ as:

$\mathrm{E}_{1}:=\begin{cases}\Phi_{1}&\text{ on }{\Lambda_{r}}\\\
\mathrm{E}^{\mathrm{ext}}&\text{ on }\Lambda_{R}\setminus\Lambda_{r}\\\
0&\text{ outside }\Lambda_{R}.\end{cases}$

The field $\mathrm{E}_{1}$ is defined piecewise but we have made it so that
the normal components of each piece agree along the common boundary.

###### Properties of the goal field.

  1. 1.

By construction, we have:

$-\mathrm{div}(\mathrm{E}^{\mathrm{ext}})=\left(\mathbf{m}^{\mathrm{ext}}-\mathrm{Leb}\right)\text{
in $\Lambda_{R}\setminus\Lambda_{r}$},$ (3.32)

with a normal component that vanishes on the outer boundary
$\partial\Lambda_{R}$ and matches the one of $\Phi_{1}$ along the inner
boundary $\partial\Lambda_{r}$. This ensures that the total goal field
$\mathrm{E}_{1}$ satisfies:

$-\mathrm{div}(\mathrm{E}_{1})=\left(\widetilde{\mathbf{X}}_{\eta_{0}}-\mathrm{Leb}\right)\text{
in }\Lambda_{R},\quad\mathrm{E}_{1}\cdot\vec{n}=0\text{ on
}\partial\Lambda_{R}.$ (3.33)

  2. 2.

We can easily estimate the energy
$\int_{\Lambda_{R}\setminus\Lambda_{r}}|\mathrm{E}^{\mathrm{ext}}|^{2}$ of
$\mathrm{E}^{\mathrm{ext}}$ using the controls mentioned above for each field.

    1. (a)

First, we have:

$\int_{\Lambda_{R}\setminus\Lambda_{r}}\left|\sum_{i\in
I}\mathrm{E}^{(i,f)}\right|^{2}=\sum_{i\in
I}\int_{H_{i}}\left|\mathrm{E}^{(i,f)}\right|^{2}=\sum_{i\in I,\partial
H_{i}\cap\partial\Lambda_{r}\neq\emptyset}\int_{H_{i}}\left|\mathrm{E}^{(i,f)}\right|^{2}.$

The first equality is due to the fact that the “flux correction” fields
$\mathrm{E}^{(i,f)}$ are zero outside “their” rectangle $H_{i}$ and thus
mutually orthogonal. The second equality is due to the fact that by
construction, we are “correcting” the density $m_{i}$ only on rectangles that
have a side included in $\partial\Lambda_{r}$ \- on the other ones, the
density is simply $1$.

Using (3.30), noting that there are $\mathcal{O}(\frac{R}{\ell})$ rectangles
$H_{i}$ in the “first layer”, and inserting (3.19), we find:

$\int_{\Lambda_{R}\setminus\Lambda_{r}}\left|\sum_{i\in
I}\mathrm{E}^{(i,f)}\right|^{2}\leq
C\ell\int_{\partial\Lambda_{r}}|\Phi_{1}|^{2}+C\times\frac{R}{\ell}\times\ell^{2}\bar{\rho}^{2}\leq
C\ell\frac{{\bar{\rho}}}{\varepsilon
R}\mathrm{W}_{2}^{2}(\mathbf{X}_{\eta_{0}},\mathbf{m}_{0})+CR\ell\bar{\rho}^{2}.$
(3.34)

    2. (b)

Secondly, we have, for the same orthogonality reasons:

$\int_{\Lambda_{R}\setminus\Lambda_{r}}\left|\sum_{i\in I}\sum_{\alpha\in
I_{i}}\mathrm{E}^{(\alpha,p)}\right|^{2}=\sum_{i\in I}\sum_{\alpha\in
I_{i}}\int_{\mathcal{R}_{\alpha}}\left|\mathrm{E}^{(\alpha,p)}\right|^{2}\leq
C\times|\Lambda_{R}\setminus\Lambda_{r}|\times 1=\mathcal{O}(\varepsilon
R^{2}).$ (3.35)

    3. (c)

Finally, we can evaluate the total energy due to the completion of charges as:

$\int_{\Lambda_{R}\setminus\Lambda_{r}}\left|\sum_{i\in I}\sum_{\alpha\in
I_{i}}\mathrm{E}^{(\alpha,c)}\right|^{2}\leq
C\times\frac{R}{\ell}\times\ell\times\bar{\rho}^{2},$ (3.36)

using the fact that there are $\mathcal{O}(\frac{R}{\ell})$ rectangles $H_{i}$
in the “first layer” and that in each of this rectangle there is at most
$\mathcal{O}(\ell)$ smaller rectangles $\mathcal{R}_{\alpha}$ whose
$n_{\alpha}\neq 0$.

In conclusion, combining (3.34) and (3.35) (the last estimate (3.36) being of
smaller order) we find:

$\int_{\Lambda_{R}\setminus\Lambda_{r}}|\mathrm{E}^{\mathrm{ext}}|^{2}\leq
C\ell\frac{{\bar{\rho}}}{\varepsilon
R}\mathrm{W}_{2}^{2}(\mathbf{X}_{\eta_{0}},\mathbf{m}_{0})+CR\ell\bar{\rho}^{2}+C\varepsilon
R^{2}.$ (3.37)

######  Remark 3.7.

For the sole purpose of Theorem 3, one could in fact stop here. Indeed, we
have constructed a point configuration $\widetilde{\mathbf{X}}$ which
approximates the initial one, has the right number of points, and for which we
can exhibit a compatible electric field ($\mathrm{E}_{1}$) which is screened
and whose energy is controlled in terms of
$\mathrm{W}_{2}^{2}(\mathbf{X}_{\eta_{0}},\mathbf{m}_{0})$. This is enough
material to bound Coulomb energies of stationary point processes in terms of
Wasserstein distances, see the proof of Proposition 3.10 below. However, for
completeness and future reference, let us proceed with the proof of our
“Wasserstein screening” as stated, which now requires that we extend our “goal
field” into a family of measures and fields indexed by the time parameter $t$.

##### Step 4. Constructing an interpolation

For each $i\in I$ we introduce the “density correction” at time $t\in[0,1]$
as:

$\mathrm{dens}_{i}(t):=\frac{1}{|H_{i}|}\int_{\partial
H_{i}\cap\partial\Lambda_{r}}\left(\Phi_{t}-t\Phi_{1}\right)\cdot\vec{n}.$
(3.38)

Again, $\mathrm{dens}_{i}$ is actually $0$ unless $H_{i}$ belongs to the
“first layer” of our subdivision. The quantity $\mathrm{dens}_{i}$ can be
positive or negative, but by the same argument as in Claim 3.4, we find that
by choosing $K$ large enough in (3.8):

$\sup_{t\in[0,1]}\sup_{i\in
I}\frac{1}{m_{i}}\left|\mathrm{dens}_{i}(t)\right|\leq\frac{1}{100}\min\left(1,\frac{\eta_{0}^{2}}{\|\chi\|_{\infty}}\right),$
(3.39)

with $\chi$ the bump function chosen for truncations (this observation will be
useful later). We can already note that for all $i\in I$ we have

$\mathrm{dens}_{i}(0)=\mathrm{dens}_{i}(1)=0.$ (3.40)

Moreover, in view of the definition (3.17) of $\Phi_{t}$:

$\frac{\mathtt{d}}{\mathtt{d}t}\mathrm{dens}_{i}(t)=\frac{1}{|H_{i}|}\int_{\partial
H_{i}\cap\partial\Lambda_{r}}\left(\mathrm{m}_{t}\mathrm{v}_{t}-\Phi_{1}\right)\cdot\vec{n}.$
(3.41)

We now consider the corresponding “field correction”, defined as the solution
to:

$-\mathrm{div}(\mathrm{E}^{\mathrm{corr},(i)}_{t})=\sum_{\alpha\in
I_{i}}\frac{1}{m_{i}}\frac{\mathtt{d}}{\mathtt{d}t}\mathrm{dens}_{i}(t)\times\delta^{(\eta_{0})}_{p_{\alpha}}\text{
in }H_{i},\\\
\mathrm{E}^{\mathrm{corr},(i)}_{t}\cdot\vec{n}=\begin{cases}(\mathrm{m}_{t}\mathrm{v}_{t}-\Phi_{1})\cdot\vec{n}&\text{
on }\partial H_{i}\cap\partial\Lambda_{r},\\\ 0&\text{ on }\partial
H_{i}\setminus\partial\Lambda_{r}.\end{cases}$ (3.42)

Note that by construction there are $|H_{i}|m_{i}$ rectangles
$\mathcal{R}_{\alpha}$ (and thus $|H_{i}|m_{i}$ points $p_{\alpha}$) within
$H_{i}$, hence we are placing a total mass $|H_{i}|\times
m_{i}\times\frac{1}{m_{i}}\times\frac{\mathtt{d}}{\mathtt{d}t}\mathrm{dens}_{i}(t)$,
which, by definition (3.41) of
$\frac{\mathtt{d}}{\mathtt{d}t}\mathrm{dens}_{i}(t)$, is compatible with the
boundary data chosen for $\mathrm{E}^{\mathrm{corr},(i)}_{t}$.

######  Claim 3.8 (Energy of the field correction).

For some $C$:

$\sum_{i\in I}\int_{H_{i}}|\mathrm{E}^{\mathrm{corr},(i)}_{t}|^{2}\leq
C\ell\int_{\partial\Lambda_{r}}\left(\bar{\rho}\mathrm{m}_{t}|\mathrm{v}_{t}|^{2}+|\Phi_{1}|^{2}\right).$
(3.43)

###### Proof of Claim 3.8.

For each $H_{i}$ in the first layer, we can decompose
$\mathrm{E}^{\mathrm{corr},(i)}_{t}$ as a sum of (using the terminology of
[PS17]):

  1. 1.

A field that “corrects the flux” i.e. solves

$-\mathrm{div}(\mathrm{E}^{(i,f)}_{t})=\frac{\mathtt{d}}{\mathtt{d}t}\mathrm{dens}_{i}(t)\text{
in
}H_{i},\quad\mathrm{E}^{(i,f)}_{t}\cdot\vec{n}=\begin{cases}(\mathrm{m}_{t}\mathrm{v}_{t}-\Phi_{1})\cdot\vec{n}&\text{
on }\partial H_{i}\cap\partial\Lambda_{r}\\\ 0&\text{ on }\partial
H_{i}\setminus\partial\Lambda_{r}.\end{cases}$

  2. 2.

Fields that “place a charge without flux creation” in each
$\mathcal{R}_{\alpha}$ ($\alpha\in I_{i}$) i.e. solve

$-\mathrm{div}(\mathrm{E}^{(\alpha,p)}_{t})=\frac{\mathtt{d}}{\mathtt{d}t}\mathrm{dens}_{i}(t)\left(\frac{1}{m_{i}}\delta^{(\eta_{0})}_{p_{\alpha}}-1\right)\text{
in }\mathcal{R}_{\alpha},\quad\mathrm{E}^{(\alpha,p)}_{t}\cdot\vec{n}=0\text{
on }\partial\mathcal{R}_{\alpha},$

which again is solvable because $|\mathcal{R}_{\alpha}|=\frac{1}{m_{i}}$.

Using the corresponding energy estimates from [PS17] (as mentioned above), we
find:

$\int_{H_{i}}|\mathrm{E}^{\mathrm{corr},(i)}_{t}|^{2}\leq C\ell\int_{\partial
H_{i}\cap\partial\Lambda_{r}}\left(|\mathrm{m}_{t}|^{2}|\mathrm{v}_{t}|^{2}+|\Phi_{1}|^{2}\right)+C\ell^{2}\times\left|\frac{\mathtt{d}}{\mathtt{d}t}\mathrm{dens}_{i}(t)\right|^{2},$
(3.44)

which we can translate into (3.43) using the uniform bound (3.15) on the
density $\mathrm{m}_{t}$. Indeed, using (3.41) and Cauchy-Schwarz’s
inequality, we see that:

$\ell^{2}\left|\frac{\mathtt{d}}{\mathtt{d}t}\mathrm{dens}_{i}(t)\right|^{2}\leq
C\frac{\ell^{2}}{\ell^{4}}\times\ell\times\int_{\partial
H_{i}\cap\partial\Lambda_{r}}\left(|\mathrm{m}_{t}|^{2}|\mathrm{v}_{t}|^{2}+|\Phi_{1}|^{2}\right),$

which is negligible compared to the first term in the right-hand side of
(3.44). Summing (3.44) over $i\in H_{i}$ yields (3.43). ∎

###### Interpolation measure and field

For $t\in[0,1]$, we let $\theta(t):=1-(1-t)^{2}$ (cf. the proof of Lemma 2.1)
and define:

  * •

The measure $\mathbf{m}^{\mathrm{ext}}_{t}$ as:

$\mathbf{m}^{\mathrm{ext}}_{t}:=t\mathbf{m}^{\mathrm{ext}}+\left(1-t\right)\mathrm{Leb}_{\Lambda_{R}\setminus\Lambda_{r}}+\sum_{i\in
I}\sum_{\alpha\in
I_{i}}\frac{1}{m_{i}}\mathrm{dens}_{i}(t)\delta^{(\eta_{0})}_{p_{\alpha}},$
(3.45)

and the measure $\widetilde{\mathbf{m}}_{t}$ (total measure at time $t$)

$\widetilde{\mathbf{m}}_{t}:=\mathbf{m}_{\theta(t)}\mathsf{1}_{\Lambda_{r}}+\mathbf{m}^{\mathrm{ext}}_{\theta(t)}\mathsf{1}_{\Lambda_{R}\setminus\Lambda_{r}},$
(3.46)

where we recall that $\mathbf{m}_{t}$ was chosen in Step 0.

  * •

The fields $\mathrm{E}^{\mathrm{ext}}_{t}$ as :

$\mathrm{E}^{\mathrm{ext}}_{t}:=\theta^{\prime}(t)\left(\mathrm{E}^{\mathrm{ext}}+\sum_{i\in
I}\mathrm{E}^{\mathrm{corr},(i)}_{\theta(t)}\right),$ (3.47)

with $\mathrm{E}^{\mathrm{ext}}$ as in (3.31) and
$\mathrm{E}^{\mathrm{corr},(i)}$ as in (3.42).

  * •

The velocity field $\widetilde{\mathrm{v}}_{t}$ as:

$\widetilde{\mathrm{v}}(t):=\begin{cases}\theta^{\prime}(t)\mathrm{v}_{\theta(t)}&\text{
on }\Lambda_{r}\\\
\frac{\mathrm{E}^{\mathrm{ext}}_{t}}{\mathrm{m}^{\mathrm{ext}}_{\theta(t)}}&\text{
on }\Lambda_{R}\setminus\Lambda_{r}.\end{cases}$ (3.48)

where $\mathrm{m}^{\mathrm{ext}}_{t}$ denotes the density of
$\mathbf{m}^{\mathrm{ext}}_{t}$, which was defined in (3.45).