#  Existence and nonexistence of minimizer for Thomas-Fermi-Dirac-von
Weizsäcker model on lattice graph

Yong Liu  Yong Liu, Department of Mathematics, University of Science and
Technology of China, Hefei, China<EMAIL_ADDRESS>,  Jun Wang  Jun Wang,
Institute of Applied System Analysis, Jiangsu University, Zhenjiang, Jiangsu,
212013, P.R. China<EMAIL_ADDRESS>,  Kun Wang  Kun Wang, Institute of
Applied System Analysis, Jiangsu University, Zhenjiang, Jiangsu, 212013, P.R.
China<EMAIL_ADDRESS>and  Wen Yang  Wen Yang, Department of
Mathematics, Faculty of Science and Technology, University of Macau, Macau,
P.R. China<EMAIL_ADDRESS>
###### Abstract.

The focus of our paper is to investigate the possibility of a minimizer for
the Thomas-Fermi-Dirac-von Weizsäcker model on the lattice graph
$\mathbb{Z}^{3}$. The model is described by the following functional:

$E(\varphi)=\sum_{y\in\mathbb{Z}^{3}}\left(|\nabla\varphi(y)|^{2}+(\varphi(y))^{\frac{10}{3}}-(\varphi(y))^{\frac{8}{3}}\right)+\sum_{x,y\in\mathbb{Z}^{3}\atop~{}\
y\neq x\hfill}\frac{{\varphi}^{2}(x){\varphi}^{2}(y)}{|x-y|},$

with the additional constraint that
$\sum\limits_{y\in\mathbb{Z}^{3}}{\varphi}^{2}(y)=m$ is sufficiently small. We
also prove the nonexistence of a minimizer provided the mass $m$ is adequately
large. Furthermore, we extend our analysis to a subset
$\Omega\subset\mathbb{Z}^{3}$ and prove the nonexistence of a minimizer for
the following functional:

$E(\Omega)=|\partial\Omega|+\sum_{x,y\in\Omega\atop~{}y\neq
x\hfill}\frac{1}{|x-y|},$

under the constraint that $|\Omega|=V$ is sufficiently large.

The research of Y. Liu is partially supported by NSFC No. 11971026 and No.
12141105. The research of J. Wang is partially supported by National Key
R$\&$D Program of China (2022YFA1005601), NSFC No. 11971202 and Outstanding
Young foundation of Jiangsu Province No. BK20200042. The research of W. Yang
is partially supported by National Key R&D Program of China 2022YFA1006800,
SRG 2023-00067-FST, NSFC No. 12171456 and NSFC No. 12271369.

Keywords: Variational methods, Lattice graph; Nonlocal term; Minimizer
solution.

##  1\. Introduction

The aim of this paper is to investigate whether a minimizer exists for two
mass-constrained variational problems on a lattice graph. Specifically, we are
interested in the Thomas-Fermi-Dirac-von Weizsäcker (TFDW) functional

$E(\varphi)=\int_{\mathbb{R}^{3}}\left(|\nabla\varphi|^{2}+\varphi^{\frac{10}{3}}-\varphi^{\frac{8}{3}}\right)dx+\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{{\varphi}^{2}(x){\varphi}^{2}(y)}{|x-y|}dxdy$
(1.1)

and the Gamow’s liquid drop functional

$E(\Omega)=|\partial\Omega|+\int_{\Omega}\int_{\Omega}\frac{1}{|x-y|}dxdy,$
(1.2)

where $\Omega\subset\mathbb{R}^{3}$ is a bounded set.

The functional (1.1) is commonly known in the Physics literature as the
Thomas-Fermi-Dirac-von Weizsäcker (TFDW) model for electrons in the absence of
an external potential(the electron density is represented by
$\varphi^{2}=\sigma$). The ionization conjecture in quantum mechanics is also
associated with this model. To introduce it, we define the minimized
variational problem

$I_{Z}^{TFDW}:=\inf\left\\{E^{TFDW}(\varphi):\varphi\in
H^{1}(\mathbb{R}^{3}),\int_{\mathbb{R}^{3}}{\varphi^{2}}dx=N\right\\},$ (1.3)

where

$\begin{split}E^{TFDW}(\varphi)&=C^{W}\int_{\mathbb{R}^{3}}|\nabla\varphi|^{2}dx+C^{TF}\int_{\mathbb{R}^{3}}\varphi^{\frac{10}{3}}-C^{D}\int_{\mathbb{R}^{3}}\varphi^{\frac{8}{3}}dx\\\
&\quad-Z\int_{\mathbb{R}^{3}}\frac{\varphi^{2}}{|x|}dx+\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{{\varphi}^{2}(x){\varphi}^{2}(y)}{|x-y|}dxdy.\end{split}$
(1.4)

The ionization conjecture states that

Ionization conjecture: The number of electrons that can be bound to an atomic
nucleus of charge $Z$ cannot exceed $Z+1$. Mathematically, if $N\leq Z+1$,
then $I_{Z}^{TFDW}$ has a minimizer.

For more details statements about this problem one can refer to [49, Problem
9] and [34, Chapter 12]. Although experimental data shows that a neutral atom
has the capacity to bind up to two extra electrons, it remains challenging to
establish a rigorous justification for this observation based on the
fundamental principles of quantum mechanics. Over the last few years, a
significant number of mathematicians have dedicated their efforts to analyzing
the problem (1.3)-(1.4). In the paper [35], the author proved the existence of
solution to the classical TFDW problem (1.3) for $N\leq Z$ and Le Bris [36]
extended it to $N\leq Z+C$ for some $C>0$. Fefferman and Seco et. al in [10,
50] proved that if $Z\rightarrow\infty$, then $N_{c}(Z)\leq Z+O(Z^{5/7})$. In
[31, 42], they proved that for all $Z\geq 1$,
$N_{c}(Z)<\min(2Z+1,1.22Z+3Z^{1/3})$ holds true. In the paper [52], the author
prove that there exists a universal constant $C>0$ such that if $N>Z+C$, then
$E^{HF}(N)$ has no minimizer, where $0\leq\gamma\leq 1$, $\gamma^{2}=\gamma$
and $Tr\gamma=N$

$E^{HF}(N)=\inf_{Tr\gamma=N}\left(Tr(-\Delta-Z|x|^{-1})\gamma)+\frac{1}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{\gamma(x;x)\gamma(y;y)-\gamma(x;y)}{|x-y|}dxdy\right\\}$

and $\gamma$ is a density matrix with eigenfunctions $u_{j}$ and corresponding
eigenvalues $\nu_{j}$ on either $L^{2}(\mathbb{R}^{3})$ or
$L^{2}(\mathbb{R}^{3};\mathbb{C}^{2})$ we shall write

$Tr[-\Delta\gamma]=\Sigma_{j}\nu_{j}\int_{\mathbb{R}^{3}}|\nabla
u_{j}|^{2}dx.$

In the recent papers [44, 12] proved that there exists a constant $C>0$ such
that for all $Z>0$, if $I_{Z}^{TFDW}$ has a minimizer, then $N\leq Z+C$. For
more results on this direction one can refer to the references [2, 11, 30, 32,
33, 43, 47, 48, 51, 54, 56] and references therein.

To the best of our knowledge the ionization conjecture was still not complete
solved in the continue case. For the special case $Z=0$, the ionization
conjecture can also be interpreted as the vacuum can only bind a finite number
of electrons in the TFDW model. The work [37] established that for a large
mass, the functional $I_{0}^{TFDW}$ given by

$I_{0}^{TFDW}=\inf_{\varphi\in
H^{1}(\mathbb{R}^{3}),\int_{\mathbb{R}^{3}}\varphi^{2}dx=m}E(\varphi)$ (1.5)

does not have a minimizer. The first objective of this paper is to examine the
existence and nonexistence of a minimizer of (1.5) on the lattice graph
$\mathbb{Z}^{3}$.

On the other hand, for $\Omega\in\mathbb{R}^{3}$, we consider the Gamow’s
liquid drop model

$E_{V}=\inf\left\\{E(\Omega):|\Omega|=V\right\\},$ (1.6)

where $E(\Omega)$ is given in (1.2). If the Coulomb term is not present in
(1.1) and (1.2), When the volume constraint $m$ is large, a Modica-Mortola
type result establishes a connection between the two functionals via Gamma
convergence, as stated in [38]. Gamow initially proposed problem (1.6) in [13]
with the intention of uncovering fundamental properties of atoms and offering
a straightforward nuclear fission model. Out of all the measurable sets with
given volume, a ball minimizes the perimeter term while maximizing the Coulomb
term. This highlights the significance of examining the interplay between the
perimeter term and Coulomb term. In [5], the authors proved the existence of a
minimizer of (1.2) as $V$ is small. The paper [28] proved the nonexistence of
the minimizer as $V$ is large by a deep techniques in geometric measure
theory. Also, Lu and Otto [37] proved the nonexistence results by using
different arguments. Furthermore, this type of model is present in the study
of diblock copolymer melt models, Thomas-Fermi type models, and various other
systems, for instance, see [1, 3, 4, 6, 7, 9, 27, 39, 40, 46] and the
references therein.

In this paper, we investigate the existence and nonexistence of a minimzer of
the energy (1.1) and (1.2) on the lattice graph $\mathbb{Z}^{3}.$ Recently,
there has been growing interest among scholars from diverse fields in graph
analysis. Specifically, researchers have explored the extension of Sobolev
inequalities and sharp constants to finite graphs, as evidenced by studies
such as [8, 23, 41, 45]. On the other hand, the study of partial differential
equations on graphs is a fascinating area of research. Notably, the papers
[15, 16, 17] have contributed significantly to this field through a series of
works. Their research employs variational methods to establish the existence
of Yamabe type equations and some nonlinear elliptic equations on graph. The
heat equation has also been a topic of interest, with several researchers
exploring its existence, uniqueness, and blow-up properties. Notable studies
in this area include [29, 57]. For further related research, interested
readers may refer to [14, 18, 19, 20, 21, 22, 26, 53] and the corresponding
references.

To present our findings lucidly, the following is a description of the basic
setting for graphs. Consider a simple, undirected, and locally finite graph
$\mathbb{G}=(\mathbb{V},\mathbb{E})$, where $\mathbb{V}$ represents the set of
vertices and $\mathbb{E}$ represents the set of edges. Vertices $x$ and $y$
are considered neighbours, written as $x\sim y$, if there is an edge
connecting them, meaning that $\\{x,y\\}\in\mathbb{E}$. A graph is considered
locally finite if each vertex has a finite number of neighbours. Denote the
set of functions on $\mathbb{G}$ by $C(\mathbb{V})$. The Laplacian on
$\mathbb{G}=(\mathbb{V},\mathbb{E})$ is defined as, for any function $u\in
C(\mathbb{V})$ and $x\in\mathbb{V}$,

$\Delta u(x)=\sum_{y\in\mathbb{V},~{}y\sim x}\left(u(x)-u(y)\right).$

For any function $u,v\in C(\mathbb{V})$, we can define the associated gradient
by

$\Gamma(u,v)(x):=\sum_{y\in\mathbb{V},~{}y\sim
x}\frac{1}{2}(u(y)-u(x))(v(y)-v(x)).$

In particular, let $\Gamma(u)=\Gamma(u,u)$ for simplicity. We denote the
length of its gradient by

$|\nabla
u|(x):=\sqrt{\Gamma(u)(x)}=\frac{1}{2}\left(\sum_{y\in\mathbb{V},~{}y\sim
x}(u(y)-u(x))^{2}\right)^{1/2}.$

Let $\mu$ be the counting measure on $\mathbb{V}$, i.e., for any subset
$A\subset\mathbb{V}$, $\mu(A)=\\#\\{x:x\in A\cap\mathbb{V}\\}$. For any
function $f$ on $\mathbb{V}$, we write

$\int_{\mathbb{V}}fd\mu:=\sum_{x\in\mathbb{V}}f(x),$

whenever it makes sense. The corresponding energy is defined as

$\mathcal{E}(u):=\int_{\mathbb{V}}|\nabla
u|^{2}d\mu=\frac{1}{2}\sum_{x\in\mathbb{V}}\sum_{y\in\mathbb{V},~{}y\sim
x}(u(y)-u(x))^{2}.$

By definition, $H^{1}(\mathbb{V})$ is a Hilbert space with the inner product
given by

$\langle
u,v\rangle:=\int_{\mathbb{V}}\left(\Gamma(u,v)+uv\right)d\mu=\frac{1}{2}\sum_{x\in\mathbb{V}}\sum_{y\in\mathbb{V},~{}y\sim
x}(u(y)-u(x))(v(y)-v(x))+\sum_{y\in\mathbb{V}}u(y)v(y).$

The lattice graph $\mathbb{Z}^{3}$ is a very representative locally finite
graph, which consists of a group of vertices that are defined as follows: the
set of vertices

$\mathbb{V}=\left\\{x=(x_{1},x_{2},x_{3}):x_{j}\in\mathbb{Z},~{}1\leq j\leq
3\right\\}$ (1.7)

and the set of edges

$\mathbb{E}=\left\\{\\{x,y\\}:x,y\in\mathbb{Z}^{3},~{}\sum_{j=1}^{3}|x_{j}-y_{j}|=1\right\\}.$
(1.8)

For a subset $\Omega$ of $\mathbb{Z}^{3}$, we can define its volume and
boundary area based on the counting measure as follows:

$|\Omega|:=\\#\\{x:x\in\Omega\\}\quad\mathrm{and}\quad|\partial\Omega|:=\left|\left\\{{y\in\Omega,~{}\exists
x\notin\Omega\,\,\,\textrm{such that}\ \\{x,y\\}\in E}\right\\}\right|.$

Let $C(\mathbb{Z}^{3})$ be the set of real-valued functions on
$\mathbb{Z}^{3}$ and $\ell^{p}(\mathbb{Z}^{3})$, $p\in[1,\infty]$, be the
space of $\ell^{p}$ summable functions on $\mathbb{Z}^{3}$ w.r.t. the counting
measure. We write $\|\cdot\|_{p}$ as the $\ell^{p}(\mathbb{Z}^{3})$ norm,
i.e.,

$\|u\|_{p}=\begin{cases}\sum\limits_{x\in\mathbb{Z}^{3}}|u(x)|^{p},\ &1\leq
p<\infty,\\\ \sup\limits_{x\in\mathbb{Z}^{3}}|u(x)|,&p=\infty.\end{cases}$

By definition, we know that

$\|u\|_{H^{1}(\mathbb{Z}^{3})}^{2}=\sum_{y\in\mathbb{Z}^{3}}(|\nabla
u(y)|^{2}+|u(y)|^{2}).$

From [24, The inequality (5)], we have the following equivalent norm

$\|u\|_{H^{1}(\mathbb{Z}^{3})}\simeq\|u\|_{\ell^{2}(\mathbb{Z}^{3})}=\sum_{y\in\mathbb{Z}^{3}}|u(y)|^{2}.$
(1.9)

In the lattice graph setting, the correponding version of energy functional
(1.1) and (1.2) are the following respectively

$E(\varphi)=\sum_{y\in\mathbb{Z}^{3}}\left(|\nabla\varphi(y)|^{2}+\varphi^{\frac{10}{3}}(y)-\varphi^{\frac{8}{3}}(y)\right)+\sum_{x,y\in\mathbb{Z}^{3}\atop~{}\
y\neq x\hfill}\frac{{\varphi}^{2}(x){\varphi}^{2}(y)}{|x-y|},$ (1.10)

and

$E(\Omega)=|\partial\Omega|+\sum_{x,y\in\Omega\atop~{}y\neq
x\hfill}\frac{1}{|x-y|}.$ (1.11)

The first two results concern the functional (1.10) by establishing the
existence of minimizer for small mass and non-existence of minimizer for large
mass.

######  Theorem 1.1.

If $m$ is small enough, then the variational problem

$I_{0}^{TFDW}=\inf_{\sum\limits_{y\in\mathbb{Z}^{3}}{\varphi}^{2}(y)=m}E(\varphi)$
(1.12)

has a minimizer.

######  Theorem 1.2.

There exists a positive constant $m_{0}$ such that the variational problem

$I_{0}^{TFDW}=\inf\limits_{\sum\limits_{y\in\mathbb{Z}^{3}}{\varphi}^{2}(y)=m}E(\varphi)$
(1.13)

does not have a minimizer if $m>m_{0}$.

######  Remark 1.3.

Drawing on the perspective from the continuous case presented in [37], Theorem
1.2 can be understood as asserting that, from a physical standpoint, the
vacuum can only bind a limited number of electrons in the TFDW model on a
lattice graph. In turn, this result provides affirmation for the ”ionization
conjecture” in the specific scenario of $Z=0$ in the lattice graph case.

Consider (1.11), we have the following non-existence result

######  Theorem 1.4.

For a subset $\Omega\subset\mathbb{Z}^{3}$, there exists a positive constant
$V_{0}$, such that the variational problem

$I_{0}^{TFDW}(\Omega)=\inf\limits_{|\Omega|=V}E(\Omega)$ (1.14)

does not have a minimizer if $V>V_{0}$.

Theorems 1.1-1.4 can be viewed as counterparts to the results presented in [5]
and [37] in the discrete setting. Although the results obtained in this work
align with those of the continuous case and are conceptually natural when
extending from Euclidean space to lattice graph. Nevertheless, the proof
methodologies differ due to the distinct characteristics of graph structures
in comparison to continuous spaces. Specifically, in the proof of Lemma 4.1,
while a similar conclusion as (4.4) can be easily established using the
sophisticated scaling function $\psi_{\eta}$ in the continuous setting, as
demonstrated in [37]. It seems that equivalent techniques are not available in
the lattice graph setting. Instead, we introduce the test function (4.5) as a
replacement for the scaling function $\psi_{\eta}$ in order to establish a
similar conclusion in the discrete case. This highlights a notable distinction
between the discrete and continuous cases.

The remainder of this paper is structured as follows. In section 2, we
introduce some fundamental notations and present the formula for the measure
of the ball on the lattice graph $\mathbb{Z}^{3}$. In section 3, we establish
the proof for Theorem 1.1. Section 4 is dedicated to proving Theorem 1.2.
Finally, in section 5, we provide the proof for Theorem 1.4.

Notation

Throughout the paper the letter $C$ will stand for positive constants which
are allowed to vary among different formulas and even within the same lines.

##  2\. Preliminary results

This section is devoted to provide some definitions of the integer lattice
graph $\mathbb{Z}^{3}$ and useful conclusions appeard in previous works. A
subset $\Omega\subset\mathbb{Z}^{3}$ is considered connected if there exists a
path between any two vertices $x$ and $y$ in $\Omega$, denoted by
$\\{x_{i}\\}_{i=1}^{l}\subset\Omega$, such that $x_{1}=x$, $x_{l}=y$ and
$x_{j}\sim x_{j+1}$ for any $1\leq j\leq l-1$. The distance $d(x,y)$ between
two vertices $x$ and $y$ in $\mathbb{Z}^{3}$ is defined as the infimum of $l$
such that $x=x_{1}\sim x_{2}\sim\cdots\sim x_{l}=y$, i.e.,

$d(x,y)=\inf\\{l:x=x_{1}\sim x_{2}\sim\cdots\sim x_{l}=y\\}.$

The diameter of the $\Omega$ is defined as

$\textrm{diam}(\Omega):=\sup_{x,y\in\mathbb{Z}^{3}}d(x,y).$

For any positive integer $R$, the ball on the integer lattice graph is defined
by

$B_{R}(0)=\left\\{x\in\mathbb{Z}^{3},~{}d(x,0)\leq R\right\\},$

and its boundary is

$\partial B_{R}(0)=\left\\{x\in\mathbb{Z}^{3},~{}d(x,0)=R\right\\}.$

In the following lemma we compute the volume and surface area of $B_{R}(0)$ in
the discrete setting

######  Lemma 2.1.

For any positive integer $R$, then we have $|\partial B_{R}|=4R^{2}+2$ and
$|B_{R}|=\frac{4R^{3}+6R^{2}+8R+3}{3}$.

###### Proof.

According to the definition of the distance for the lattice graph, we see that
for any $x=(x_{1},x_{2},x_{3})\in\partial B_{R}(0)$, it holds that

$|x_{1}|+|x_{2}|+|x_{3}|=R.$

Since $\partial B_{R}(0)$ is symmetric, we divide it into three parts as
depicted in Figure 1.

Figure 1. The sketch of $\partial B_{R}$.

($i$) If $x$ lies at the intersection of the axis and the boundary of the
circle with radius $R$, as presented in Figure 1 $(a)$, then there is exactly
one value of $x_{i}$ that is not equal to zero. There are six possible
options:

$(\pm R,0,0),(0,\pm R,0),(0,0,\pm R).$

($ii$) If $x$ lies at the intersection of the plane and the boundary of the
sphere with radius $R$, as illustrated in Figure 1 $(b)$, then $x_{i}$ can
only be zero for one particular index $i$. Without loss of generality, we
assume that $x_{1}=0$, which implies $|x_{2}|+|x_{3}|=R$ and $x_{2},x_{3}\neq
0$. It can be easily demonstrated that there are $4C_{R-1}^{1}=4(R-1)$
possible choices in this case. While if $x_{2}=0$ or $x_{3}=0$, and the
remaining two components are nonzero, we can obtain an additional $8(R-1)$
points.

($iii$) If $x$ is located at the intersection of the quadrant and the boundary
of the sphere with radius $R$, as depicted in Figure 1 $(c)$, none of the
components $x_{i}$ can be zero. In other words, we have
$|x_{1}|+|x_{2}|+|x_{3}|=R$, where $x_{1},x_{2},x_{3}\neq 0$. In this case,
there are $8C_{R-1}^{2}=4(R^{2}-3R+2)$ points that satisfy this condition.

Summarize the discussion from ($i$) to ($iii$), we have

$\begin{split}|\partial
B_{R}|=\left|\\{x\in\mathbb{Z}^{3}\mid|x_{1}|+|x_{2}|+|x_{3}|=R\\}\right|=6+12(R-1)+4(R^{2}-3R+2)=4R^{2}+2.\end{split}$

Therefore, one gets

$\begin{split}|B_{R}|&=1+\sum\limits_{\rm{i}=1}^{R}|\partial
B_{i}|=1+\sum\limits_{\rm{i}=1}^{R}{(4i^{2}+2)}=1+2R+4\sum\limits_{\rm{i}=1}^{R}i^{2}=\frac{4R^{3}+6R^{2}+8R+3}{3}.\end{split}$

This completes the proof. ∎

Next let us recall the following two discrete Brézis-Lieb type Lemmas.

######  Lemma 2.2.

$($[23, Lemma 10]$)$ Let $\Omega\subset\mathbb{Z}^{N}$ be a domain and
$\\{u_{n}\\}\subset\ell^{q}(\Omega)$ with $1\leq q<\infty$. If $\\{u_{n}\\}$
is bounded in $\ell^{q}(\Omega)$ and $u_{n}\rightarrow u$ pointwise on
$\Omega$ as $n\to\infty$, then

$\lim_{n\rightarrow\infty}\left(\|u_{n}\|_{\ell^{q}(\Omega)}^{q}-\|u_{n}-u\|_{\ell^{q}(\Omega)}^{q}\right)=\|u\|_{\ell^{q}(\Omega)}^{q}.$
(2.1)

######  Lemma 2.3.

$($[55, Lemma 3.7]$)$ Let $1\leq p<\infty$ and the sequence $\\{u_{n}\\}$ is
bounded in $\ell^{\frac{2Np}{N+\alpha}}(\mathbb{Z}^{N})$. Suppose that
$u_{n}\rightarrow u$ pointwise on ${\mathbb{Z}^{N}}$ as $n\to\infty$, then

$\begin{split}\lim\limits_{n\to\infty}\left(\sum_{x,y\in{\mathbb{Z}^{N}}\hfill\atop~{}\
y\neq
x\hfill}\frac{|u_{n}(x)|^{p}|u_{n}(y)|^{p}}{|x-y|^{N-\alpha}}-\sum_{x,y\in{\mathbb{Z}^{N}}\hfill\atop~{}\
y\neq
x\hfill}\frac{|u_{n}(x)-u(x)|^{p}|u_{n}(y)-u(y)|^{p}}{|x-y|^{N-\alpha}}\right)=\sum_{x,y\in{\mathbb{Z}^{N}}\hfill\atop~{}\
y\neq x\hfill}\frac{|u(x)|^{p}|u(y)|^{p}}{|x-y|^{N-\alpha}}.\end{split}$ (2.2)

The subsequent two important inequalities play a crucial role in our argument

######  Lemma 2.4.

$($[23, Lemma 7], [25, Lemma 2.11]$)$ If $u\in\ell^{p}(\mathbb{Z}^{N})$, then
for all $q\geq p$,

$\|u\|_{\ell^{q}(\mathbb{Z}^{N})}\leq\|u\|_{\ell^{p}(\mathbb{Z}^{N})}.$ (2.3)

######  Lemma 2.5.

$($[25, Theorem 1.1]$)$ $($Discrete Hardy-Littlewood-Sobolve Inequality$)$ Let
$0<\alpha<N$, $1<r,s<\infty$ and
$\frac{1}{r}+\frac{1}{s}+\frac{{N-\alpha}}{N}\geq 2$. Assume
$f\in\ell^{r}(\mathbb{Z}^{N})$ and $g\in\ell^{s}(\mathbb{Z}^{N})$. Then there
exists a positive constant $C$ depending only on $r,s,\alpha$ such that

$\sum_{x,y\in\mathbb{Z}^{N}\hfill\atop~{}\ y\neq
x\hfill}{\frac{{f(x)g(y)}}{|x-y|^{N-\alpha}}}\leq
C{\|f\|_{{\ell^{r}}(\mathbb{Z}^{N})}}{\|g\|_{{\ell^{s}}(\mathbb{Z}^{N})}}.$
(2.4)

##  3\. The existence of a minimizer for (1.10)

This section is dedicated to presenting the proof of Theorem 1.1. For the
convenience, we denote

$I_{m}:=\inf\left\\{E(\varphi):\varphi\in
H^{1}(\mathbb{Z}^{3}),\sum_{y\in\mathbb{Z}^{3}}{\varphi^{2}(y)}=m\right\\}$

and

$F(s):=s^{\frac{5}{3}}-s^{\frac{4}{3}},\,\,\,D(f,g)=\sum_{x,y\in\mathbb{Z}^{3}\atop~{}\
x\neq y\hfill}\frac{f^{2}(x)g^{2}(y)}{|x-y|}.$

We begin by establishing a property of strict subadditivity for $I_{m}$.

######  Lemma 3.1.

If $m>0$ is small enough, then for any $m_{1}\in(0,m)$ it holds that

$I_{m}<I_{m_{1}}+I_{m-m_{1}}.$ (3.1)

###### Proof.

First, we claim that if $\theta>1$, then the equality

$I_{\theta}m<\theta I_{m}$ (3.2)

holds for $m>0$ is small enough. Indeed, let $\varphi$ satisfying $I_{m}\leq
E(\varphi)\leq I_{m}+\varepsilon$ and $\|\varphi\|_{\ell^{2}}^{2}=m$. By (2.3)
and the Hardy-Littlewood-Sobolve inequality (2.4), we have

$\begin{split}I_{\theta m}\leq
E\left(\sqrt{\theta}\varphi\right)=~{}&\theta\sum_{y\in\mathbb{Z}^{3}}|\nabla\varphi(y)|^{2}+\theta^{\frac{5}{3}}\sum_{y\in\mathbb{Z}^{3}}\varphi^{\frac{10}{3}}(y)-\theta^{\frac{4}{3}}\sum_{y\in\mathbb{Z}^{3}}\varphi^{\frac{8}{3}}(y)+\theta^{2}D(\varphi^{2},\varphi^{2})\\\
=~{}&\theta\sum_{y\in\mathbb{Z}^{3}}\left(|\nabla\varphi(y)|^{2}+\varphi^{\frac{10}{3}}(y)-\varphi^{\frac{8}{3}}(y)\right)+\theta
D(\varphi^{2},\varphi^{2})+(\theta^{\frac{5}{3}}-\theta)\sum_{y\in\mathbb{Z}^{3}}\varphi^{\frac{10}{3}}(y)\\\
&+(\theta-\theta^{\frac{4}{3}})\sum_{y\in\mathbb{Z}^{3}}\varphi^{\frac{8}{3}}(y)+(\theta^{2}-\theta)D(\varphi^{2},\varphi^{2})\\\
\leq~{}&\theta(I_{m}+\varepsilon)+(\theta^{\frac{5}{3}}-\theta)\|\varphi\|_{\ell^{2}(\mathbb{Z}^{3})}^{\frac{10}{3}}+(\theta-\theta^{\frac{4}{3}})\|\varphi\|_{\ell^{2}(\mathbb{Z}^{3})}^{\frac{8}{3}}+(\theta^{2}-\theta)C\|\varphi\|_{\ell^{\frac{12}{5}}(\mathbb{Z}^{3})}^{4}\\\
\leq~{}&\theta(I_{m}+\varepsilon)+(\theta^{\frac{5}{3}}-\theta)m^{\frac{5}{3}}+(\theta-\theta^{\frac{4}{3}})m^{\frac{4}{3}}+(\theta^{2}-\theta)Cm^{2}.\end{split}$

where $C$ is a positive constant. Then we can choose $m$ small enough such
that

$(\theta^{\frac{5}{3}}-\theta)m^{\frac{5}{3}}+(\theta-\theta^{\frac{4}{3}})m^{\frac{4}{3}}+(\theta^{2}-\theta)Cm^{2}<0.$

As a consequence, we get that

$I_{\theta m}<\theta I_{m}.$

We take $\theta>1$ such that $\theta m_{1}=m$. Then
$m-m_{1}=\frac{\theta-1}{\theta}m$. From (3.2) we get that

$I_{m}=I_{\theta m_{1}}<\theta I_{m_{1}}$ (3.3)

and

$I_{m}=I_{\frac{\theta}{\theta-1}(m-m_{1})}<\frac{\theta}{\theta-1}I_{m-m_{1}}.$
(3.4)

Combining with (3.3) and (3.4), it follows that

$I_{m}=\frac{1}{\theta}I_{m}+\frac{\theta-1}{\theta}I_{m}<\frac{1}{\theta}\theta
I_{m_{1}}+\frac{\theta-1}{\theta}\frac{\theta}{\theta-1}I_{m-m_{1}}=I_{m_{1}}+I_{m-m_{1}}.$

Thus we complete the proof. ∎

Now we are ready to give the proof of Theorem 1.1.

###### Proof of Theorem 1.1..

Let $\\{\varphi_{n}\\}\in{\ell^{2}}(\mathbb{Z}^{3})$ be a minimizing sequence
such that

$\|\varphi_{n}\|_{\ell^{2}(\mathbb{Z}^{3})}^{2}=m\quad\mathrm{and}\quad\mathop{\lim}\limits_{n\rightarrow\infty}E(\varphi_{n})=I_{m}.$

Since the sequence $\\{\varphi_{n}\\}$ is bounded in
${H^{1}}(\mathbb{Z}^{3})$, it follows that
${\varphi_{n}}\rightharpoonup\varphi$ in ${H^{1}}(\mathbb{Z}^{3})$ and
${\varphi_{n}}\to\varphi$ pointwise on $\mathbb{Z}^{3}$. Since the $\ell^{2}$
norm is weakly lower semi-continuous, it leads to

$0\leq\|\varphi\|_{\ell^{2}(\mathbb{Z}^{3})}^{2}\leq m.$

Next, we divide the proof into the following three steps

Step 1: The sequence $\\{\varphi_{n}\\}$ is non-vanishing. By (1.9) and the
inequality (2.3), then

$E(\varphi_{n})\geq\sum_{y\in\mathbb{Z}^{3}}|\nabla\varphi_{n}(y)|^{2}-\sum_{y\in\mathbb{Z}^{3}}\varphi^{\frac{8}{3}}_{n}(y)\geq
Cm-m^{\frac{4}{3}},$

where $C$ is some positive constant. Since
$E^{\prime}(\varphi_{n})\varphi_{n}=0$, it follows that

$\sum_{y\in\mathbb{Z}^{3}}|\nabla\varphi_{n}(y)|^{2}+\frac{10}{3}\sum_{y\in\mathbb{Z}^{3}}\varphi^{\frac{10}{3}}_{n}(y)-\frac{8}{3}\sum_{y\in\mathbb{Z}^{3}}\varphi^{\frac{8}{3}}_{n}(y)+4D(\varphi_{n}^{2},\varphi_{n}^{2})=0.$

Thus, by the interpolation inequality, we get

$\begin{split}Cm-m^{\frac{4}{3}}&\leq\lim_{n\to\infty}E(\varphi_{n})=\lim_{n\to\infty}\left(-\frac{7}{3}\sum_{y\in\mathbb{Z}^{3}}\varphi^{\frac{10}{3}}_{n}(y)+\frac{5}{3}\sum_{y\in\mathbb{Z}^{3}}\varphi^{\frac{8}{3}}_{n}(y)-3D(\varphi_{n}^{2},\varphi_{n}^{2})\right)\\\
&\leq\lim_{n\to\infty}\frac{5}{3}\sum_{y\in\mathbb{Z}^{3}}\varphi^{\frac{8}{3}}_{n}(y)\leq\lim_{n\to\infty}\frac{5}{3}\|\varphi_{n}\|_{\ell^{2}(\mathbb{Z}^{3})}^{2}\|\varphi_{n}\|_{\ell^{\infty}(\mathbb{Z}^{3})}^{\frac{2}{3}}=\lim_{n\to\infty}\frac{5}{3}m\|\varphi_{n}\|_{\ell^{\infty}(\mathbb{Z}^{3})}^{\frac{2}{3}},\end{split}$

which implies

$\lim\limits_{n\to\infty}\|\varphi_{n}\|_{\ell^{\infty}(\mathbb{Z}^{3})}>\frac{1}{2}C>0$

provided $m$ is small enough. Since
$\|\varphi_{n}\|_{\ell^{2}(\mathbb{Z}^{3})}^{2}=m$, it follows that the
maximum of $|\varphi_{n}|$ is attainable in a bounded subset of
$\mathbb{Z}^{3}$. We set $\widetilde{\varphi}_{n}(y):=\varphi_{n}(y-{y_{n}})$,
where $y_{n}$ is chosen such that
$|\varphi_{n}(y_{n})|=\max\limits_{y}|\varphi_{n}(y)|$. Taking a subsequence
if necessary, we have

$\widetilde{\varphi}_{n}\to\widetilde{\varphi}\,\,\,\,\,\,\text{pointwise}\,\,\,\text{in}\,\,\,\mathbb{Z}^{3}$

and

$|\widetilde{\varphi}(0)|=\lim_{n\to\infty}\|\widetilde{\varphi}_{n}\|_{\ell^{\infty}(\mathbb{Z}^{3})}>0.$

By translation invariance, we can choose a minimizing sequence such that the
limit $\varphi\neq 0$.

Step 2: We claim that $\|\varphi\|_{\ell^{2}(\mathbb{Z}^{3})}^{2}=m$. Using a
contradiction argument, suppose that

$0<\|\varphi\|_{\ell^{2}(\mathbb{Z}^{3})}^{2}=m_{1}<m.$

By the Lemma 2.2, we get that

$\|\varphi_{n}-\varphi\|_{\ell^{2}(\mathbb{Z}^{3})}^{2}=m-m_{1}.$

From Lemmas 2.2-2.3 and [23, Corollary 11] that

$I_{m}=\lim_{n\rightarrow\infty}E(\varphi_{n})=\lim_{n\rightarrow\infty}E(\varphi_{n}-\varphi)+E(\varphi)\geq
I_{m_{1}}+I_{m-m_{1}}.$

This is in contradiction to the inequality (3.1).

Step 3: $\varphi$ is the minimizer of the functional $E$, i.e.
$E(\varphi)=I_{m}$. According to the infimum of the $I_{m}$ and
$\|\varphi\|_{\ell^{2}(\mathbb{Z}^{3})}^{2}=m$, then $E(\varphi)\geq I_{m}$.
By the inequality (2.3), we get

$\|\varphi_{n}-\varphi\|_{\ell^{\frac{8}{3}}(\mathbb{Z}^{3})}\leq\|\varphi_{n}-\varphi\|_{\ell^{2}(\mathbb{Z}^{3})}.$

Then it follows that
$\mathop{\lim}\limits_{n\to\infty}\|\varphi_{n}\|_{\ell^{\frac{8}{3}}(\mathbb{Z}^{3})}=\|\varphi\|_{\ell^{\frac{8}{3}}(\mathbb{Z}^{3})}$.
Using the Fatou’s Lemma we have

$\begin{split}I_{m}=\lim_{n\rightarrow\infty}E(\varphi_{n})&=\lim_{n\rightarrow\infty}\left(\sum_{y\in\mathbb{Z}^{3}}\left(|\nabla\varphi_{n}(y)|^{2}+F(\varphi_{n}^{2}(y))\right)+D(\varphi_{n}^{2},\varphi_{n}^{2})\right)\\\
&\geq\sum_{y\in\mathbb{Z}^{3}}\left(|\nabla\varphi(y)|^{2}+F(\varphi^{2}(y))\right)+D(\varphi^{2},\varphi^{2})=E(\varphi).\end{split}$

This implies that $E(\varphi)=I_{m}$. Therefore, $\varphi$ is a minimizer of
the functional $E$ provided the mass $m$ is small enough. Hence, we finish the
proof. ∎

##  4\. The nonexistence of a minimizer for (1.10)

In this section, we will establish the proof for Theorem 1.2. The key idea is
to demonstrate that the minimizer $\varphi$ of (1.10) must have a uniform
bound in the $l^{2}(\mathbb{Z}^{3})$ norm, i.e.,
$m:=\sum\limits_{y\in\mathbb{Z}^{3}}\varphi^{2}<C$, where $C$ is a finite
positive constant. Consequently, if $m>C$, the functional (1.10) cannot have
any minimizer. To begin with, we will examine this fact by analyzing the
following elementary property of $\varphi$.

######  Lemma 4.1.

For every $x$ in $\mathbb{Z}^{3}$, then we have
$\varphi(x)\in[0,(\frac{4}{5})^{\frac{3}{2}}]$ and $-\varphi^{2}(x)\lesssim
F({\varphi}^{2})\leq 0$.

###### Proof.

As $E(|\varphi|)\leq E(\varphi)$, we can focus our analysis on the case where
$\varphi\geq 0$. Consider the function $F(s)=s^{\frac{5}{3}}-s^{\frac{4}{3}}$,
we have $\min\limits_{s\geq 0}F(s)=F\left((\frac{4}{5})^{3}\right)$. We define
$\varphi_{1}=\min\left\\{\varphi,(\frac{4}{5})^{\frac{3}{2}}\right\\}$ and
$m_{1}:=\sum\limits_{y\in\mathbb{Z}^{3}}\varphi_{1}^{2}(y)$, where it is
evident that $m_{1}\leq m$. Next, we will proceed with our argument in the
following three steps:

Step 1. We claim that

$E(\varphi)\leq E(\varphi_{1})+I_{m-m_{1}}.$ (4.1)

Indeed, for any $\varepsilon>0$, choosing $\psi$ with
$\sum\limits_{y\in\mathbb{Z}^{3}}\psi^{2}(y)=m-m_{1}$ and $I_{m-m_{1}}\leq
E(\psi)\leq I_{m-m_{1}}+\varepsilon$. We then define

$\psi_{L}(x):=C_{L}[\varphi_{1}(x)+\psi(x+Le_{1})],$

where $e_{1}:=(1,0,0)\in\mathbb{Z}^{3}$ and $C_{L}$ is a normalization
constant such that

$\sum\limits_{y\in\mathbb{Z}^{3}}\psi_{L}^{2}(y)=\sum\limits_{y\in\mathbb{Z}^{3}}\varphi^{2}(y)=m.$

Moreover, it is not difficult to check that

$\begin{split}\sum\limits_{y\in\mathbb{Z}^{3}}\psi_{L}^{2}(y)&=C_{L}^{2}\left(\sum\limits_{y\in\mathbb{Z}^{3}}\varphi_{1}^{2}(y)+\sum\limits_{y\in\mathbb{Z}^{3}}\psi^{2}(y)+2\sum\limits_{y\in\mathbb{Z}^{3}}\varphi_{1}(y)\psi\left(y+Le_{1}\right)\right)\\\
&=C_{L}^{2}\left(m+2\sum_{y\in\mathbb{Z}^{3}}\varphi_{1}(y)\psi\left(y+Le_{1}\right)\right).\end{split}$

On the other hand, it is worth noting that $C_{L}\to 1$ as $L$ approaches
infinity. Consequently

$\limsup_{L\to\infty}E(\psi_{L})=E(\varphi_{1})+E(\psi)\leq
E(\varphi_{1})+I_{m-m_{1}}+\varepsilon.$ (4.2)

Using the fact that $\varphi$ is the minimizer of (1.10), we get that
$E(\varphi)\leq E(\psi_{L})$ for any $L$. Sending $\varepsilon\to 0$ in (4.2),
we arrive at (4.1).

Step 2. In this step we shall show that

$E(\varphi)\leq E(\varphi_{1}).$ (4.3)

In order to achieve this inequality, it suffices to verify that

$I_{m-m_{1}}=\inf_{\sum\limits_{y\in\mathbb{Z}^{3}}\psi^{2}(y)=m-m_{1}}E(\psi)\leq
0.$ (4.4)

Now we define the function $\psi_{n}$ in the following way:

$\psi_{n}(x):=\begin{cases}(n-\ell+1)d,&\ell\leq n,\\\ 0,&\ell>n,\end{cases}$
(4.5)

where

$\ell=d(x,0)\quad\text{and}\quad
d^{2}=\frac{m-m_{1}}{\sum\limits_{\ell=1}^{n}(n-\ell+1)^{2}(4\ell^{2}+2)+(n+1)^{2}}.$

To obtain a more precise comprehension of the meaning of $\psi_{n}$, we
provide a schematic diagram in Figure 3.

Figure 2. The sketch of $\psi_{n}$

With the help of Lemma 2.1, we can easily compute that

$\displaystyle\sum\limits_{y\in\mathbb{Z}^{3}}\psi_{n}^{2}(y)={}$
$\displaystyle(n+1)^{2}d^{2}+\sum_{\ell=1}^{n}\sum_{y\in\partial
B_{\ell}(0)\cap\mathbb{Z}^{3}}\psi_{n}^{2}(y)$ $\displaystyle={}$
$\displaystyle\left[(n+1)^{2}+\sum_{\ell=1}^{n}(n-\ell+1)^{2}(4\ell^{2}+2)\right]d^{2}=m-m_{1}.$

Regarding the gradient term, using Lemma 2.1 we have

$\displaystyle\sum_{y\in\mathbb{Z}^{3}}|\nabla\psi_{n}|^{2}$
$\displaystyle=\frac{1}{2}\sum_{y\in\mathbb{Z}^{3}}\sum_{z\sim y,\atop
z\in\mathbb{Z}^{3}}(\psi_{n}(x)-\psi_{n}(y))^{2}\leq
Cd^{2}\sum_{\ell=1}^{n+1}(4\ell^{2}+2)$ (4.6)
$\displaystyle=\frac{C(m-m_{1})\sum\limits_{\ell=1}^{n+1}(4\ell^{2}+2)}{\sum\limits_{\ell=1}^{n}(n-\ell+1)^{2}(4\ell^{2}+2)+(n+1)^{2}}.$

After direct computation, we get that

$a_{n}:=\sum_{\ell=1}^{n+1}(4\ell^{2}+2)=\frac{2}{3}(2n^{3}+9n^{2}+16n+9),$

and

$b_{n}:=\sum_{\ell=1}^{n}(n-\ell+1)^{2}(4\ell^{2}+2)+(n+1)^{2}=\frac{1}{15}(2n^{5}+10n^{4}+30n^{3}+50n^{2}+43n+15).$

Therefore, we have $\lim\limits_{n\to+\infty}\frac{a_{n}}{b_{n}}=0$ and it
implies that

$\lim_{n\to+\infty}\sum_{y\in\mathbb{Z}^{3}}|\nabla\psi_{n}(y)|^{2}=0.$ (4.7)

While for the nonlinear term, after a straightforward calculation we obtain
that

$\displaystyle\sum_{y\in\mathbb{Z}^{3}}\psi_{n}^{\frac{10}{3}}(y)=$
$\displaystyle\sum_{\ell=1}^{n}(n-\ell+1)^{\frac{10}{3}}d^{\frac{10}{3}}(4\ell^{2}+2)+(n+1)^{\frac{10}{3}}d^{\frac{10}{3}}$
$\displaystyle=$
$\displaystyle~{}d^{\frac{10}{3}}\sum_{\ell=1}^{n}\left((n-\ell+1)^{2}(4\ell^{2}+2)^{\frac{3}{5}}\right)^{\frac{5}{3}}+(n+1)^{\frac{10}{3}}d^{\frac{10}{3}}$
$\displaystyle\leq$
$\displaystyle\left(d^{2}\sum_{\ell=1}^{n}(n-\ell+1)^{2}(4\ell^{2}+2)^{\frac{3}{5}}\right)^{\frac{5}{3}}+\left((n+1)^{2}d^{2}\right)^{\frac{5}{3}}$
$\displaystyle=$
$\displaystyle\left(\frac{(m-m_{1})\sum\limits_{\ell=1}^{n}(n-\ell+1)^{2}(4\ell^{2}+2)^{\frac{3}{5}}}{\sum\limits_{\ell=1}^{n}(n-\ell+1)^{2}(4\ell^{2}+2)+(n+1)^{2}}\right)^{\frac{3}{5}}+\left(\frac{(m-m_{1})(n+1)^{2}}{\sum\limits_{\ell=1}^{n}(n-\ell+1)^{2}(4\ell^{2}+2)+(n+1)^{2}}\right)^{\frac{5}{3}}.$

We set

$c_{n}:=\sum_{\ell=1}^{n}(n-\ell+1)^{2}(4\ell^{2}+2)^{\frac{3}{5}}\quad\text{and}\quad
d_{n}=(n+1)^{2}.$

We notice that

$c_{n}\leq(4n^{2}+2)^{\frac{3}{5}}\sum_{\ell=1}^{n}\ell^{2}\leq\frac{1}{3}(n+1)^{3}(4n^{2}+2)^{\frac{3}{5}}.$

Then one can easily see that

$\lim_{n\to\infty}\frac{c_{n}}{b_{n}}=0\quad\mathrm{and}\quad\lim_{n\to\infty}\frac{d_{n}}{b_{n}}=0.$
(4.8)

This implies that

$\sum_{y\in\mathbb{Z}^{3}}\psi_{n}^{\frac{10}{3}}(y)\to 0\quad\mathrm{as}\quad
n\to\infty.$

Similarly, we can prove

$\sum_{y\in\mathbb{Z}^{3}}\psi_{n}^{\frac{8}{3}}(y)\to 0\quad\mathrm{as}\quad
n\to\infty.$

Therefore,

$\lim_{n\to\infty}\sum_{y\in\mathbb{Z}^{3}}F(\psi_{n}^{2}(y))=0.$ (4.9)

Finally, we study the nonlocal term. From the Hardy-Littlewood Soblve
inequality (2.4) we derive that

$\sum_{y\in\mathbb{Z}^{3}}\sum_{x\neq y,\atop
x\in\mathbb{Z}^{3}}\frac{\psi_{n}^{2}(y)\psi_{n}^{2}(x)}{|y-x|}\leq
C\|\psi_{n}\|_{\ell^{\frac{12}{5}}(\mathbb{Z}^{3})}^{4}.$

By a direct computation, one has

$\begin{split}\sum_{y\in\mathbb{Z}^{3}}\psi_{n}^{\frac{12}{5}}(y)=&\left(\sum_{\ell=1}^{n}(n-\ell+1)^{\frac{12}{5}}(4\ell^{2}+2)+(n+1)^{\frac{12}{5}}\right)d^{\frac{12}{5}}\\\
\leq&\left(\sum_{\ell=1}^{n}(n-\ell+1)^{2}(4\ell^{2}+2)^{\frac{5}{6}}d^{2}\right)^{\frac{6}{5}}+\left((n+1)^{2}d^{2}\right)^{\frac{6}{5}}\\\
=&\left(\frac{(m-m_{1})\sum\limits_{\ell=1}^{n}(n-\ell+1)^{2}(4\ell^{2}+2)^{\frac{5}{6}}}{\sum\limits_{\ell=1}^{n}(n-\ell+1)^{2}(4\ell^{2}+2)+(n+1)^{2}}\right)^{\frac{6}{5}}+\left(\frac{(m-m_{1})(n+1)^{2}}{\sum\limits_{\ell=1}^{n}(n-\ell+1)^{2}(4\ell^{2}+2)+(n+1)^{2}}\right)^{\frac{6}{5}}.\end{split}$

We define

$e_{n}:=\sum_{\ell=1}^{n}(n-\ell+1)^{2}(4\ell^{2}+2)^{\frac{5}{6}}.$

Obviously, we have

$e_{n}\leq(4n^{2}+2)^{\frac{5}{6}}\sum_{\ell=1}^{n}\ell^{2}\leq\frac{1}{3}(n+1)^{3}(4n^{2}+2)^{\frac{5}{6}},$

then

$\lim_{n\to\infty}\frac{e_{n}}{b_{n}}=0.$

Together with (4.8) we have

$\lim_{n\to\infty}\sum_{y\in\mathbb{Z}^{3}}\sum_{x\neq y,\atop
x\in\mathbb{Z}^{3}}\frac{\psi_{n}^{2}(y)\psi_{n}^{2}(x)}{|x-y|}=0.$ (4.10)

From (4.7)-(4.10), one deduces that $E(\psi_{n})\to 0$ as $n\to\infty$, which
together with (4.1) implies that

$E(\varphi)\leq E(\varphi_{1}).$ (4.11)

Step 3. We claim that

$\varphi(x)\leq\left(\frac{4}{5}\right)^{\frac{3}{2}},\quad\forall
x\in\mathbb{Z}^{3}.$

We prove the above claim by contradiction. Suppose that there exists a point
$x_{0}\in\mathbb{Z}^{3}$ such that
$\varphi(x_{0})>(\frac{4}{5})^{\frac{3}{2}}$. Then, one has

$\displaystyle\sum_{y\in\mathbb{Z}^{3}}F(\varphi_{1}^{2}(y))$
$\displaystyle=\sum_{\left\\{y\in\mathbb{Z}^{3}:\varphi(y)>(4/5)^{\frac{3}{2}}\right\\}}F(\varphi_{1}^{2}(y))+\sum_{\left\\{y\in\mathbb{Z}^{3}:\varphi(y)\leq(4/5)^{\frac{3}{2}}\right\\}}F(\varphi_{1}^{2}(y))$
$\displaystyle=\sum_{\left\\{y\in\mathbb{Z}^{3}:\varphi(y)>(4/5)^{\frac{3}{2}}\right\\}}F((4/5)^{3})+\sum_{\left\\{y\in\mathbb{Z}^{3}:\varphi(y)\leq(4/5)^{\frac{3}{2}}\right\\}}F(\varphi^{2}(y))$
$\displaystyle<\sum_{\left\\{y\in\mathbb{Z}^{3}:\varphi(y)>(4/5)^{\frac{3}{2}}\right\\}}F(\varphi^{2}(y))+\sum_{\left\\{y\in\mathbb{Z}^{3}:\varphi(y)\leq(4/5)^{\frac{3}{2}}\right\\}}F(\varphi^{2}(y))$
$\displaystyle=\sum_{y\in\mathbb{Z}^{3}}F(\varphi^{2}(y)).$

Given that $|\nabla\varphi_{1}|\leq|\nabla\varphi|$ and
$D(\varphi_{1}^{2},\varphi_{1}^{2})\leq D(\varphi^{2},\varphi^{2})$, we can
conclude that $E(\varphi)>E(\varphi_{1})$. However, this contradicts the
inequality (4.11). Thus, it proves
$0\leq\varphi\leq(\frac{4}{5})^{\frac{3}{2}}$, and the remaining part of the
lemma can be easily verified. ∎

The following lemma presents a crucial estimate of the minimizer $\varphi$.

######  Lemma 4.2.

Given any $x\in\mathbb{Z}^{3}$ and two radii $r,R>0$ such that $R>r+1$, the
following holds:

$\sum_{y\in(B_{r+1}(x)\cup
B_{r}(x))\cap\mathbb{Z}^{3}}\varphi^{2}(y)\geq\frac{1}{3(R+r)}\sum_{y\in
B_{r}(x)\cap\mathbb{Z}^{3}}\varphi^{2}(y)\sum_{y\in(B_{R}(x)\setminus
B_{r+1}(x))\cap\mathbb{Z}^{3}}\varphi^{2}(y).$

###### Proof.

Without loss of generality, we can assume that $x=0$ due to the translational
invariance. Let us define the cutoff functions
$\chi_{1},\chi_{2}:\mathbb{Z}^{3}\to[0,1]$ as

$\displaystyle\chi_{1}(x)=\begin{cases}1,~{}&d(x,0)\leq r,\\\ 0,~{}&d(x,0)\geq
r+1,\end{cases}\quad\chi_{2}(x)=\begin{cases}0,~{}&d(x,0)\leq r,\\\
1,~{}&d(x,0)\geq r+1.\end{cases}$

Let $\varphi_{1}=\chi_{1}\varphi$ and $\varphi_{2}=\chi_{2}\varphi.$ Similar
to the proof of (4.1), we can obtain the following equality:

$E(\varphi)\leq E(\varphi_{1})+E(\varphi_{2}).$

This implies

$\begin{split}&\sum_{y\in\mathbb{Z}^{3}}\left(|\nabla\varphi_{1}(y)|^{2}+|\nabla\varphi_{2}(y)|^{2}-|\nabla\varphi(y)|^{2}\right)+\sum_{y\in\mathbb{Z}^{3}}\left(F({\varphi^{2}_{1}(y)})+F({\varphi^{2}_{2}(y)})-F({\varphi^{2}(y)})\right)\\\
&\geq
D(\varphi^{2},\varphi^{2})-D(\varphi_{1}^{2},\varphi_{1}^{2})-D(\varphi_{2}^{2},\varphi_{2}^{2}).\end{split}$
(4.12)

We shall analyze the equation (4.12) by studying the gradient term first. If
$d(x,0)<r$ or $d(x,0)>r+1$. Then

$|\nabla\varphi_{1}|^{2}+|\nabla\varphi_{2}|^{2}-|\nabla\varphi|^{2}=0.$

If $d(x,0)=r$,

$\begin{split}|\nabla\varphi_{1}|^{2}+|\nabla\varphi_{2}|^{2}-|\nabla\varphi|^{2}=\sum_{y\sim
x,\atop
d(y,0)=r+1}\frac{1}{2}\left(\varphi^{2}(x)+\varphi^{2}(y)-(\varphi(y)-\varphi(x))^{2}\right)=\sum_{y\sim
x,\atop d(y,0)=r+1}\varphi(x)\varphi(y).\end{split}$

If $d(x,0)=r+1$,

$\begin{split}|\nabla\varphi_{1}|^{2}+|\nabla\varphi_{2}|^{2}-|\nabla\varphi|^{2}=\sum_{y\sim
x,\atop
d(y,0)=r}\frac{1}{2}\left(\varphi^{2}(y)+\varphi^{2}(x)-(\varphi(y)-\varphi(x))^{2}\right)=\sum_{y\sim
x,\atop d(y,0)=r}\varphi(x)\varphi(y).\end{split}$

Therefore, from the above three equations we get that

$\begin{split}\sum_{y\in\mathbb{Z}^{3}}\left(|\nabla\varphi_{1}(y)|^{2}+|\nabla\varphi_{2}(y)|^{2}-|\nabla\varphi(y)|^{2}\right)=&\sum_{y\in\partial
B_{r}\cap\mathbb{Z}^{3}}\sum_{x\sim y,\atop x\in{\partial
B}_{r+1}\cap\mathbb{Z}^{3}\hfill}\varphi(x)\varphi(y)\\\ &+\sum_{y\in\partial
B_{r+1}\cap\mathbb{Z}^{3}}\sum_{x\sim y,\atop x\in{\partial
B}_{r}\cap\mathbb{Z}^{3}\hfill}\varphi(x)\varphi(y)\\\
\leq&~{}\sum_{y\in(B_{r+1}\setminus
B_{r-1})\cap\mathbb{Z}^{3}}6\varphi^{2}(y),\end{split}$ (4.13)

where the number $6$ appears due to the fact that each vertex point has 6
neighboring points at most. For the second term on the left hand side of
(4.12), one can easily check that

$\begin{split}F(\varphi_{1}^{2})+F(\varphi_{2}^{2})-F(\varphi^{2})=0.\end{split}$

As a consequence,

$\sum_{y\in\mathbb{Z}^{3}}\left(F(\varphi_{1}^{2}(y))+F(\varphi_{2}^{2}(y))-F(\varphi^{2}(y))\right)=0.$
(4.14)

Consider the term on the right hand side of (4.12), after a direct computation
we get that

$\begin{split}D(\varphi^{2},\varphi^{2})-D(\varphi_{1}^{2},\varphi_{1}^{2})-D(\varphi_{2}^{2},\varphi_{2}^{2})=~{}&2\sum_{x\in\mathbb{Z}^{3}}\sum_{y\neq
x,\atop
y\in\mathbb{Z}^{3}\hfill}\frac{\varphi_{1}^{2}(x)\varphi_{2}^{2}(y)}{|x-y|}=2\sum_{x\in\mathbb{Z}^{3}}\sum_{y\neq
x,\atop
y\in\mathbb{Z}^{3}\hfill}\frac{\varphi_{1}^{2}(x)\varphi_{2}^{2}(y)}{|x-y|}\\\
\geq~{}&2\sum_{x\in B_{r}\cap\mathbb{Z}^{3}}\sum_{y\neq x,\atop
y\in(B_{R}\setminus
B_{r+1})\cap\mathbb{Z}^{3}\hfill}\frac{\varphi_{1}^{2}(x)\varphi_{2}^{2}(y)}{|x-y|}\\\
=~{}&2\sum_{x\in B_{r}\cap\mathbb{Z}^{3}}\sum_{y\neq x,\atop
y\in(B_{R}\setminus
B_{r+1})\cap\mathbb{Z}^{3}\hfill}\frac{\varphi^{2}(x)\varphi^{2}(y)}{|x-y|}\\\
\geq~{}&\frac{2}{r+R}\sum_{x\in
B_{r}\cap\mathbb{Z}^{3}}\varphi^{2}(x)\sum_{y\in(B_{R}\setminus
B_{r+1})\cap\mathbb{Z}^{3}}\varphi^{2}(y).\end{split}$ (4.15)

By (4.12) and (4.13)-(4.15), we complete the proof. ∎

The following lemma provides an upper bound for the mass
$m=\sum\limits_{y\in\mathbb{Z}^{3}}\varphi^{2}(y)$.

######  Lemma 4.3.

For any ball $B_{R_{0}}(x)$ with $R_{0}\geq 4$ and $\sum\limits_{y\in
B_{R_{0}}(x)}\varphi^{2}(y)\geq C_{0}$, there exists a universal constant
$C_{0}$ such that:

$m=\sum_{y\in\mathbb{Z}^{3}}\varphi^{2}(y)\leq 2\sum_{y\in
B_{2R_{0}}(x)\cap\mathbb{Z}^{3}}\varphi^{2}(y).$

###### Proof.

Assuming $x=0$ for simplicity of notation, we consider the balls $B_{r}$ that
are concentric with $B_{R_{0}}$. For $r\geq 1$, we define

$S(r)=\sum_{y\in B_{r}(x)\cap\mathbb{Z}^{3}}\varphi^{2}(y).$

From Lemma 4.2, we can deduce that for all $R>r+1$:

$S(r+1)-S(r-1)\geq\frac{S(r)(S(R)-S(r+1))}{3(R+r)}.$

Therefore, for $R\geq 16$ and $r\in[\frac{R}{4}+1,\frac{R}{2}-1]$ (so that
$r<R-2)$, by monotonicity of $S(r)$, one has

$\displaystyle
S(r+1)-S(r-1)\geq\frac{S(r)(S(R)-S(r+1))}{6R}\geq\frac{S\left(\frac{R}{4}\right)\left(S(R)-S\left(\frac{R}{2}\right)\right)}{6R}.$

Integrating $r$ from $\frac{R}{4}+1$ to $\frac{R}{2}-1$, we get

$\displaystyle\sum_{\frac{R}{4}+1}^{\frac{R}{2}-1}S(r+1)-\sum_{\frac{R}{4}+1}^{\frac{R}{2}-1}S(r-1)=~{}$
$\displaystyle
S\left(\frac{R}{2}\right)+S\left(\frac{R}{2}-1\right)-S\left(\frac{R}{4}+1\right)-S\left(\frac{R}{4}\right)$
$\displaystyle\leq~{}$ $\displaystyle
2S\left(\frac{R}{2}\right)-2S\left(\frac{R}{4}\right).$

Hence we obtain

$S\left(\frac{R}{2}\right)-S\left(\frac{R}{4}\right)\geq\frac{S\left(\frac{R}{4}\right)}{C_{1}}\left(S(R)-S\left(\frac{R}{2}\right)\right)$

for some universal constant $C_{1}$. The inequality holds for any $R\geq 16$.
In particular, for any $k\geq 0$, we have

$\begin{split}S\left(2^{k-1}R\right)-S\left(2^{k-2}R\right)&\geq\frac{S\left(2^{k-2}R\right)}{C_{1}}\left(S\left(2^{k}R\right)-S\left(2^{k-1}R\right)\right)\\\
&\geq\frac{S\left(\frac{R}{4}\right)}{C_{1}}\left(S\left(2^{k}R\right)-S\left(2^{k-1}R\right)\right).\\\
\end{split}$

Let us choose $C_{0}=2C_{1}$. Then, by taking $R=4R_{0}\geq 16$, we can state
that for all $k\geq 0$:

$\begin{split}S\left(2^{k}R\right)-S\left(\frac{R}{2}\right)&=\sum_{i=0}^{k}\left(S(2^{i}R)-S(2^{i-1}R)\right)\leq\sum_{i=0}^{k}\left(\frac{C_{1}}{S(\frac{R}{4})}\right)^{i+1}\left(S\left(\frac{R}{2}\right)-S\left(\frac{R}{4}\right)\right)\\\
&\leq S\left(\frac{R}{2}\right)-S\left(\frac{R}{4}\right).\end{split}$

In the last inequality, we utilized the fact that
$S\left(\frac{R}{4}\right)=S\left(R_{0}\right)\geq C_{0}=2C_{1}$. Then we get
that

$2S\left(2R_{0}\right)=2S\left(\frac{R}{2}\right)\geq
S\left(2^{k}R\right)=\sum_{y\in B_{2^{k}R}\cap\mathbb{Z}^{3}}\varphi^{2}(y)$

for any $k$. The lemma is proved by sending $k\to\infty$. ∎

Now we are ready to give the proof of Theorem 1.2.

###### Proof of Theorem 1.2..

Let

$R_{0}:=\inf_{R}\left\\{\,\exists x\in\mathbb{Z}^{3}\,\textrm{such
that}\sum_{y\in B_{R}(x)\cap\mathbb{Z}^{3}}\varphi^{2}(y)\geq C_{0}\right\\},$
(4.16)

where $C_{0}$ is the constant in Lemma 4.3. It is clear that
$R_{0}\in(2,\infty)$ for $m>C_{0}$. By Lemma 4.3, we obtain

$m\leq 2\sum_{y\in B_{2R_{0}}(x)\cap\mathbb{Z}^{3}}\varphi^{2}(y).$

For the ball $B_{2R_{0}}$, we can find $n$ points $x_{1},...,x_{l}\in
B_{2R_{0}}$ such that

$B_{2R_{0}}(x)\subset\bigcup_{i=1}^{l}B_{\left([R_{0}]-1\right)}(x_{i}),$

where $n$ is independent of $R_{0}$. It follows that

$m\leq 2\sum_{y\in B_{2R_{0}}(x)\cap\mathbb{Z}^{3}}\varphi^{2}(y)=2\sum_{y\in
B_{2R_{0}}(x)\cap\mathbb{Z}^{3}}\varphi^{2}(y)\leq 2\sum_{i=1}^{l}\sum_{y\in
B_{\left([R_{0}]-1\right)}(x_{i})\cap\mathbb{Z}^{3}}\varphi^{2}(y).$

One deduces from (4.16) that $m\leq 2nC_{0}<\infty$. ∎

##  5\. The nonexistence of a minimizer for (1.11)

In this section, we will prove Theorem 1.4 by contradiction. Suppose that
there exists a minimizer $\Omega$, we first demonstrate that the energy
satisfies the estimate:

$\inf\limits_{{|\Omega|=V}}E(\Omega)\lesssim V.$ (5.1)

On the other hand, we shall establish an estimate for the nonlocal term as
follows:

$\sum\limits_{x,y\in\Omega\hfill\atop y\neq x\hfill}{\frac{1}{|x-y|}}\gtrsim
V\log V.$ (5.2)

By combining (5.1) and (5.2), we arrive at a contradiction and it leads to the
original conclusion. To this end we define

$E(V):=\inf\\{E(\Omega)\mid|\Omega|=V\\}.$ (5.3)

The following lemma presents an inequality of subadditivity for $E(V)$.

######  Lemma 5.1.

Assuming $V_{0}>1$, $V_{1}>1$, and $V=V_{0}+V_{1}$, it follows that

$E(V)\leq E(V_{0})+E(V_{1}).$ (5.4)

###### Proof.

By the definition of infimum, for any given $\epsilon\geq 0$, there are
bounded sets ${\Omega}_{i},i=0,1$, such that

$E({\Omega}_{i})\leq E(V_{i})+\epsilon.$

Additionally, there exists a large ball $B_{R}$ such that ${\Omega}_{i}\subset
B_{R}(0)$. Let ${\widetilde{\Omega}}$ be defined as
${\Omega}_{0}+(d+{\Omega}_{1})$, where $d$ is a shift vector with $|d|>2R$. It
should be noted that when $x\in{\Omega}_{0}$ and $y\in d+{\Omega}_{1}$, we
have $|x-y|>d-2R>0$. It is evident that
$|\partial\widetilde{\Omega}|=|\partial{\Omega}_{0}|+|\partial{\Omega}_{1}|$
and $|\widetilde{\Omega}|=|{\Omega}_{0}|+|{\Omega}_{1}|$ due to the definition
of $\widetilde{\Omega}$. Thus, based on the definition of $E(V)$, we have:

$\begin{split}E(V)&\leq
E(\widetilde{\Omega})=|\partial{\Omega}_{0}|+|\partial{\Omega}_{1}|+\sum_{x,y\in{\Omega}_{0}\hfill\atop~{}\
y\neq x\hfill}{\frac{1}{|x-y|}}+\sum_{x,y\in d+{\Omega}_{1}\hfill\atop~{}\
~{}\ y\neq x\hfill}{\frac{1}{|x-y|}}+2\sum_{x\in{\Omega}_{0},~{}y\in
d+{\Omega}_{1}\hfill\atop~{}\ \quad y\neq x\hfill}{\frac{1}{|x-y|}}\\\ &\leq
E({\Omega}_{0})+E({\Omega}_{1})+\frac{2V_{0}V_{1}}{|d|-2R}\\\ &\leq
E(V_{0})+E(V_{1})+2\epsilon+\frac{2V_{0}V_{1}}{|d|-2R}.\end{split}$

As $|d|$ tends to infinity and $\epsilon$ tends to zero, the inequality is
established. ∎

The subsequent lemma establishes an estimate for $E(V)$ and elucidates the
properties of the minimizer $\Omega$.

######  Lemma 5.2.

Let $V\gg 1$, then

  1. $(i)$

We have $E(V)\lesssim V$.

  2. $(ii)$

If $\Omega$ is a minimizer of (1.11), then $\Omega$ must be a connected set.

###### Proof.

$(i).$ By Lemma 2.1,

$|\partial\Omega|\sim V^{\frac{2}{3}}~{}\mbox{for a ball}~{}\Omega~{}\mbox{of
volume}~{}V.$

Using the fact that $|x-y|\geq 1$ we have

$\sum\limits_{x,y\in{\Omega}\hfill\atop~{}y\neq
x\hfill}{\frac{1}{|x-y|}}\lesssim V^{2}.$

Therefore

$E(\Omega)\lesssim V^{2}+V^{\frac{2}{3}},$

which implies

$E(V)\lesssim 1\quad\mbox{for}\quad 1\leq|V|<2.$ (5.5)

Through the process of induction and utilizing the subadditivity inequality
(5.4), we are able to derive that

$E(V)\lesssim N\quad\textrm{for}\quad N\leq V<(N+1).$

$(ii).$ When $\Omega$ is not a connected domain, we can divide it into two
separate regions $\Omega_{1}$ and $\Omega_{2}$ in such a way that no point
$x\in\Omega_{1}$ is connected to any point $y\in\Omega_{2}$ (i.e., $x\nsim
y$). Let us define

$\widehat{\Omega}={\Omega}_{1}+(d+{\Omega}_{2}),$

where $d$ is a shift vector with $|d|>1$. It is clear that

$|\partial\widehat{\Omega}|=|\partial{\Omega}_{1}|+|\partial{\Omega}_{2}|,\quad|\widehat{\Omega}|=|{\Omega}_{1}|+|{\Omega}_{2}|=V,$

and

$\sum\limits_{x,y\in{\widehat{\Omega}}\hfill\atop~{}y\neq
x\hfill}\frac{1}{|x-y|}<\sum\limits_{x,y\in{\Omega}\hfill\atop~{}y\neq
x\hfill}\frac{1}{|x-y|}.$

Therefore, we have $E(\widehat{\Omega})<E(\Omega)$, which contradicts the
definition of $\Omega$. ∎

The following lemma provides lower bound for the nonlocal expression
$\sum\limits_{x,y\in{\Omega}\hfill\atop~{}y\neq x\hfill}\frac{1}{|x-y|}$.

######  Lemma 5.3.

Suppose that $V\gg 1$ and $\Omega$ minimizes (1.11), then the following
estimate holds for the nonlocal expression:

$\sum\limits_{x,y\in{\Omega}\hfill\atop~{}y\neq x\hfill}\frac{1}{|x-y|}\gtrsim
V\log V.$ (5.6)

###### Proof.

Our first assertion is that

$|\Omega\cap B_{R}(x)|\gtrsim R\,\,\,\,\textrm{for all
vertics}\,\,x\in\Omega\,\,\textrm{and
radii}\,1<R<\frac{1}{2}\textrm{diam}(\Omega).$ (5.7)

Let $x$ be a fixed point in $\Omega$ and
$1<R<\frac{1}{2}\textrm{diam}(\Omega)$ be a chosen radius. By the definition
of the diameter of $\Omega$, we see that $\Omega$ is not entirely contained
within $B_{R}(x)$. Using the fact that $\Omega$ is connected, we are able to
find a path $\\{x_{i}\\}_{i=1}^{R}\in\Omega\cap B_{R}(x)$ such that $x=x_{1}$.
This implies that $|\Omega\cap B_{R}(x)|\geq R$ and it proves the assertion
(5.7). As depicted in Figure (3).

Figure 3. The sketch of $B_{R}$ and $\Omega$.

Now we claim

$\sum\limits_{x,y\in{\Omega}\hfill\atop~{}y\neq x\hfill}\frac{1}{|x-y|}\gtrsim
V\log\textrm{diam}(\Omega).$

It is worth noting that, for all $1<t<\frac{1}{2}\textrm{diam}(\Omega)$,

$\begin{split}A(t):=\big{|}\\{(x,y)\in\Omega\times\Omega||x-y|<t\\}\big{|}=\sum_{x\in\Omega}\big{|}\\{y\in\Omega||y-x|<t\\}\big{|}=\sum_{x\in\Omega}\big{|}\Omega\cap
B_{t}(x)\big{|}\gtrsim t|\Omega|.\end{split}$

Hence we get from Lemma 2.1 ($\textrm{diam}\,\Omega\gtrsim V^{\frac{1}{3}}$)
that

$\displaystyle\sum\limits_{x,y\in{\Omega}\hfill\atop~{}y\neq
x\hfill}\frac{1}{|x-y|}$
$\displaystyle\geq\sum_{t=1}^{\infty}\frac{1}{t}\big{|}\\{(x,y)\in\Omega\times\Omega|t-1\leq|x-y|<t\\}\big{|}$
$\displaystyle\geq\sum_{t=1}^{\infty}\frac{1}{t}(A(t)-A(t-1))\geq\sum_{t=1}^{\frac{1}{2}\textrm{diam}(\Omega)}\frac{1}{t(t+1)}A(t)$
$\displaystyle\geq\sum_{t=1}^{\frac{1}{2}\textrm{diam}(\Omega)}\frac{1}{2t}|\Omega|\gtrsim|\Omega|\log|\textrm{diam}(\Omega)|\gtrsim
V\log V.$

Therefore, we obtain the estimate (5.6). ∎

###### Proof of Theorem 1.4.

Theorem 1.4 follows from Lemma 5.2 and Lemma 5.3. ∎

## References

  * [1] G. Alberti, R. Choksi, F. Otto,  Uniform energy distribution for an isoperimetric problem with long-range interactions.  J. Amer. Math. Soc., 22 (2009), 569-605.
  * [2] R. Benguria, H. Brézis, E. H. Lieb,  Thomas-Fermi-von Weizsäcker theory of atoms and molecules.  Comm. Math. Phys., 79 (1981), 167-180.
  * [3] C.-M. Care, N. H. March,  Electron crystallization.  Adv. Phys., 24 (1975), 101-116.
  * [4] L.-Q. Chen, A. G. Khachaturyan,  Dynamics of simultaneous ordering and phase separation and effect of long-range Coulomb interactions.  Phys. Rev. Lett., 70 (1993), 1477-1480.
  * [5] R. Choksi, M.-A. Peletier,  Small volume fraction limit of the diblock copolymer problem: I. Sharp-interface functional.  SIAM J. Math. Anal., 42 (2010), 1334-1370.
  * [6] R. Choksi, M. A.Peletier,  Small volume-fraction limit of the diblock copolymer problem: II. Diffuse-interface functional.  SIAM J. Math. Anal., 43 (2011), 739-763.
  * [7] M. Cicalese, E. Spadaro,  Droplet minimizers of an isoperimetric problem with long-range interactions.  Comm. Pure Appl. Math., 66 (2013), 1298-1333.
  * [8] F. R. K. Chung and S.-T. Yau,  Eigenvalues of graphs and Sobolev inequalities.  Combin. Probab. Comput., 4 (1995), 11-25.
  * [9] R.-L. Frank, E.-H. Lieb,  A compactness lemma and its application to the existence of minimizers for the liquid drop model.  SIAM J. Math. Anal., 47 (2015), 4436-4450.
  * [10] C.-L. Fefferman, L. A. Seco,  Asymptotic neutrality of large ions.  Comm. Math. Phys., 128 (1990), 109-130.
  * [11] E. Fermi,  Un metodo statistico per la determinazione di alcune prioprietà dell’atomo.  Rend. Accad. Naz. Lincei., 32 (1927), 602-607.
  * [12] R.-L. Frank, P.-T. Nam , H. Van Den Bosch,  The ionization conjecture in Thomas-Fermi-Dirac-von Weizsäcker theory.  Comm. Pure Appl. Math., 71 (2018), 577-614.
  * [13] G. Gamow,  Mass defect curve and nuclear constitution.  Proc. R. Soc. Lond. A, 126 (1930), 632-644.
  * [14] A. Grigor’yan,  Introduction to analysis on graphs, University Lecture Series, 71. American Mathematical Society, Providence, RI, 2018.
  * [15] A. Grigor’yan, Y. Lin, Y.-Y. Yang,  Yamabe type equations on graphs.  J. Differ. Equ., 261 (2016), 4924-4943.
  * [16] A. Grigor’yan, Y. Lin, Y.-Y. Yang,  Kazdan-Warner equation on graph.  Calc. Var. Partial Differ. Equ., 55 (2016), 1-13.
  * [17] A. Grigor’yan, Y. Lin, Y.-Y. Yang,  Existence of positive solutions to some nonlinear equations on locally finite graphs.  Sci. China Math., 60 (2017), 1311-1324.
  * [18] C. Gong, Y. Lin, S. Liu, S.-T. Yau,  Li-Yau inequality for unbounded Laplacian on graphs.  Adv. Math., 357 (2019), 106822, 23 pp.
  * [19] P. Horn, Y. Lin,Liu, S. Liu, S.-T. Yau,  Volume doubling, Poincaré inequality and Gaussian heat kernel estimate for non-negatively curved graphs.  J. Reine Angew. Math., 757 (2019), 89-130.
  * [20] A. Huang, Y. Lin, S.-T. Yau,  Existence of solutions to mean field equations on graphs.  Commun. Math. Phys., 377 (2020) 613-621.
  * [21] H.-Y. Huang, J. Wang, W. Yang,  Mean field equation and relativistic Abelian Chern-Simons model on finite graphs.  J. Funct. Anal., 281 (2021), 109218.
  * [22] B. B. Hua, D. Mugnolo,  Time regularity and long-time behavior of parabolic p-Laplace equations on infinite graphs.  J. Differ. Equ., 259 (2015), 6162-6190.
  * [23] B. Hua, R. Li,  The existence of extremal functions for discrete Sobolev inequalities on lattice graphs.  J. Differ. Equ., 305 (2021), 224-241.
  * [24] B. Hua, W. D. Xu,  Existence of ground state solutions to some Nonlinear Schrödinger equations on lattice graphs.  Calc. Var. Partial Differential Equations, 62 (2023), 127\.
  * [25] G.-G. Huang, C.-M. Li, X. Yin,  Existence of the maximizing pair for the discrete Hardy-Littlewood-Sobolev inequality,  Discrete Contin. Dyn. Syst., 35 (2015), 935-942.
  * [26] M. Keller, M. Schwarz,  The Kazdan-Warner equation on canonically compactifiable graphs.  Calc. Var. Partial Differ. Equ., 57 (2018), 1-18.
  * [27] H. Knüpfer, C. B. Muratov,  On an isoperimetric problem with a competing nonlocal term I: the planar case.  Comm. Pure Appl. Math., 66 (2013), 1129-1162.
  * [28] H. Knüpfer, C. B. Muratov,  On an isoperimetric problem with a competing nonlocal term II: the general case.  Comm. Pure Appl. Math., 67 (2014), 1974-1994.
  * [29] W.-J. Liu, K.-W. Chen, J. Yu,  Extinction and asymptotic behavior of solutions for the $\omega$-heat equation on graphs with source and interior absorption,  J. Math. Anal. Appl., 435 (2016), 112-132.
  * [30] E.-H, Lieb,  Thomas-Fermi and related theories of atoms and molecules.  Rev. Modern Phys., 53 (1981), 603-641.
  * [31] E.-H. Lieb,  Bound on the maximum negative ionization of atoms and molecules.  Phys. Rev. A., 29 (1984), 3018-30.
  * [32] E. H. Lieb, I. M. Sigal, B. Simon, W. Thirring,  Asymptotic neutrality of large-$Z$ ions.  Phys. Rev. Lett., 52 (1984), 994-996.
  * [33] E.-H. Lieb, I. M. Sigal, B. Simon, W. Thirring,  Approximate neutrality of large-Zions.  Comm. Math. Phys., 116 (1988), 635-644.
  * [34] E.-H. Lieb, R. Seiringer,  The stability of matter in quantum mechanics, Cambridge University Press, Cambridge, 2010.
  * [35] P.-L. Lions,  Solutions of Hartree-Fock equations for Coulomb systems.  Comm. Math. Phys., 109 (1987), 33-97.
  * [36] C. Le Bris,  Some results on the Thomas-Fermi-Dirac-von Weizsäcker model.  Differ. Integral Equa., 6 (1993), 337-353.
  * [37] F.-J. Lu, F. Otto,  Nonexistence of a minimizer for Thomas-Fermi-Dirac-von Weizsäcker model.  Comm. Pure Appl. Math., 67 (2014), 1605-1617.
  * [38] L. Modica,  The gradient theory of phase transitions and the minimal interface criterion.  Arch. Rational Mech. Anal., 98 (1987), 123-142.
  * [39] C.-B. Muratov, Theory of domain patterns in systems with long-range interactions of Coulomb type, Phys. Rev. E., 66 (2002), 066108.
  * [40] C.-B. Muratov,  Droplet phases in non-local Ginzburg-Landau models with Coulomb repulsion in two dimensions,  Comm. Math. Phys., 299 (2010), 45-87.
  * [41] A. Nagai, Y. Kametaka, K. Watanabe,  The best constant of discrete Sobolev inequality.  J. Phys. A., 42 (2009), 454014.
  * [42] P.-T. Nam,  New bounds on the maximum ionization of atoms.  Comm. Math. Phys., 312 (2012), 427-445.
  * [43] P.-T. Nam,  The ionization problem.  Eur. Math. Soc. Newsl. No. 118 (2020), 22-27.
  * [44] P.-T. Nam, H. Van Den Bosch,  Nonexistence in Thomas-Fermi-Dirac-von Weizsäcker theory with small nuclear charges.  Math. Phys. Anal. Geom., 20 (2017), Art. 6, 32 pp.
  * [45] M.-I. Ostrovskii,  Sobolev spaces on graphs.  Quaest. Math., 28 (2005) 501-523.
  * [46] T. Ohta, K. Kawasaki,  Equilibrium morphology of block copolymer melts.  Macromolecules, 19 (2009), 2621-2632.
  * [47] D. Ruiz,  On the Schrödinger-Poisson-Slater system: behavior of minimizers, radial and nonradial cases.  Arch. Ration. Mech. Anal., 198 (2010), 349-368.
  * [48] O. Sánchez, J. Soler,  Long-time dynamics of the Schrödinger-Poisson-Slater system.  J. Statist. Phys., 114 (2004), 179-204.
  * [49] B. Simon,  Schrödinger operators in the twenty-first century, Mathematical physics 2000, 283-288, Imp. Coll. Press, London, 2000.
  * [50] L.-A. Seco, I.-M. Sigal, and J.-P. Solovej,  Bound on the ionization energy of large atoms.  Commun. Math. Phys., 131 (1990), 307-315.
  * [51] J.-P. Solovej,  Proof of the ionization conjecture in a reduced Hartree-Fock model.  Invent. Math., 104 (1991), 291-311.
  * [52] J.-P. Solovej,  The ionization conjecture in Hartree-Fock theory.  Ann. of Math., 158 (2003), 509-576.
  * [53] L.-L. Sun, L.-Q. Wang,  Brouwer degree for Kazdan-Warner equations on a connected finite graph.  Adv. Math., 404 (2022), part B, Paper No. 108422, 29 pp.
  * [54] L.-H. Thomas,  The calculation of atomic fields.  Proc. Cambridge Phil. Soc., 23 (1927), 542-548.
  * [55] J. Wang, Y.-N. Zhu, K. Wang,  Existence and asymptotical behavior of the ground state solution for the Choquard equation on lattice graphs.  Electron, Res. Arch., 31 (2023), 812-839.
  * [56] C.-F. Weizsäcker,  Zur Theorie der Kernmassen,  Z. Phys., 96 (1935), 431–458.
  * [57] R.-K. Wojciechowski,  Heat kernel and essential spectrum of infinite graphs.  Indiana Univ. Math. J., 58 (2009), 1419–1441.


# Identification of Probability weighted ARX models with arbitrary domains

Alessandro Brusaferria,b, Matteo Matteuccib, Stefano Spinellia,b aCNR-
Institute of Intelligent Industrial Technologies and Systems for Advanced
Manufacturing, Milan, Italy
bPolitecnico di Milano - Department of Electronics, Informatics and
Bioengineering, Milan, Italy
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Hybrid system identification is a key tool to achieve reliable models of
Cyber-Physical Systems from data. PieceWise Affine models guarantees universal
approximation, local linearity and equivalence to other classes of hybrid
system. Still, PWA identification is a challenging problem, requiring the
concurrent solution of regression and classification tasks. In this work, we
focus on the identification of PieceWise Auto Regressive with eXogenous input
models with arbitrary regions (NPWARX), thus not restricted to polyhedral
domains, and characterized by discontinuous maps. To this end, we propose a
method based on a probabilistic mixture model, where the discrete state is
represented through a multinomial distribution conditioned by the input
regressors. The architecture is conceived following the Mixture of Expert
concept, developed within the machine learning field. To achieve nonlinear
partitioning, we parametrize the discriminant function using a neural network.
Then, the parameters of both the ARX submodels and the classifier are
concurrently estimated by maximizing the likelihood of the overall model using
Expectation Maximization. The proposed method is demonstrated on a nonlinear
piece-wise problem with discontinuous maps.

###### Index Terms:

Hybrid systems identification, NPWARX, Neural Network, Probabilistic mixture,
Expectation Maximization

##  I Introduction

The identification of hybrid dynamics is fundamental to achieve reliable
models of Cyber-Physical Systems (CPS) from data [1]. Hybrid systems provide a
unified framework to represent the heterogeneous interactions between control
logic/rules and the continuous dynamics of CPS, including sharp changes in
operating points, plants regimes, constraints on values/variations of system
inputs/outputs, etc. Several models and identification approaches have been
proposed in the control community during the last decades (see e.g, [2],[3]
for detailed reviews). A major distinction is in the characterization of the
discrete state evolution, which is often defined as autonomous or due to
partitioning of the continuous state-input domain, leading to switching and
piecewise models.

In this work, we focus on the latter, representing an important class -
subject of a continuous interest - due to its universal approximation
properties, local linearity and equivalence to other hybrid system modeling
approaches, thus fostering analysis and deployment in advanced control
applications [4].

Considering input-output models, PWARX (i.e., PieceWise Auto Regressive
systems with eXogenous input) are structured by a set of affine maps and
through a partition of the continuous feature space, determining the discrete
mode [2]. The identification of PWARX models is a challenging problem,
requiring the concurrent estimation of the sub-model dynamics and the
partitioning function through intertwined regression and classification tasks
[3]. Besides, as logic conditions in hybrid systems can cause sudden changes
in the dynamics, non-smooth function approximations have to be addresses
(e.g., using discontinuous PWARX mappings), which usually results challenging
using classical nonlinear regression techniques assuming continuous patterns
[5].

Mixed-integer programming approaches have been proposed to achieve global
optimal solutions to PWARX identification (see e.g., [6] and references
therein). As the problem is NP-hard, these approaches result affordable only
for very small data sets [6]. Hence, a lot of research studies have been
dedicated to approximate methods and heuristics aimed to find good sub-
optional solution in reasonable time. In this context, the most popular
methods include the bounded error approach [7], the Bayesian approach [8],
sparse optimization [9] and the broad family of clustering-based approaches
(see e.g., [10],[11],[12]). A detailed review is reported in [2]. Despite the
different techniques employed, the common workflow starts from local model
estimation and clustering by relaxing the piecewise constraints, followed by
partitioning through linear discrimination techniques [3]. Hence, as the
obtained clusters are not guaranteed to be linearly separable, a post-
processing phase is usually required to address misclassifications [6].

A probabilistic generalization of PWARX (i.e., PrARX) is proposed in [13],
resulting in a Mixture of Expert (MoE) architecture [14], enabling concurrent
estimation of the ARXs and parameters of the boundaries in a single
optimization problem. Then, the corresponding PWARX model is obtained using a
straighforward procedure. A variational inference based technique is
introduced [15], supporting model structure selections in cases of small data-
set size.

Nevertheless, the vast majority of existing studies focus on piecewise models
characterized by polyhedral domains in the regression space. The introduction
of arbitrary regions leads to the nonlinearly piecewise models (i.e, NPWARX),
a direct extension of PWARX still not deeply studied and solved in the
literature [16]. As a matter of fact, when the target system is characterized
by these arbitrary partitions, forcing linear separability through PWARX
models leads to multiple identical submodels in different regions governed by
the same dynamics, while NPWARX models enable a unique representation.
Besides, NPWARX supports better estimations of submodels, that would be
assigned by PWARX to different modes, particularly in sub-regions with fewer
samples [16]. However, NPWARX introduce a more challenging nonlinear
classification problem within the identification procedure. To address this
issue, authors in [16] propose a technique based on kernel regression and
Support Vector Machines.

In this work, we investigate a different approach, leveraging on tools from
the probabilistic mixture field. Specifically, we extend the PrARX model by
introducing a nonlinear partitioning function parametrized by a neural
network. Estimation is performed through an integrated optimization, by
maximizing the overall model likelihood. Then, we tackle this challenging
problem through Expectation Maximization (EM), disentangling classification
and regression sub-tasks easier to solve within each iteration.

The paper is organized as follows: Section 2 states the identification problem
and details the PrARX model; Section 3 reports the developed model and EM-
based estimation procedure; Section 4 describes the adopted case study and
summarizes the results achieved.

##  II Background

###  II-A PWARX model identification

Consider discrete-time dynamical systems of the form:

$y_{k}=f(x_{k})+e_{k}$ (1)

where $y_{k},e_{k}\in\mathbb{R}$ are, respectively, the observed system output
and error term at discrete time $k\in\mathbb{Z}$ while $f(.)$ is a
deterministic mapping from the regression vector:

$x_{k}=[y_{k-1},...,y_{k-n_{a}},u^{T}_{k-1},...,u^{T}_{k-n_{b}}]^{T}$ (2)

$u_{k}\in\mathbb{R}^{q}$ represent exogenous inputs and
$n_{a},n_{b}\in\mathbb{N}^{+}$ fixed input and output lags (i.e., model
orders).

PieceWise affine AutoRegressive eXogenous (PWARX) models are defined by a
piecewise affine mapping as follows:

$\displaystyle\begin{dcases}y_{k}&=f_{1}(\varphi_{k})=\theta_{1}^{T}\varphi_{k}+e_{k},\quad\text{if}\quad
x_{k}\in\mathcal{X}_{1}\\\ &...\\\
y_{k}&=f_{s}(\varphi_{k})=\theta_{s}^{T}\varphi_{k}+e_{k},\quad\text{if}\quad
x_{k}\in\mathcal{X}_{s}\\\ \end{dcases}$ (3)

where $\varphi_{k}\equiv[x_{k}\ 1]^{T}$ is the extended regression vector,
$S\in\mathbb{N}^{+}$ is the finite number of discrete state (or modes), with
affine submodels parameters $\theta_{s}\in\mathbb{R}^{r}$ for
$s\in\left\\{1,...,S\right\\}$, with $r=n_{a}+q\cdot n_{b}+1$. The switching
mechanism between the discrete states is determined by a complete partition of
the regressor domain $\mathcal{X}\subseteq\mathbb{R}^{r}$ in the collection of
convex polyhedra (i.e, regions) $\left\\{\mathcal{X}_{s}\right\\}_{s=1}^{S}$
described by:

$\mathcal{X}_{s}=\left\\{x_{k}\in\mathbb{R}^{r}:H_{s}\varphi_{k}\preceq_{[s]}0\right\\}$
(4)

where the matrices $H_{s}\in\mathbb{R}^{\mu_{s}\times r}$ represent a set of
hyperplanes defining the regions by $\mu_{s}$ linear inequalities and
’$\preceq_{[s]}$’ denotes a $\mu_{s}$-dimensional vector whose elements can be
the symbols $\leq$ and $<$. Hence, the active discrete state $z_{k}$ is given
by the region to which the regressor $x_{k}$ belongs at $k$:

$z_{k}=s\Leftrightarrow x_{k}\in\mathcal{X}_{s},\ s=1,...,S$ (5)

The general PWARX identification problem include the estimation of the ARX
parameters $\left\\{\theta_{s}\right\\}_{s=1}^{S}$ and the regions
$\left\\{\mathcal{X}_{s}\right\\}_{s=1}^{S}$, as well as model orders
$n_{a},n_{b}$ and the number of modes $S$, given a collection of input/output
data $\mathcal{D}=\left\\{(x_{k},y_{k})\right\\}_{k=1}^{N}$.
In this work, we adopt the common assumption of fixed modes/orders, often
estimated in practice by preliminary data analysis, cross validation and model
order selection techniques [3]. The investigation of automatic feature
selection techniques (e.g., including dedicated regularization terms) is left
to future extensions. Besides, we assume that in $\mathcal{D}$ all modes are
sufficiently excited and inputs and outputs are bounded [17].

###  II-B Nonlinearly PieceWise ARX model

As introduced above, NPWARX models extend PWARX by considering nonlinear
boundaries between arbitrary regions, thus not restricted to set of
hyperplanes defining polyhedral domains [16]. Formally, NPWARX are defined as
in (3) while relaxing the polyhedrical assumption in (4).

###  II-C Probability weighted ARX model

In the PrARX model proposed in [13], the deterministic partition of PWARX is
replaced by probabilistic boundaries, obtained by a softmax function over a
linear transformation of the extended regressor (i.e, linear gate):

$\displaystyle y_{k}$ $\displaystyle=f_{Pr}(x_{k})+e_{k},\ \ \
f_{Pr}(x_{k})=\sum_{s=1}^{S}p_{s}\theta_{s}^{T}\varphi_{k}$ (6) $\displaystyle
p_{s}$
$\displaystyle=\frac{\exp\left(\eta_{s}^{T}\varphi_{k}\right)}{1+\sum_{j=1}^{S-1}\exp\left(\eta_{j}^{T}\varphi_{k}\right)}$
(7)

where $\eta_{j}\in\mathbb{R}^{r}$ and $p_{s}$ denotes the probability that
$\varphi_{k}$ belongs to the discrete state $s$.
Parameters $\eta$ are directly related to the hyperplanes determining the
regions. Indeed, the PrARX model can be straightforwardly transformed to a
PWARX model by:

$H_{s}=\left[(\eta_{1}-\eta_{s})...(\eta_{S}-\eta_{s})\right]^{T},\ s=1,...,S$
(8)

Besides being a probabilistic generalization of PWARX, the PrARX represents a
particular form of the MoE architecture, including linear experts and linear
gates [14],[18]. Due to this simplified parametrization, the steepest descent
method is adopted in [13] to maximize the likelihood during estimation.

##  III Methods

###  III-A Nonlinearly PrARX model

We augment the PrARX model, following the Mixture of Expert approach, by
replacing the linear gate with a neural network aimed to learn arbitrary
region partitioning the regressor space into discrete states. Hence, we obtain
a probabilistic generalization of NPWARX models, as PrARX are for PWARX [15].
We briefly denote it NPrARX hereinafter.

Specifically, the discrete state $z_{k}$, at discrete time $k\in\mathbb{Z}$,
is represented by a latent categorical random variable taking values in the
set $\left\\{1,...,S\right\\}$, generated according to a $x_{k}$-conditioned
multinomial distribution:

$z_{k}|x_{k}\sim Mult(1;p(z_{k}=1|x_{k}),...,p(z_{k}=S|x_{k}))$ (9)

The output variable $y_{k}$ is generated according to a set of emission
distributions, conditioned by $x_{k}$, selected by the active discrete mode:

$y_{k}|z_{k}=s,x_{k}\ \sim\ p(y_{k}|x_{k},z_{k}=s)$ (10)

Hence, data generation follows a hierarchical process, starting from the mode
sampling followed by output emission, both conditioned on the regressor.

Considering an ARX model form in each discrete state, and assuming a Gaussian
distribution with zero mean and standard deviation $\sigma_{s}$ for $e_{k}$,
we obtain a set of $s$ noisy linear models:

$y_{k}|z_{k}=s,x_{k}\ \sim
N(.;\theta_{s,0}+\theta_{s,1:}^{T}x_{k},\sigma_{s}^{2})$ (11)

It is worth noting that, despite the Gaussian noise assumption adopted in this
work, alternative noise patterns can be considered, including specific
distributions for each mode [19].
The overall model results in a non-Markov switching probabilistic mixture
[15], with semi-parametric density defined as:

$\displaystyle\begin{aligned}
p(y_{k}|x_{k},\Theta)&=\sum_{s=1}^{S}p(z_{k}=s,y_{k}|x_{k})\\\
&=\sum_{s=1}^{S}p(z_{k}=s|x_{k})p(y_{k}|x_{k},z_{k}=s)\\\
&=\sum_{s=1}^{S}g_{s}(x_{k},\mathcal{W})p(y_{k}|x_{k},\theta_{s},\sigma_{s})\end{aligned}$
(12)

where
$\Theta=\left\\{\mathcal{W},\left\\{\mathcal{\theta}_{s},\sigma_{s}\right\\}_{s=1}^{S}\right\\}$
summarizes the parameters to be estimated and
$p(y_{k}|x_{k},\theta_{s},\sigma_{s})$ is the $s$-th mode density:

$p(y_{k}|x_{k},\theta_{s},\sigma_{s})=\frac{1}{\sqrt{2\pi}\sigma_{s}}e^{-\frac{1}{2\sigma_{s}^{2}}\left(y_{k}-\theta_{s}^{T}\varphi_{k}\right)^{2}}$
(13)

We structure the function $g_{s}(x_{k},\mathcal{W})$, that provides the
probability of each discrete state given the regressor, as follows:

$g_{s}(x_{k},\mathcal{W})=\frac{\exp\left(\alpha_{s}(x_{k},\mathcal{W})\right)}{1+\sum_{j=1}^{S-1}\exp\left(\alpha_{j}(x_{k},\mathcal{W})\right)}$
(14)

where the logits $\alpha_{j}(x_{k},\mathcal{W})$ are defined using a neural
network processing the regressor vector $x_{k}$. In the probabilistic mixture
literature, this is often referred to as gate. Besides, the parameters of one
of the $S$ modes is forced to the null vector to support identifiability, as
suggested in [14].
While several network architectures can be adopted, depending on the specific
characteristics of the application at hand, in this work we employ a feed-
forward form, defined as:

$\alpha_{s}(x_{k},\mathcal{W})=\sum\limits_{l=1}^{n_{h}}W_{l,s}^{(2)}f_{l}\Bigg{[}\sum_{i=1}^{r-1}W_{i,l}^{(1)}x_{n}(i)+W_{0,l}^{(1)}\Bigg{]}+W_{0,s}^{(2)}$
$\small f_{l}(x)=tanh(x)=\frac{e^{x}-e^{-x}}{e^{x}-e^{-x}}$ (15)

where we reported a single hidden layer of $n_{h}\in\mathbb{N}^{+}$ units to
lighten notation. It it worth noting that alternative nonlinear activation
functions can be adopted and that the number of hidden units and layers
constitute hyperparameters to be configured. To this end, in this work, we
exploit cross-validation as detailed in section IV.

As introduced above, conditioning the gating mechanism through flexible models
such as neural networks provides trainable nonlinear decision boundaries for
discrete state classification in the regressor space, as opposed to the
polyhedral partitioning of PrARX. On the other hand, such flexibility comes at
the cost of a more complex optimization problem for parameters estimation
[19]. To address this issues, we employ Expectation Maximization (EM), as
detailed in the next section.

###  III-B Parameters estimation by Expectation Maximization

A major feature of the Probabilistic ARX approach is the concurrent learning
of classification (i.e., gating function) and regression (i.e., ARX submodels)
tasks, while pursuing parameters estimation.

To this end, two major families of approaches can be followed, namely the
frequentist and the Bayesian. In this work we adopt the former, leaving the
investigation of the second one to future extension, e.g., through the
exploitation of Bayesian Neural Networks. Hence, we target a Maximum
Likelihood Estimation (MLE) of model’s parameters.

Specifically, given a data set of i.i.d. observations
$\mathcal{D}=\left\\{(x_{k},y_{k})\right\\}_{k=1}^{N}$, the likelihood
function factorizes as follows:

$\mathcal{L}(\Theta)=\prod_{k=1}^{N}\left[\sum_{s=1}^{S}g_{s}(x_{k},\mathcal{W})p(y_{k}|x_{k},\theta_{s},\sigma_{s})\right]$
(16)

leading to the following objective:

$\hat{\Theta}\small{=}\underset{\Theta}{\small{argmax}}\left[\sum_{k=1}^{N}ln\left[\sum_{s=1}^{S}g_{s}(x_{k},\mathcal{W})p(y_{k}|x_{k},\theta_{s},\sigma_{s})\right]\right]$
(17)

The maximization of this function is a complex task, due to the inclusion of
the hidden discrete state. This is a common issue of latent variable models,
from probabilistic mixtures to Hidden Markov Models. In this fields, EM is
typically adopted to derive the MLE. In fact, it provides a more effective
alternative to the direct likelihood maximization (e.g, though steepest
descent, Newton-Raphson, etc.), specially when dealing with complex models and
conditioning functions [19].

Representing a general purpose algorithm for estimation with missing data, EM
has several attractive features. First of all, as opposed to the other
methods, it has been proved that the observed log-likelihood increase at each
iteration [19]. Moreover, it provides a natural way to tackle the challenging
sum inside the logarithm form of (17), and enables the exploitation of
effective closed form solution in the maximization sub-problems during the
iterations.

Considering a set of random indicator variables assigning a generating mode
$s$ to each sample $k$ in the observation set, $Z=\left\\{z_{k,s};\
k=1,...,N,s=1,...,S\right\\}$, we get the joint likelihood over the complete
data (i.e, including the missing indicator variables), factorizing as follows:

$\mathcal{L}_{EM}(\Theta)=\prod_{k=1}^{N}\prod_{s=1}^{S}\left[g_{s}(x_{k},\mathcal{W})p(y_{k}|x_{k},\theta_{s},\sigma_{s})\right]^{z_{k,s}}$
(18)

thus leading to a much more tractable form:

$\small
ln\mathcal{L}_{EM}=\sum_{k=1}^{N}\sum_{s=1}^{S}z_{k,s}\left[ln(g_{s}(x_{k},\mathcal{W}))+ln(p(y_{k}|x_{k},\theta_{s},\sigma_{s}))\right]$
(19)

It is worth noting that the indicator variables $z_{k,s}$ filter out all but
the related discrete state in the overall model.
Then, the missing indicator variables are averaged out by iteratively
computing the expectation of the complete data likelihood over them
$\mathbb{E}_{Z}[ln(p(X,Y,Z|\Theta))]$

$\displaystyle\begin{aligned}
\mathcal{Q}(\Theta^{i},\Theta^{i-1})&=\sum_{k=1}^{N}\sum_{s=1}^{S}\xi_{k,s}^{\Theta^{i-1}}ln[g_{s}(x_{k},\mathcal{W}^{i})]\\\
&+\xi_{k,s}^{\Theta^{i-1}}ln[p_{k}(y_{k}|x_{k},\theta^{i}_{s},\sigma_{s}^{i})]\\\
&=\mathcal{Q}^{c}(\mathcal{W}^{i})+\sum_{s=1}^{K}\mathcal{Q}_{s}^{r}(\theta^{i}_{s},\sigma_{s})\end{aligned}$
(20)

where $i=1,..,M$ represents the number of iteration of the EM algorithm.
Notably, the expectation over the binary indicator variables $z_{k,s}$ assume
continuous values, thus simplifying the optimization problem.
$\xi_{k,s}^{\Theta^{i-1}}$ represents the expected value of $z_{k,s}$, defined
as follows:

$\displaystyle\begin{aligned}
\xi_{k,s}^{\Theta^{i-1}}&=\mathbb{E}[z_{k,s}|X,Y,\Theta^{i-1}]=p(z_{k}=s|x_{k},y_{k},\Theta^{i-1})\\\
&=\frac{p(z_{k}=s,y_{k}|x_{k},\Theta^{i-1})p(y_{k}|x_{k},\theta^{i-1}_{s},\sigma_{s}^{i-1})}{p(y_{k}|x_{k})}\\\
&=\frac{g_{s}(x_{k},\mathcal{W}^{i-1})p(y_{k}|x_{k},\theta^{i-1}_{s},\sigma_{s}^{i-1})}{\sum_{j=1}^{S}g_{j}(x_{k},\mathcal{W}^{i-1})p(y_{k}|x_{k},\theta^{i-1}_{j},\sigma_{j}^{i-1})}\end{aligned}$
(21)

constituting the posterior obtained by the normalized product of the mode-wise
likelihood times the prior $g_{s}(x_{k},\mathcal{W})$.

The training algorithm proceed by alternating between posterior
$\xi_{k,s}^{\Theta^{i-1}}$estimation (i.e, E-step) and maximization of the
$\mathcal{Q}$-function (i.e, M-step), devoted to update the estimation of
$\Theta$. The former employs the parameters computed by the previous M-step,
starting from an initial guess.

Notably, this approach disentangles the estimation of regression (ARX) and
classification (gate) parameters within the iteration, resulting in sub-
problems easier to tackle. Specifically, the former results in a set of $S$
weighted least squared problems
$\mathcal{Q}_{s}^{r}(\theta^{i}_{s},\sigma_{s})$, easily solvable by:

$\theta_{s}^{i}=\underset{\theta_{s}}{argmin}\sum_{k=1}^{N}\xi_{k,s}^{\Theta^{i-1}}\left\\{-\frac{1}{2(\sigma_{s}^{i-1})^{2}}\left(y_{k}-\theta_{s}^{T}\varphi_{k}\right)\right\\}$
(22) $\theta_{s}^{i}=(\Phi^{T}\Xi_{k}\Phi)^{-1}\Phi^{T}\Xi_{s}y\ ,\
\Xi_{s}=diag(\xi_{k,:}^{\Theta^{i-1}})$ (23)

The latter, i.e., $\mathcal{Q}^{c}(\mathcal{W}^{i})$, constitutes a 1-of-K
classification task over soft labels $\xi_{k,s}^{\Theta^{i-1}}$, where the
gate learns to approximate the posterior while the ARX sub-models compete to
get responsibility over samples, gaining specialization within different
regions of the regression space. During prediction, the sub-model with highest
probability is selected. It is worth noting that, as opposed to separated
clustering, such an integrated discriminant strategy foster the exploitation
of network flexibility on the decision surface between the classes rather than
on the overall distribution [20], often resulting in improved classification
performances. Despite a closed form solution is not available for this
problem, standard network training algorithms can be exploited in order to
increase $\mathcal{Q}$, thus leading to a Generalized-EM approach which still
guarantee convergence [19]. Specifically, we deploy the classification task
from logits, increasing numerical stability with reference to the conventional
cross-entropy over softmax [21]. Then, the network weights are updated using
the Adam algorithm, conceived to tackle noisy and sparse gradients [22].

To avoid potential variance collapses [23], we adopt the Maximum A Posteriori
estimate proposed in [24], given by:

$\sigma_{s}^{2^{i}}=\frac{v^{2}/S+V_{s}^{i}}{\upsilon_{0}+\sum_{k=1}^{N}\xi_{k,s}^{\Theta^{i-1}}+D+2}$
(24)

where $v^{2}=(1/N)\sum_{k=1}^{N}(y_{k}-\bar{y})^{2}$,
$\bar{y}=(1/N)\sum_{k=1}^{N}y_{k}$,
$V^{i}_{s}=\sum_{k=1}^{N}\xi_{k,s}^{\Theta^{i-1}}\left(y_{k}-\theta_{s}^{iT}\varphi_{k}\right)\left(y_{k}-\theta_{s}^{iT}\varphi_{k}\right)^{T}$,
including the weakest prior $\upsilon_{0}=D+2$ with data dimension $D$. Since
not subjected to potential collapse, an overall standard deviation parameter
can be estimated by the usual MLE form:
$\sigma^{2^{i}}=\sum_{s=1}^{S}\sum_{k=1}^{N}\xi_{k,s}^{\Theta^{i-1}}\left(y_{k}-\theta_{s}^{iT}\varphi_{k}\right)^{2}/\sum_{s=1}^{S}\sum_{k=1}^{N}\xi_{k,s}^{\Theta^{i-1}}$.

The training algorithm proceed by alternating E/M-steps until convergence,
tested by stopping conditions on the likelihood or maximum number of
iterations. As common for EM, several initialization are required, by sampling
different values of the parameters to tackle eventual convergence to poor
local minima. As shown in the result section, we found useful experimentally
to sample the initial parameters of the network from relatively large values
and to initialize the sub-models biases using k-means [23], to avoid modes
collapse.

##  IV Numerical results

To test the proposed method, we considered a slight modification of the widely
adopted case study proposed in [7], assumed here to introduce nonlinear
partitioning requirements. Besides, it has been conceived as an extension to
the simpler nonlinearly piecewise affine map estimation problem considered in
[16]. Specifically, the system is constituted by the following components:

$\displaystyle\begin{dcases}\text{if}&4y_{k-1}-u_{k-1}+10<0\ :\\\
&y_{k}=-0.4y_{k-1}+u_{k-1}+1.5+e_{k}\\\ \text{if}&4y_{k-1}-u_{k-1}+10\geq 0\
\text{and}\ 5y_{k-1}+u_{k-1}-6\leq 0\ :\\\
&y_{k}=0.5y_{k-1}-u_{k-1}-0.5+e_{k}\\\ \text{if}&5y_{k-1}+u_{k-1}-6>0\ :\\\
&y_{k}=-0.4y_{k-1}+u_{k-1}+1.5+e_{k}\end{dcases}$

with $e_{k}\sim N(0,0.2^{2})$ and input sampled from a uniform distribution
$u_{k}\sim U[-4,4]$. The aim is to infer the same dynamics occurring in two
parts of the regressor space, by reconstructing two modes. We might remark
that, even if this system is defined by two affine maps, identification
methods based on linear partition of the regression space would require three
modes as the regions are not linearly separable.

We generated a sequence 6000 samples. The first 5000 has been devoted to
training/validation and the last 1000 represents the test set. To configure
the hyperparameters, we used a cross-validation procedure, leaving 1000
samples to validation.

Both the overall model architecture and the EM-based estimation algorithm have
been developed in Numpy-1.18. The weighted least squares subproblems are
implemented using the Linear Regression model of Scikit-learn-0.23 while the
neural network using Tensorflow-2.2 API.

By cross-validation, we found experimentally that a neural network with 10
$tanh$ hidden units is sufficient to properly fit the problem at hand. The
maximum number of iterations has been set to 500, with a stopping tolerance on
the log-likelihood variation of 1e-4. The Adam algorithm starts with a
learning rate of 0.01, running for 3 epochs mini-batches of 100 samples in
each M-iteration. The network starts from random normal weights with zero-
mean. Specifically, we observed that the common weight initialization to small
(i.e., approximately zero) random values frequently results in convergence to
poor local solutions, where the modes collapse to a common (i.e., average)
mode. To avoid this issue, we set the standard deviation of the network
initializer to a high value (i.e., 10), to stimulate the system to start from
initial quite different local configurations, thus limiting mode collapse.
Using such configuration, we found that a maximum number of 5 random
executions of the estimation algorithm are sufficient to reach good local
solutions.

To evaluate the accuracy of the ARX models reconstruction and of the estimated
mode sequence, we employed [25]:

$\small\mathcal{F}_{\theta}=\frac{1}{S}\sum_{s=1}^{S}\left(1-\frac{\lVert\theta_{s}-\hat{\theta}_{s}\lVert}{\lVert\theta_{s}\lVert}\right),\
\mathcal{F}_{s}=\frac{1}{N_{test}}\sum_{k=1}^{N_{test}}\mathbf{1}_{\hat{s}_{k}=s_{k}}$
(25)

namely the parameters and mode fit indexes, where $\hat{\theta}_{s}$,
$\hat{s}_{k}$ indicates the estimated ARX parameters and state respectively.

As the order of the discrete states in the model do not necessarily match that
of the target system, post-processing is required before evaluation [19].
Following [25], without loss of generality, we reordered the sub-models
estimate by the Euclidean norm to the target parameters.

The obtained estimates for the ARX submodels are reported in Table I.
Statistics are computed over 100 runs of the estimation algorithm using the
same hyperparameters configuration.

TABLE I: Identified ARX parameters (100 trials) | $True$ | $Mean$ | $Std$
---|---|---|---
$\theta_{1}$ | [1.5, -0.4, 1] | [1.503, -0.400, 0.999] | [2.6e-4, 7.9e-5, 8.1e-5]
$\theta_{2}$ | [-0.5, 0.5, -1] | [-0.499, 0.499, -0.995] | [4.3e-4, 3.7e-4, 4.2e-5]
$\sigma^{2}$ | 0.2 | 0.201 | 3.1e-4

Figure 1 reports the evolution of the estimated ARX parameters and covariance
during the iterations of 5 consecutive runs of the training algorithm. Under
the adopted configuration, convergence starting from random initial conditions
have been observed in approximately 60 iterations on average. The average
computation time of one iteration for this case study on an PC with
CPU-i7-2.5-GHz-RAM-8Gb has been 0.42 seconds.

Figure 1: ARX parameters during iterations

The obtained performance indexes are: $\mathcal{F}_{\theta}=0.997$,
$\mathcal{F}_{s}=0.996$, while Figure 2 reports the residuals over the test
set.

Notably, we obtained high performances over the test set. The few visible
spikes in the residual plot are due to small errors in the estimated switching
surface for the target discontinuous PWARX, which are in general inevitable as
identification is performed over finite datasets [7]. The estimated partitions
of the regressor space are shown in Figure 3, where the different colors
represents the states (i.e, the maximum gate network activation for each
sample) and the black lines the true partitions of the target system.

Finally, we might remark that the same estimation algorithm can be applied on
a simplified discriminant function, including a linear parametrization. In
this case, an EM-based estimation of the PrARX model is obtained, enabling
transformation to the equivalent PWARX. However, three discrete modes would be
required to enable linear separation.

Figure 2: Residuals over the test data set Figure 3: Discrete states
classification in the regressor space

##  V Conclusion and Next Steps

In this work we focused on the identification of Nonlinarly Piecewise Affine
systems, characterized by arbitrary regions not restricted to polyhedral
partitions. To this end, we augmented the Probability weighted ARX model,
following a Mixture of Expert approach, by parametrizing the discriminant
function through a neural network conditioned by the regressors. As the direct
maximization of the likelihood function is particularly difficult when complex
models and conditioning functions are employed, we leveraged on Expectation
Maximization. Hence, ARX submodels fitting and classification tasks are
disentangled within the iterations, enabling a closed form solution of the
regression problems while increasing the log-likelihood at each iteration. By
application to a test case we showed the capability of the proposed approach
to achieve accurate parameters estimation and discrete state prediction in out
of sample conditions.

Next developments will include the integration of automatic techniques for
both model orders and number of modes selection, the investigation of a
Bayesian inference approach and the application to a real application case.

## References

  * [1] Y. Yuan, X. Tang, W. Zhou, W. Pan, X. Li, H.-T. Zhang, H. Ding, and J. Goncalves, “Data driven discovery of cyber physical systems,” Nature Communications, vol. 10, no. 1, p. 4894, 2019.
  * [2] A. Garulli, S. Paoletti, and A. Vicino, “A survey on switched and piecewise affine system identification,” 16th IFAC Symposium on System Identification, vol. 45, no. 16, 2012.
  * [3] F. Lauer and G. Bloch, Hybrid system identification: Theory and algorithms for learning switching models, vol. 478.  Springer, 2018.
  * [4] J. Lunze and F. Lamnabhi-Lagarrigue, Handbook of Hybrid Systems Control: Theory, Tools, Applications.  Cambridge University Press, 2009.
  * [5] S. Paoletti, A. L. Juloski, G. Ferrari-Trecate, and R. Vidal, “Identification of hybrid systems a tutorial,” Eur.Jour.of Control, pp. 242 – 260, 2007\.
  * [6] S. Paoletti, I. Savelli, A. Garulli, and A. Vicino, “A bilevel programming framework for piecewise affine system identification,” in 2019 IEEE 58th Conference on Decision and Control (CDC), pp. 7376–7381, 2019\.
  * [7] A. Bemporad, A. Garulli, S. Paoletti, and A. Vicino, “A bounded-error approach to piecewise affine system identification,” IEEE Transactions on Automatic Control, vol. 50, no. 10, pp. 1567–1580, 2005.
  * [8] A. L. Juloski, S. Weiland, and W. P. M. H. Heemels, “A bayesian approach to identification of hybrid systems,” IEEE Transactions on Automatic Control, vol. 50, no. 10, pp. 1520–1533, 2005.
  * [9] H. Ohlsson and L. Ljung, “Identification of switched linear regression models using sum-of-norms regularization,” Automatica, vol. 49, no. 4, pp. 1045 – 1050, 2013.
  * [10] G. Ferrari-Trecate, M. Muselli, D. Liberati, and M. Morari, “A clustering technique for the identification of piecewise affine systems,” Automatica, vol. 39, no. 2, pp. 205 – 217, 2003.
  * [11] H. Nakada, K. Takaba, and T. Katayama, “Identification of piecewise affine systems based on statistical clustering technique,” Automatica, vol. 41, no. 5, pp. 905 – 913, 2005.
  * [12] V. Breschi, D. Piga, and A. Bemporad, “Piecewise affine regression via recursive multiple least squares and multicategory discrimination,” Automatica, vol. 73, pp. 155–162, 11 2016.
  * [13] S. Taguchi, T. Suzuki, S. Hayakawa, and S. Inagaki, “Identification of probability weighted multiple arx models and its application to behavior analysis,” in Proceedings of the 48h IEEE Conference on Decision and Control (CDC), pp. 3952–3957, 2009.
  * [14] S. E. Yuksel, J. N. Wilson, and P. D. Gader, “Twenty years of mixture of experts,” IEEE Tr. on Neural Networks and Learning Systems, 2012.
  * [15] N. Ahmed and M. Campbell, “Variational learning of autoregressive mixtures of experts for fully bayesian hybrid system identification,” in Proceedings of the 2011 American Control Conf., pp. 139–144, 2011.
  * [16] F. Lauer and G. Bloch, “Switched and piecewise nonlinear hybrid system identification,” in Hybrid Systems: Computation and Control, 2008.
  * [17] M. Mejari, V. Breschi, and D. Piga, “Recursive bias-correction method for identification of piecewise affine output-error models,” IEEE Control Systems Letters, vol. 4, no. 4, pp. 970–975, 2020.
  * [18] A. S. Weigend, M. Mangeas, and A. N. Srivastava, “Nonlinear gated experts for time series: discovering regimes and avoiding overfitting,” International journal of neural systems, vol. 6 4, pp. 373–99, 1995.
  * [19] S. Fruhwirth-Schnatter, G. Celeux, and C. P. Robert, Handbook of Mixture Analysis.  Chapman and Hall/CRC, 2019.
  * [20] Y. Bengio and P. Frasconi, “An input output hmm architecture,” Advances in Neural Information Processing Systems, vol. 7, 12 1995.
  * [21] A. Géron, Hands-on machine learning with Scikit-Learn and TensorFlow.  O’Reilly Media, 2017.
  * [22] D. Kingma and J. Ba, “Adam: A method for stochastic optimization,” International Conference on Learning Representations, 12 2014.
  * [23] K. P. Murphy, Machine learning: a probabilistic perspective.  MIT, 2012.
  * [24] C. Fraley and A. Raftery, “Bayesian regularization for normal mixture estimation and model-based clustering,” Journal of Classification, vol. 24, pp. 155–181, 02 2007.
  * [25] G. Pillonetto, “A new kernel-based approach to hybrid system identification,” Automatica, vol. 70, pp. 21 – 31, 2016.