# On Density Functional Theory models for one-dimensional homogeneous
materials

Bouchra Bensiali  Ecole Centrale Casablanca, Bouskoura, Ville Verte, P.B.
27182, Morocco<EMAIL_ADDRESS>,  Salma Lahbabi
EMAMI, LRI, ENSEM, UHIIC, 7 Route d’El Jadida, B.P. 8118 Oasis, Casablanca;
UM6P-College of Computing, Lot 660, Hay Moulay Rachid Ben Guerir, 43150,
Morocco<EMAIL_ADDRESS>,  Abdallah Maichine  LAMA, Department of
Mathematics, Faculty of Sciences, Mohammed V University in Rabat, 4 Av. Ibn
Battouta, 1014, Rabat, Morocco<EMAIL_ADDRESS>and  Othmane Mirinioui
Ecole Centrale Casablanca, Bouskoura, Ville Verte, P.B. 27182, Morocco
<EMAIL_ADDRESS>
###### Abstract.

This paper studies DFT models for homogeneous 1D materials in the 3D space. It
follows our previous work about DFT models for homogeneous 2D materials in 3D.
We show how to reduce the problem from a 3D energy functional to a _2D_ energy
functional. The kinetic energy is treated as in [7, 8] by diagonilizing
admissible states, and writing the kinetic energy as the infimum of a modified
kinetic energy functional on reduced states. Besides, we treat here the
Hartree interaction term in 2D, and show how to properly define the mean-field
potential, through Riesz potential. We then show the well posdness of the
reduced model and present some numerical illustrations.

###### Contents

  1. 1 Introduction
  2. 2 2D Hartree interaction
  3. 3 Model reduction
  4. 4 Well-posdness of the reduced model
  5. 5 Numerical simulation
    1. 5.1 Thomas-Fermi model
    2. 5.2 Numerical schemes
    3. 5.3 Numerical results

##  1\. Introduction

In recent years, low-dimensional materials, such as 2D sheets, 1D nanowires
and quantum dots have become of great importance in condensed matter physics
due to their distinctive features, which are relevant in many applications [6,
19, 22, 3, 9, 1, 26]. In order to investigate their electronic properties,
such as electron’s mobility and band gap, Density Functional Theory (DFT)
calculations are often used [4, 24, 11, 18, 20, 25]. Our objective is to prove
and study DFT models for low-dimensional systems in order to allow low
computational cost. In [2], the Thoma-Fermi-von Weisäcker model for 1D and 2D
materials is derived by means of thermodynamic limit procedure. The works [7,
8] have initiated the study of DFT models for low dimensional materials that
are homogeneous, and we were particularly interested in 2D materials embedded
in the 3D space. In this paper, we focus on one–dimensional materials that
have an atomic length in two dimensions of the space, and a crystalline
structure in the remaining dimension.

In [7], 2D materials are considered with a nuclear charge distribution of the
form

$\mu_{2D}(x_{1},x_{2},x_{3})=\mu_{2D}(x_{3})$

that depends only on the orthogonal direction $x_{3}$. The electrons are
described by one body density matrices $\gamma$ that satisfy Pauli principle
$0\leq\gamma\leq 1$. As the problem is convex, $\gamma$ share the same
symmetry as $\mu$,

(1.1) $U_{R}\gamma=\gamma U_{R},\quad\forall R\in\mathbb{R}^{2},$

where $U_{R}f(x)=f(x_{1}-R_{1},x_{2}-R_{2},x_{3})$ is the translation operator
along the first two variables. Note that the density associated to a $\gamma$
satisfying (1.1) depends only on the third variable
$\rho_{\gamma}(x)=\rho_{\gamma}(x_{3})$ and the trace per unit surface of
$\gamma$ can be defined by $\underline{{\rm
Tr}}_{2D}(\gamma)=\int_{\mathbb{R}}\rho_{\gamma}$. In the reduced Hartree-Fock
(rHF) model, the energy of the system is of the form

$\underline{{\mathcal{E}}}_{2D}(\gamma)=\frac{1}{2}\underline{{\rm
Tr}}_{2D}\left(-\Delta_{3}\gamma\right)+{\mathcal{D}}_{1}(\rho_{\gamma}-\mu),$

where $\Delta_{d}$ denotes the the Laplacian operator on
$L^{2}(\mathbb{R}^{d})$, and ${\mathcal{D}}_{1}$ is the 1D Hartree
interaction. One main result of [7] is a thorough discussion of this term and
the definition of the mean-field potential (see [7, Proposition 3.3]). We
prove that this problem is equivalent to a one-dimensional problem set on
operators $G$ acting on $L^{2}(\mathbb{R})$, which are positive $G\geq 0$, but
which no longer need to satisfy Pauli principle. Instead, an additional term
appears in the energy: the three-dimensional kinetic energy
$\frac{1}{2}\underline{{\rm Tr}}_{2D}(-\Delta_{3}\gamma)$ is replaced by a
one-dimensional kinetic energy of the form

$\frac{1}{2}{\rm Tr}(-\Delta_{1}G)+\pi{\rm Tr}(G^{2}).$

Although we restrict ourselves to the reduced Hartree-Fock model for
simplicity, similar derivations can be performed for general Kohn–Sham models
and one can consider models of the form

$\inf\left\\{\frac{1}{2}{\rm Tr}(-\Delta_{1}G)+\pi{\rm
Tr}(G^{2})+{\mathcal{D}}_{1}(\rho_{G}-\mu)+E^{\rm xc}(\rho_{G})\right\\}.$

In [8], we have considered a uniform magnetic field which is perpendicular to
the 2D material axis. The reduction here is obtained considering operators
commuting with _magnetic translations_. We first provide a decomposition of
such operators and then state the equivalent one-dimensional problem in which
the kinetic energy per unit surface $\frac{1}{2}\underline{{\rm
Tr}}_{2D}(\mathrm{L}_{\mathrm{A}}\gamma)$ has been repalaced by

$\frac{1}{2}{\rm Tr}(-\Delta_{1}G)+{\rm Tr}(F(b,G)),$

where $\mathrm{L}_{\mathrm{A}}$ denotes the three-dimensional Landau operator,
$b$ the magnetic field strenghten, and $F(b,\cdot)$ a suitable piecewise
linear function, see for instance [8, Theorem 3.1].

In the present paper, we consider homogeneous 1D materials, i.e. nuclear
charge distributions of the form

$\mu(x_{1},x_{2},x_{3})=\mu(x_{1},x_{2})$

that do not depend on the third variable $x_{3}$, the axis of the material. As
we consider convex models, we assume that the electronic density matrices
follow the same symmetries

(1.2) $\tau_{R}\gamma=\gamma\tau_{R},\quad\forall R\in\mathbb{R},$

where $\tau_{R}f(x_{1},x_{2},x_{3})=f(x_{1},x_{2},x_{3}-R)$ denotes the
translation operator by the vector $Re_{3}$. Again, states satisfying the
symmetry (1.2) have densities that depend only on $x_{1}$ and $x_{2}$ and the
trace per unit length of $\gamma$ can be defined by $\underline{{\rm
Tr}}(\gamma)=\int_{\mathbb{R}^{2}}\rho_{\gamma}$. In the reduced Hartree-Fock
model, the energy of the system is

(1.3) $\underline{{\mathcal{E}}}(\gamma)=\frac{1}{2}\underline{{\rm
Tr}}\left(-\Delta_{3}\gamma\right)+\tilde{{\mathcal{D}}}_{2}(\rho_{\gamma}-\mu),$

where $\tilde{{\mathcal{D}}}_{2}$ is the 2D Hartree interaction. Formally, the
interaction term ${\mathcal{D}}_{2}$ is defined by

(1.4)
$\tilde{{\mathcal{D}}}_{2}(f,g)=-2\int_{\mathbb{R}^{2}\times\mathbb{R}^{2}}{f(x)}\log(\left|x-y\right|)g(y).$

This definition is only valid for functions that decay fast enough at
infinity. Besides, it is unclear whether this bilinear form
$\tilde{{\mathcal{D}}}_{2}$ is bounded from below or not. One key result of
this paper is a regularization of this expression adapted to neutral systems;
and a definition of the corresponding mean-field potential (see Section 2).

For the kinetic part, following the same arguments as in [7], we prove that
this problem is equivalent to a two-dimensional problem, where the three-
dimensional kinetic energy $\frac{1}{2}\underline{{\rm
Tr}}(-\Delta_{3}\gamma)$ is replaced by a two-dimensional kinetic energy of
the form

$\frac{1}{2}{\rm Tr}(-\Delta_{2}G)+\frac{\pi}{6}{\rm Tr}(G^{3}),$

defined on ${\mathcal{G}}:=\\{G\in{\mathbb{S}}(L^{2}(\mathbb{R}^{2})):G\geq
0,\;{\rm Tr}(G)<{\infty}\\}$. The reduced 2D rHF energy then reads

(1.5) ${\mathcal{E}}(G)=\frac{1}{2}{\rm Tr}(-\Delta_{1}G)+\frac{\pi}{6}{\rm
Tr}(G^{3})+{\mathcal{D}}_{2}(\rho_{G}-\mu).$

We show that this problem of minimizing the energy (1.5) over a suitable set
of admissible states admits a unique minimizer, called ground state. We also
derive the self-consistent equation that this ground state needs to satisfy.

Finally, we perform numerical simulations on the simpler Thomas–Fermi model,
where the kinetic energy is expressed as an explicit functional of the density
$\int\rho^{5/3}$. This allows us to exhibit the properties of the regularized
2D Hartree interaction with respect to the classical expression (1.4).

This article is organized as follow. Section 2 is devoted to the
regularization of the 2D Hartree interaction and the definition of the
corresponding mean-field potential. In Section 3, we show how to reduce the 3D
problem (1.3) into the 2D model (1.5), that we study in Section 4. Finally, in
Section 5, we consider the orbital free Thomas-Fermi model as a semi–classical
limit of the rHF model in order to perform elementary numerical illustrations.
Acknowledgments. The research leading to these results has received funding
from OCP grant AS70 “Towards phosphorene based materials and devices”. S.
Lahbabi thanks the CEREMADE for hosting her during the final writing of this
article.

##  2\. 2D Hartree interaction

For 1D crystals, the Hartree interaction kernel has been derived in [2] by
means of thermodynamic limit procedure. It is given by

$G(x)=-2\pi\log(\left|(x_{1},x_{2})\right|)+\sum_{k\in\mathbb{Z}}\left(\frac{1}{\left|x-ke_{3}\right|}-\int_{-1/2}^{1/2}\frac{dt}{\left|x-(y,k)e_{3}\right|}\right),\quad
x\in\mathbb{R}^{3},$

and the Hartree interaction energy is formally given by

${\mathcal{D}}_{\rm per}(f,g)=\int_{\Gamma\times\Gamma}f(x)G(x-y)g(y),$

where $\Gamma=\mathbb{R}^{2}\times\left[-1/2,1/2\right]$. For functions $f,g$
which do not depend on the third variable $x_{3}$, this interaction reduces to

(2.1)
${\mathcal{D}}_{2}(f,g)=-2\int_{\mathbb{R}^{2}\times\mathbb{R}^{2}}{f(x)}\log(\left|x-y\right|)g(y).$

The integral in (2.1) is only valid for functions decaying fast enough and the
map $f\mapsto{\mathcal{D}}_{2}(f)$ is not convex in general. The aim of this
section is to give a regularization of this expression suitable for neutral
charge distributions, which, in our context, will be $f=\rho-\mu$. We
therefore adopt the following definition

(2.2)
${\mathcal{D}}(f,g):=4\pi\int_{\mathbb{R}^{2}}\dfrac{\overline{{\mathcal{F}}_{2}(f)(k)}{\mathcal{F}}_{2}(g)(k)}{|k|^{2}}{\mathrm{d}}k,$

where $f$ and $g$ belong to the following Coulomb space

${\mathcal{C}}:=\left\\{f\in
L^{1}(\mathbb{R}^{2}),\;k\mapsto\frac{{\mathcal{F}}_{2}(f)(k)}{|k|}\in
L^{2}(\mathbb{R}^{2})\right\\}.$

We have denoted by ${\mathcal{F}}_{2}$ the Fourier transform in
$\mathbb{R}^{2}$ defined by
${\mathcal{F}}_{2}f(k)=\frac{1}{2\pi}\int_{\mathbb{R}^{2}}f(x)e^{-ikx}dx$. We
note that the integrability of
$k\mapsto\frac{|{\mathcal{F}}_{2}(f)(k)|^{2}}{|k|^{2}}$ over $\mathbb{R}^{2}$
implies that

${\mathcal{F}}_{2}(f)(0)=\int_{\mathbb{R}^{2}}f(x){\mathrm{d}}x=0,$

which means that elements of ${\mathcal{C}}$ are neutral distributions.
Conversely, for neutral distributions that decay fast enough,
${\mathcal{D}}_{2}$ and ${\mathcal{D}}$ coincide.

######  Lemma 2.1.

If $f\in L^{1}(\mathbb{R}^{2})$ is such that $\int f=0$ and
$\log(1+\left|\cdot\right|)f\in L^{1}(\mathbb{R}^{2})$, then

${\mathcal{D}}(f)={\mathcal{D}}_{2}(f).$

The proof of this lemma can be found in [5].

For $f\in{\mathcal{C}}$, let us introduce ${\mathrm{W}}_{f}$ as

${\mathrm{W}}_{f}:={\mathcal{F}}_{2}^{-1}\left(k\mapsto\frac{{\mathcal{F}}_{2}(f)(k)}{|k|}\right)=(-\Delta_{2})^{-\frac{1}{2}}f\quad\in
L^{2}(\mathbb{R}^{2}).$

Thanks to Parseval identity, one can rewrite the interaction ${\mathcal{D}}$
as

${\mathcal{D}}(f,g)=4\pi\int_{\mathbb{R}^{2}}\overline{{\mathrm{W}}_{f}(x)}{W_{g}(x)}{\mathrm{d}}x,\qquad\forall
f,g\in{\mathcal{C}}.$

An explicit expression of ${\mathrm{W}}_{f}$, for $f\in{\mathcal{C}}$, is
given through the Riesz potential as follows (see for instance [23, Chapter
V]).

######  Lemma 2.2.

Let $f\in{\mathcal{C}}$. Then,

(2.3)
${\mathrm{W}}_{f}(x)=\frac{1}{2\pi}\int_{\mathbb{R}^{2}}\frac{f(y)}{|x-y|}{\mathrm{d}}y.$

The integral in (2.3) is absolutely convergent whenever $f\in
L^{p}(\mathbb{R}^{2})$, for some $1\leq p<2$ (see for instance [23, Theorem 1
of Chapter V]). Furthermore, if $p\in(1,2)$ and $q=\frac{2p}{2-p}$. Then,

(2.4) $\|{\mathrm{W}}_{f}\|_{q}\leq C\|f\|_{p},\qquad\forall f\in
L^{p}(\mathbb{R}^{2}),$

for a suitable constant $C$ which is independent of $f$.

We turn now to the definition of a _mean-field potential_ ${\mathrm{V}}_{f}$
associated to a charge distribution $f\in{\mathcal{C}}$. It needs to satisfy

${\mathcal{D}}(f,g)=\int_{\mathbb{R}^{2}}{\mathrm{V}}_{f}(x)g(x){\mathrm{d}}x,\qquad\forall
g\in{\mathcal{C}}$

and

$-\Delta{\mathrm{V}}_{f}=4\pi f$

in a weak sense. One would consider ${\mathrm{V}}_{f}=4\pi(-\Delta)^{-1}f=4\pi
W_{W_{f}}$. However, the integral in (2.3) is not necessarily convergent in
the critical case of $p=2$. We suggest an alternative way to define the
potential ${\mathrm{V}}_{f}$, for $f\in{\mathcal{C}}$ through a modification
of the Riesz potential suggested in [13]. Let us define the following two
functions

$q_{1}(x)=|x|\left(1+\log^{+}\left(\frac{1}{|x|}\right)\right)\quad\text{and}\quad
q_{2}(x)=|x|\left(1+\log^{+}\left(|x|\right)\right),$

where

$\log^{+}(|x|):=\log(|x|){\mathds{1}}_{\\{|x|\geq
1\\}}=\begin{cases}\log(|x|)\quad\text{if}\quad|x|\geq 1\\\
0\quad\text{else}\end{cases}.$

In the following proposition, we define and state some properties of the
modified Riesz potential. For the proof and further details, we refer to [13,
Corollary 3.23].

######  Proposition 2.3.

For $g\in L^{2}(\mathbb{R}^{2})$, define

(2.5)
$\mathcal{I}_{1}^{g}(x):=\frac{1}{2\pi}\int_{|y|<1}\frac{g(y)}{|x-y|}\;{\mathrm{d}}y\quad\text{and}\quad\mathcal{I}_{2}^{g}(x):=\frac{1}{2\pi}\int_{|y|\geq
1}\left(\frac{1}{|x-y|}-\frac{1}{|y|}\right)g(y)\;{\mathrm{d}}y.$

One has the following

  1. a)

${\mathcal{I}}_{1}^{g}$ and ${\mathcal{I}}_{2}^{g}$ are well-defined, i.e.,
the integrals in (2.5) are convergent.

  2. b)

There exist two real numbers $C_{1}$ and $C_{2}$ such that

(2.6)
$\displaystyle\left\|\frac{\mathcal{I}_{1}^{g}}{q_{1}}\right\|_{L^{2}}\;\;\leq\;C_{1}\left\|g\right\|_{L^{2}}\quad
and\quad\left\|\frac{\mathcal{I}_{2}^{g}}{q_{2}}\right\|_{L^{2}}\;\;\leq\;C_{2}\left\|g\right\|_{L^{2}}.$

As a consequence of Proposition 2.3, we have $\mathcal{I}_{i}^{f}\in
L^{1}_{\rm loc}(\mathbb{R}^{2})$, for $i\in\left\\{1,2\right\\}$. Now, for
$f\in{\mathcal{C}}$, we define the mean field potential as

(2.7)
${\mathrm{V}}_{f}:=4\pi\left(\mathcal{I}_{1}^{{\mathrm{W}}_{f}}+\mathcal{I}_{2}^{{\mathrm{W}}_{f}}\right).$

The following theorem is one of the main results of the paper and it allows to
rewrite the potential energy ${\mathcal{D}}$ in terms of the mean–field
potential $V_{f}$.

######  Theorem 2.4.

Let $f\in{\mathcal{C}}$. Then

(2.8)
$\mathcal{D}(f,h)=\int_{\mathbb{R}^{2}}{\mathrm{V}}_{f}(x)\,h(x){\mathrm{d}}x\qquad\forall
h\in\mathcal{C},$

and ${\mathrm{V}}_{f}$ satisfies

(2.9) $-\Delta{\mathrm{V}}_{f}=4\pi f$

in a distributional sense.

###### Proof.

Let us first prove (2.8) for $h$ in the Schwartz space of fast decaying
functions $\mathcal{S}(\mathbb{R}^{2})$. Since
${\mathcal{I}}_{j}^{{\mathrm{W}}_{f}}h=(q_{j}^{-1}{\mathcal{I}}_{j}^{{\mathrm{W}}_{f}})\,q_{j}h\in
L^{1}(\mathbb{R}^{2})$, for $j\in\\{1,2\\}$. Then, one can write

$\displaystyle\frac{1}{4\pi}\int_{\mathbb{R}^{2}}{\mathrm{V}}_{f}(x)\,h(x)\,{\mathrm{d}}x:=$
$\displaystyle\int_{\mathbb{R}^{2}}\mathcal{I}_{1}^{{\mathrm{W}}_{f}}(x)h(x)\,{\mathrm{d}}x+\int_{\mathbb{R}^{2}}\mathcal{I}_{2}^{{\mathrm{W}}_{f}}(x)h(x)\,{\mathrm{d}}x$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{\mathbb{R}^{2}}\int_{\\{|y|<1\\}}\frac{{\mathrm{W}}_{f}(y)}{|x-y|}h(x)\,{\mathrm{d}}y{\mathrm{d}}x$
(2.10) $\displaystyle+\frac{1}{2\pi}\int_{\mathbb{R}^{2}}\int_{\\{|y|\geq
1\\}}\left(\frac{1}{|x-y|}-\frac{1}{|y|}\right){\mathrm{W}}_{f}(y)h(x)\,{\mathrm{d}}y{\mathrm{d}}x.$

In order to apply Fubini theorem and change the order of the integrals in the
above double intgrals, we have to check the absolute integrability. Regarding
the first term, one can note that
$(|q_{1}|)^{-1}\mathcal{I}_{1}^{|{\mathrm{W}}_{f}|}$ is in
$L^{2}(\mathbb{R}^{2})$, thanks to (2.6), and so is $q_{1}h$. Therefore, one
has

$\displaystyle\frac{1}{2\pi}\int_{\mathbb{R}^{2}}\int_{\\{|y|<1\\}}\frac{{\mathrm{W}}_{f}(y)}{|x-y|}h(x)\,{\mathrm{d}}y{\mathrm{d}}x$
$\displaystyle=\frac{1}{2\pi}\int_{|y|<1}\int_{\mathbb{R}^{2}}\frac{{\mathrm{W}}_{f}(y)}{|x-y|}h(x)\,{\mathrm{d}}x{\mathrm{d}}y$
(2.11)
$\displaystyle=\int_{\\{|y|<1\\}}{\mathrm{W}}_{f}(y){\mathrm{W}}_{h}(y)\;{\mathrm{d}}y.$

Taking into the account that the estimate in (2.6) holds true if one
substitutes the kernel defining $\mathcal{I}_{2}$ by its absolute value (see
for instance [13, Theorem 3.19]), then the last double intergal of (2) is
absolutely convergent as well. Thus,

$\displaystyle\frac{1}{2\pi}\int_{\mathbb{R}^{2}}\int_{|y|\geq
1}\left(\frac{1}{|x-y|}-\frac{1}{|y|}\right){\mathrm{W}}_{f}(y)h(x)\,{\mathrm{d}}y{\mathrm{d}}x$
$\displaystyle\qquad\qquad=\frac{1}{2\pi}\int_{|y|\geq
1}\int_{\mathbb{R}^{2}}\left(\frac{1}{|x-y|}-\frac{1}{|y|}\right){\mathrm{W}}_{f}(y)h(x)\,{\mathrm{d}}x{\mathrm{d}}y$
$\displaystyle\qquad\qquad=\int_{|y|\geq
1}{\mathrm{W}}_{f}(y)\left({\mathrm{W}}_{h}(y)-\frac{1}{2\pi|y|}\int_{\mathbb{R}^{2}}h\right)\,{\mathrm{d}}y$
(2.12) $\displaystyle\qquad\qquad=\int_{|y|\geq
1}{\mathrm{W}}_{f}(y){\mathrm{W}}_{h}(y)\,{\mathrm{d}}y.$

By adding (2.11) et (2.12), we obtain

$\int_{\mathbb{R}^{2}}{\mathrm{V}}_{f}(x)\,h(x)\,dx\>=4\pi\int_{\mathbb{R}^{2}}{\mathrm{W}}_{f}(x){\mathrm{W}}_{h}(x)\,dx.$

Now, let $h\in{\mathcal{C}}$. Then, $h\in L^{1}(\mathbb{R}^{2})$ and
${\mathrm{W}}_{h}\in L^{2}(\mathbb{R}^{2})$. By density of the Schwartz space
${\mathcal{S}}$, one can find $(h_{n})_{n}\subset{\mathcal{S}}$ with
$\int_{\mathbb{R}^{2}}h_{n}\,{\mathrm{d}}x=0$ such that $h_{n}\to h$ in
$L^{1}(\mathbb{R}^{2})$ and
$\frac{{\mathcal{F}}_{2}(h_{n})}{|\cdot|}\to\frac{{\mathcal{F}}_{2}(h)}{|\cdot|}$
in $L^{2}(\mathbb{R}^{2})$. Hence, one obtains (2.8), for $h\in{\mathcal{C}}$,
by letting $n\to{\infty}$.

We move now to the proof of (2.9). If one denotes by $(\cdot,\cdot)$ the
distributional brakets, then, for $\phi\in C^{\infty}_{c}(\mathbb{R}^{2})$,

$\displaystyle(-\Delta_{2}V_{f},\phi)=(V_{f},-\Delta_{2}\phi)=\int_{\mathbb{R}^{2}}V_{f}\Delta_{2}\phi\,{\mathrm{d}}x.$

One knows that $h:=-\Delta\phi$ satisfies $\int h=0$. Thus, by the previous
result

$\displaystyle(-\Delta_{2}V_{f},\phi)={\mathcal{D}}(f,h)$
$\displaystyle=4\pi\int_{\mathbb{R}^{2}}\frac{\overline{{\mathcal{F}}_{2}(f)(k)}{\mathcal{F}}_{2}(h)(k)}{\left|k\right|^{2}}{\mathrm{d}}k$
$\displaystyle=4\pi\int_{\mathbb{R}^{2}}\overline{{\mathcal{F}}_{2}(f)(k)}{\mathcal{F}}_{2}(\phi)(k){\mathrm{d}}k=\langle
4\pi f,\phi\rangle.$

∎

##  3\. Model reduction

In this section, we show that the 3D reduced Hartree Fock model for
homogeneous 1D crystal is equivalent to a 2D model. Other DFT models can be
handled similarly. We recall that the nuclear charge distribution $\mu\in
L^{1}_{\rm loc}(\mathbb{R}^{3})$ satisfies

(3.1) $\mu\geq
0,\quad\text{and}\quad\mu(x_{1},x_{2},x_{3})=\mu(x_{1},x_{2},x_{3}-R),\quad\forall
R\in\mathbb{R}.$

We assume that the charge per unit length is finite, that is,

$Z:=\int_{\mathbb{R}^{2}}\mu(x_{1},x_{2},0)\;{\mathrm{d}}x_{1}\,{\mathrm{d}}x_{2}<{\infty}.$

In Hartree–Fock type models, the state of electrons is described by a bounded
self–adjoint operator $\gamma$ satisfying Pauli exclusion principle
$0\leq\gamma\leq 1$. In our particular case, the system has the following
features.

  * •

It is equidistributed in one direction of the three–dimensional space,

  * •

it is confined in the remaining two directions with a finite number of
particles.

The electronic states $\gamma$ share the aforementioned invariance as follows

(3.2) $\gamma\tau_{R}=\tau_{R}\gamma,\qquad\forall R\in\mathbb{R},$

where $\tau_{R}$ denotes the translation operator by the vector $Re_{3}$
($\tau_{R}f(x)=f(x_{1},x_{2},x_{3}-R)$). The 3D rHF energy of a state $\gamma$
in

${\mathcal{K}}:=\\{\gamma\in\underline{{\mathfrak{S}}}(L^{2}(\mathbb{R}^{3})),\;0\leq\gamma\leq
1,\;\tau_{R}\gamma=\gamma\tau_{R},\;\forall R\in\mathbb{R}\\}.$

is

$\underline{{\mathcal{E}}}(\gamma)=\frac{1}{2}\underline{{\rm
Tr}}\left(-\Delta_{3}\gamma\right)+{\mathcal{D}}(\rho_{\gamma}-\mu),$

where we have denoted by $\underline{{\mathfrak{S}}}(L^{2}(\mathbb{R}^{3}))$
the space of self–adjoint, locally trace class operators acting on
$L^{2}(\mathbb{R}^{3})$, and ${\mathcal{D}}$ is the regularized 2D Hartree
interaction introduced in the previous section. We now define the 2D energy
functional for a state in

(3.3) ${\mathcal{G}}:=\\{G\in{\mathfrak{S}}(L^{2}(\mathbb{R}^{2})),\;G\geq
0\\},$

where we have denoted by ${{\mathfrak{S}}}(L^{2}(\mathbb{R}^{3}))$ the space
of self–adjoint, trace class operators acting on $L^{2}(\mathbb{R}^{2})$, by

(3.4) ${\mathcal{E}}(G)=\frac{1}{2}{\rm Tr}(-\Delta_{1}G)+\frac{\pi}{6}{\rm
Tr}(G^{3})+{\mathcal{D}}(\rho_{G}-\mu).$

The main result of the section is the following theorem.

######  Theorem 3.1.

The ground state energies

$\underline{I}=\inf\left\\{\underline{{\mathcal{E}}}(\gamma),\;\gamma\in{\mathcal{K}}\right\\}$

and

(3.5) $I=\inf\left\\{{\mathcal{E}}(G),\;G\in{\mathcal{G}}\right\\}$

are equal and the minimizers of both energies share the same densities, which
does not depend on $x_{3}$.

The proof of theorem (3.1) follows the same line as the one of [7, Theorem
2.8]. We highlight here some arguments in the 2D case.

_Diagonalization of admissible states._ We start by showing that states in
${\mathcal{K}}$ can be diagonalized. Let us introduce the partial Fourier
transform

$\displaystyle{\mathcal{U}}:L^{2}(\mathbb{R}^{3})$ $\displaystyle\to
L^{2}(\mathbb{R},L^{2}(\mathbb{R}^{2}))$ $\displaystyle f$
$\displaystyle\mapsto{\mathcal{F}}_{1}(f(x_{1},x_{2},\cdot))(k),$

where ${\mathcal{F}}_{1}$ refers to the Fourier transform in
$L^{1}(\mathbb{R})$
${\mathcal{F}}_{1}(f(x_{1},x_{2},\cdot))(k)=\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}{\mathrm{e}}^{-{\mathrm{i}}kt}f(x_{1},x_{2},t){\mathrm{d}}t$.
According to the Floquet–Bloch decomposition (see [21, Section XIII–16]), for
any $\gamma\in{\mathcal{K}}$ there exists a unique family of self–adjoint
operators $(\gamma_{k})_{k}$ on $L^{2}(\mathbb{R}^{2})$, called the fibers of
$\gamma$, such that

(3.6)
${\mathcal{U}}\gamma{\mathcal{U}}^{-1}=\int_{\mathbb{R}}^{\oplus}\gamma_{k}{\mathrm{d}}k,$

that is, for any $f=(f_{k})_{k\in\mathbb{R}}\in
L^{2}(\mathbb{R},L^{2}(\mathbb{R}^{2}))$,

$({\mathcal{U}}\gamma{\mathcal{U}}^{-1}f)_{k}=\gamma_{k}f_{k}.$

We point out that, if $\gamma\in{\mathcal{K}}$, then its fibers $\gamma_{k}$
satisfy Pauli principle $0\leq\gamma_{k}\leq 1$. The following proposition
links other properties of $\gamma\in{\mathcal{K}}$ to those of its fibers.

######  Proposition 3.2.

Let $\gamma=\int_{\mathbb{R}}^{\oplus}\gamma_{k}\in{\mathcal{K}}$ and suppose
that $\gamma$ is locally trace class. One has

  * •

$\gamma$ admits an integral kernel
$\gamma(\cdot;\cdot):(\mathbb{R}^{3})^{2}\to\mathbb{R}$ and

(3.7)
$\gamma(x_{1},x_{2},x_{3};y_{1},y_{2},y_{3})=\frac{1}{{2\pi}}\int_{\mathbb{R}}{\mathrm{e}}^{{\mathrm{i}}k(x_{3}-y_{3})}\gamma_{k}(x_{1},x_{2};y_{1},y_{2}){\mathrm{d}}k,$

where, for every $k\in\mathbb{R}$,
$\gamma_{k}(\cdot;\cdot):\mathbb{R}^{2}\to\mathbb{R}$ is the integral kernel
of $\gamma_{k}$.

  * •

Let $\rho_{\gamma}$ denote the density of $\gamma$, given by
$\rho_{\gamma}(x)=\gamma(x,x)$ for all $x\in\mathbb{R}^{3}$. Then,

$\rho(x_{1},x_{2},x_{3})=\frac{1}{{2\pi}}\int_{\mathbb{R}}\rho_{k}(x_{1},x_{2};x_{1},x_{2}){\mathrm{d}}k,$

where $\rho_{k}$ refers to the density of $\gamma_{k}$.
In particular, $\rho_{\gamma}$ does not depend on the third variable $x_{3}$.

  * •

The _trace per unit length_ of $\gamma$ is given by

$\underline{\rm
Tr}(\gamma):=\int_{\mathbb{R}^{2}}\rho_{\gamma}=\frac{1}{{2\pi}}\int_{\mathbb{R}}{\rm
Tr}(\gamma_{k}){\mathrm{d}}k.$

###### Proof.

Since $\gamma$ is locally trace class, then $\gamma_{k}$ is trace class, for
all $k\in\mathbb{R}$. By straightforward computation, one gets for a smooth
enough $f$ and $(x_{1},x_{2},k)\in\mathbb{R}^{3}$

$\displaystyle({\mathcal{U}}\gamma f)_{k}(x_{1},x_{2})$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}{\mathrm{e}}^{-{\mathrm{i}}kt}\,(\gamma
f)(x_{1},x_{2},t){\mathrm{d}}t$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}\int_{\mathbb{R}^{3}}{\mathrm{e}}^{-{\mathrm{i}}kt}\gamma(x_{1},x_{2},t;y_{1},y_{2},\tilde{t})f(y_{1},y_{2},\tilde{t})\;{\mathrm{d}}y_{1}{\mathrm{d}}y_{2}{\mathrm{d}}\tilde{t}{\mathrm{d}}t,$

On the other hand,

$\displaystyle(\gamma_{k}({\mathcal{U}}f)_{k})(x_{1},x_{2})$
$\displaystyle=\int_{\mathbb{R}^{2}}\gamma_{k}(x_{1},x_{2};y_{1},y_{2})({\mathcal{U}}f)_{k}(y_{1},y_{2}){\mathrm{d}}y_{1}{\mathrm{d}}y_{2}$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}^{2}}\int_{\mathbb{R}}{\mathrm{e}}^{-{\mathrm{i}}k\tilde{t}}\gamma_{k}(x_{1},x_{2};y_{1},y_{2})f(y_{1},y_{2},\tilde{t}){\mathrm{d}}y_{1}{\mathrm{d}}y_{2}{\mathrm{d}}\tilde{t}$

By identification,

$\gamma_{k}(x_{1},x_{2};y_{1},y_{2})=\int_{\mathbb{R}}{\mathrm{e}}^{-{\mathrm{i}}k(t-\tilde{t})}\gamma(x_{1},x_{2},t;y_{1},y_{2},\tilde{t}){\mathrm{d}}t.$

Therefore, applying an inverse Fourier transform, one obtains (3.7). The claim
on the density follows immediately by (3.7). ∎

_Diagonalization of the kinetic energy_ Let
$\gamma=\int_{\mathbb{R}}^{\oplus}\gamma_{k}\in{\mathcal{K}}$. Let
$f=(f_{k})_{k\in\mathbb{R}}\in L^{2}(\mathbb{R},L^{2}(\mathbb{R}^{2}))$ such
that $f_{k}$ is smooth enough for any $k$. Straightforward computation yields

$({\mathcal{U}}(-\Delta_{3})\gamma{\mathcal{U}}^{-1}f)_{k}=(-\Delta_{2}+k^{2})\gamma_{k}f_{k}.$

Hence, the kinetic energy per unit length of $\gamma$ can be written

$\frac{1}{2}\underline{\rm
Tr}(-\Delta_{3}\gamma)=\frac{1}{{4\pi}}\int_{\mathbb{R}}{\rm
Tr}((-\Delta_{2}+k^{2})\gamma_{k})\,{\mathrm{d}}k.$

_Reduced states_ Now, we map any admissible state $\gamma$ to a reduced 2D
state $G$ as follows

$\displaystyle\Lambda:{\mathcal{K}}$ $\displaystyle\to{\mathcal{G}}$
$\displaystyle\gamma$
$\displaystyle\mapsto\Lambda(\gamma):=\frac{1}{2\pi}\int_{\mathbb{R}}\gamma_{k}{\mathrm{d}}k.$

We notice that elements of ${\mathcal{G}}$, unlike those of ${\mathcal{K}}$,
are compact operators. We have the following result.

######  Proposition 3.3.

The map $\Lambda$ is onto. Furthermore, for any $\gamma\in{\mathcal{K}}$,

$\rho_{\gamma}(x_{1},x_{2},x_{3})=\rho_{\Lambda(\gamma)}(x_{1},x_{2}),\qquad\forall(x_{1},x_{2})\in\mathbb{R}^{2}.$

###### Proof.

Let $G\in{\mathcal{G}}$. One can write
$G=\sum_{j}g_{j}|\varphi_{j}\rangle\langle\varphi_{j}|$, with
$(\varphi_{j})_{j}$ an orthonormal family in $L^{2}(\mathbb{R}^{2})$ and
$g_{j}\geq 0$ such that $\sum_{j}g_{j}={\rm Tr}(G)<{\infty}$. We construct
$\gamma$ satisfying $\Lambda(\gamma)=G$ through its fibers $(\gamma_{k})_{k}$.
Let us define, for any $k\in\mathbb{R}$, $\gamma_{k}$ as

$\gamma_{k}=\sum_{j}\alpha_{j}^{(k)}|\varphi_{j}\rangle\langle\varphi_{j}|\quad\text{with}\quad\alpha_{j}^{(k)}=1_{\\{|k|\leq
s_{j}\\}},$

where, for any $j\in\mathbb{N}$, $s_{j}\geq 0$ will be determined later. One
has $0\leq\gamma_{k}\leq 1$ and
$\gamma_{k}\in{\mathfrak{S}}(L^{2}(\mathbb{R}^{2}))$. Thus,
$\gamma\in{\mathcal{K}}$. Moreover,
$\rho_{\gamma_{k}}=\sum_{j}\alpha_{j}^{(k)}|\varphi_{j}|^{2}$. Therefore,

$\displaystyle\rho_{\gamma}(x)$
$\displaystyle=\frac{1}{2\pi}\int_{\mathbb{R}}\rho_{\gamma_{k}}{\mathrm{d}}k$
$\displaystyle=\frac{1}{2\pi}\sum_{j}\int_{\mathbb{R}}1_{\\{|k|\leq
s_{j}\\}}|\varphi_{j}(x_{1},x_{2})|^{2}{\mathrm{d}}k=\frac{1}{\pi}\sum_{j}s_{j}|\varphi_{j}(x_{1},x_{2})|^{2}.$

Finally, with the choice $s_{j}=\pi{g_{j}}$, one obtains
$\rho_{\gamma}=\rho_{G}$, which proves the claim. ∎

_Reduction of the kinetic energy_ We now rewrite the kinetic energy as a
functional of elements of ${\mathcal{G}}$.

######  Theorem 3.4.

Let $G\in{\mathcal{G}}$. Then,

(3.8) $\inf_{\Lambda(\gamma)=G}\underline{{\mathcal{E}}}^{\rm
kin}(\gamma)=\frac{1}{2}{\rm Tr}(-\Delta_{2}G)+\frac{\pi^{2}}{6}{\rm
Tr}(G^{3}).$

###### Proof.

Let $G\in{\mathcal{G}}$ and
$\gamma=\int^{\oplus}_{\mathbb{R}}\gamma_{k}\in{\mathcal{K}}$ with
$\Lambda(\gamma)=G$. Write
$G=\sum_{j}g_{j}|\varphi_{j}\rangle\langle\varphi_{j}|$. We have

$\displaystyle\underline{{\mathcal{E}}}^{\rm kin}(\gamma)$
$\displaystyle=\frac{1}{2}\underline{\rm
Tr}(-\Delta_{3}\gamma)=\frac{1}{4\pi}\int_{\mathbb{R}}{\rm
Tr}((-\Delta_{2}+k^{2})\gamma_{k})\,{\mathrm{d}}k$
$\displaystyle=\frac{1}{2}{\rm
Tr}(-\Delta_{2}G)+\frac{1}{4\pi}\int_{\mathbb{R}}\sum_{j}k^{2}m_{j,k}{\mathrm{d}}k,$

where $m_{j,k}:=\langle\gamma_{k}\varphi_{j},\varphi_{j}\rangle$. Notice that

(3.9) $0\leq m_{j,k}\leq
1\quad\text{and}\quad\displaystyle\int_{\mathbb{R}}m_{j,k}\,{\mathrm{d}}k=2\pi
g_{j}.$

Hence,

$\displaystyle\underline{{\mathcal{E}}}^{\rm kin}(\gamma)\geq\frac{1}{2}{\rm
Tr}(-\Delta_{2}G)+\frac{1}{4\pi}\sum_{j}\inf\left\\{\int_{\mathbb{R}}k^{2}m(k){\mathrm{d}}k,\;m(k)\;\text{satisifying}\;\eqref{eq
bathtube cond}\right\\}.$

According to the Bath-tub principle, [14, Theorem 1.14], the infinimum in the
RHS of the above inequality is achieved by $m_{j}^{*}(k):=1_{\\{|k|<\pi
g_{j}\\}}$, the same state exhibited in the proof of Proposition 3.3.
Therefore,

$\displaystyle\underline{{\mathcal{E}}}^{\rm kin}(\gamma)$
$\displaystyle\geq\frac{1}{2}{\rm Tr}(-\Delta_{2}G)+\frac{1}{6\pi}\sum_{j}(\pi
g_{j})^{3}=\frac{1}{2}{\rm Tr}(-\Delta_{2}G)+\frac{\pi^{2}}{6}{\rm
Tr}(G^{3}),$

with equality for $\gamma^{*}\in{\mathcal{K}}$ with
$\gamma_{k}^{*}=\sum_{j}m_{j}^{*}(k)|\varphi_{j}\rangle\langle\varphi_{j}|$. ∎

######  Remark 3.5.

A consequence of Theorem 3.4 is that for a representable charge density
$\rho\in L^{1}(\mathbb{R}^{2})$, one could associate the kinetic energy
functional as follows

(3.10) ${\mathcal{E}}^{\rm kin}(\rho)=\inf\left\\{{\rm
Tr}(-\Delta_{2}G)+\frac{2\pi^{2}}{3}{\rm
Tr}(G^{3}),\;G\in{\mathcal{G}}_{0},\;\rho_{G}=\rho\right\\}.$

Finding an explicit expression for the above functional
$\rho\mapsto{\mathcal{E}}^{\rm kin}(\rho)$ is highly challenging. It is only
possible to associate a lower bound to this energy as follows

(3.11) $C\left(\int_{\mathbb{R}^{2}}\rho_{G}^{5/3}\right)\leq\frac{1}{2}{\rm
Tr}(-\Delta_{2}G)+\frac{\pi^{2}}{6}{\rm Tr}(G^{3})={\mathcal{E}}^{\rm
kin}(\rho).$

for an appropriate constant $C$. This is a consequence of the Lieb–Thirring
inequality, see for instance [15, 16] and [7, Appendix A] (see Section 5.1).
In particular, this yields that $\rho_{G}\in L^{5/3}$, once the kinetic enegry
of $G$ is finite. Actually, the two dimensional Lieb–Thiring inequality
applied to $G/Z$ implies that $\rho_{G}\in L^{2}(\mathbb{R}^{2})$ and the
Hoffman-Ostenhof inequality claims that $\sqrt{\rho_{G}}\in
H^{1}(\mathbb{R}^{2})$, see [10].

The proof of Theorem 3.1 can be easily deduced from Theorem 3.4 and
Proposition 3.3.

##  4\. Well-posdness of the reduced model

This section is devoted to the analysis of the reduced two-dimensional problem
(3.5). In particular, we prove well-posdness and show that the ground state
satisfies a self-consistent equation.

######  Theorem 4.1.

Assume that $\mu\in L^{2}(\mathbb{R}^{2})\cap L^{1}(\mathbb{R}^{2})$. Then the
energy ${\mathcal{E}}$ in (3.4) admits a unique minimizer $G_{*}$ in
${\mathcal{G}}$. Furthermore, one has

(4.1) $\begin{cases}\quad
G_{*}=\frac{\sqrt{2}}{\,\pi}\left(\lambda_{*}-H_{*}\right)_{+}^{\frac{1}{2}}\\\
\quad H_{*}=-\Delta_{2}+{\mathrm{V}}_{\rho_{*}-\mu}\end{cases},$

where $\rho_{*}:=\rho_{G_{*}}$, $x_{+}:=\max(0,x)$ for $x\in\mathbb{R}$ and
$\left(\lambda_{*}-H_{*}\right)_{+}^{\frac{1}{2}}$ should be unterstood in
terms of functional calculus for the self–adjoint operator $H_{*}$.

###### Proof.

The proof of this theorem follows the same lines as the one of [7, Theorem
2.8]. We detail here the main arguments in our 2D case. _Existence and
uniqueness of the minimizer_. The set ${\mathcal{C}}$ is convex and the
functional ${\mathcal{E}}$ is convex and positive. Then, the existence of a
minimizer is reduced to the lower semi-continuity of ${\mathcal{E}}$. Let let
$(G_{n})_{n}\subset{\mathcal{G}}$ be a minimizing sequence of ${\mathcal{E}}$
and denote by $\rho_{n}:=\rho_{G_{n}}$. One has

${\mathcal{E}}(G_{n})=\frac{1}{2}{\rm
Tr}(-\Delta_{2}G_{n})+\frac{\pi^{2}}{6}{\rm
Tr}(G_{n}^{3})+{\mathcal{D}}(\rho_{n}-\mu)\leq C,$

where $C$ is a positive constant. In particular, $(G_{n})_{n}$ and $(-\Delta
G_{n})_{n}$ are bounded sequences in the Schatten space of trace class
operators ${\mathfrak{S}}(L^{2}(\mathbb{R}^{2}))$ and $(G_{n})_{n}$ is bounded
in the Schatten space of cubic trace class operators
${\mathfrak{S}}_{3}(L^{2}(\mathbb{R}^{2}))$. Hence, taking into the account
that $-\Delta_{2}$ is closed, there exists
$G_{*}\in{\mathfrak{S}}(L^{2}(\mathbb{R}^{2}))\cap{\mathfrak{S}}_{3}(L^{2}(\mathbb{R}^{2}))$
such that

$(G_{n},-\Delta_{2}G_{n})\rightharpoonup(G_{*},-\Delta_{2}G_{*})\quad\text{in}\quad\left({\mathfrak{S}}(L^{2}(\mathbb{R}^{2}))\right)^{2},$
$G_{n}\rightharpoonup
G_{*}\quad\text{in}\quad{\mathfrak{S}}_{3}(L^{2}(\mathbb{R}^{2})),$

and

$\frac{1}{2}{\rm Tr}(-\Delta_{2}G_{*})+\frac{2\pi^{2}}{3}{\rm
Tr}(G_{*}^{3})\leq\liminf_{n}\left(\frac{1}{2}{\rm
Tr}(-\Delta_{2}G_{n})+\frac{2\pi^{2}}{3}{\rm Tr}(G_{n}^{3})\right).$

By the Lieb-Thirring inequality [15], $(\rho_{n})_{n}$ is bounded in
$L^{2}(\mathbb{R}^{2})$. Thus, $(\rho_{n})_{n}$ converges weakly in
$L^{2}(\mathbb{R}^{2})$. Moreover, its weak limit is $\rho_{*}:=\rho_{G_{*}}$.
Indeed, for any function $h\in C_{c}(\mathbb{R}^{2})$ one has

$\displaystyle\int_{\mathbb{R}^{2}}h\rho_{*}{\mathrm{d}}x$ $\displaystyle={\rm
Tr}\left((1-\Delta)^{-1}{\mathcal{M}}_{h}G_{*}(1-\Delta)\right)$
$\displaystyle=\lim_{n}{\rm
Tr}\left((1-\Delta)^{-1}{\mathcal{M}}_{h}G_{n}(1-\Delta)\right)=\lim_{n}\int_{\mathbb{R}^{2}}h\rho_{n},$

where ${\mathcal{M}}_{h}$ stands for the multiplication operator by $h$ in
$L^{2}(\mathbb{R}^{2})$, and the convergence is obtained by the fact that
$(1-\Delta)^{-1}{\mathcal{M}}_{h}$ is a compact operator. Besides,
$\|{\mathrm{W}}_{\rho_{n}-\mu}\|_{2}={\mathcal{D}}(\rho_{n}-\mu)\leq C$, for
all $n\in\mathbb{N}$. Hence,

$W_{\rho_{n}-\mu}\rightharpoonup W_{*}\text{ in
}L^{2}(\mathbb{R}^{2})\quad\text{and}\quad\|W_{*}\|_{2}\leq\liminf_{n}\|{\mathrm{W}}_{\rho_{n}-\mu}\|_{2}.$

Let us prove that $W_{*}=W_{\rho_{*}-\mu}$. Let $\varphi$ be a test function
in the Schwartz space ${\mathcal{S}}(\mathbb{R}^{2})$ such that
$\int_{\mathbb{R}^{2}}\varphi\,{\mathrm{d}}x=0$. By Parseval’s identity, one
gets

$\displaystyle\int_{\mathbb{R}^{2}}{\mathrm{W}}_{\rho_{n}-\mu}(x)\varphi(x){\mathrm{d}}x=\left\langle\frac{{\mathcal{F}}_{2}(\rho_{n}-\mu)}{|k|},{\mathcal{F}}_{2}(\varphi)\right\rangle$
$\displaystyle=\left\langle{\mathcal{F}}_{2}(\rho_{n}-\mu),\frac{{\mathcal{F}}_{2}(\varphi)}{|k|}\right\rangle$
$\displaystyle=\langle\rho_{n}-\mu,{\mathrm{W}}_{\varphi}\rangle.$

Letting $n$ go to $\infty$, and noting that $\rho_{n}-\mu\in
L^{2}(\mathbb{R}^{2})$, we obtain

$\int_{\mathbb{R}^{2}}{\mathrm{W}}_{*}(x)\varphi(x){\mathrm{d}}x=\int_{\mathbb{R}^{2}}\left(\rho_{*}-\mu\right)(x)\,{\mathrm{W}}_{\varphi}(x)\,{\mathrm{d}}x=\int_{\mathbb{R}^{2}}{\mathrm{W}}_{\rho_{*}-\mu}(x)\varphi(x){\mathrm{d}}x.$

This actually identifies ${\mathrm{W}}_{*}$ with
${\mathrm{W}}_{\rho_{*}-\mu}$. Summarizing, one has $G_{*}\in{\mathcal{G}}$
and

$I\leq{\mathcal{E}}(G_{*})\leq\liminf_{n}{\mathcal{E}}(G_{n})=I,$

which proves that $G_{*}$ is actually a minimizer of ${\mathcal{E}}$. The
uniqueness comes from the strict convexity of ${\rm Tr}\left(G^{3}\right)$.

_Proof of the Euler Lagrange equation._ Let $G_{*}$ denote the minimizer of
${\mathcal{E}}$, $\rho_{*}:=\rho_{G_{*}}$ its density and
$V_{*}:=V_{\rho_{*}-\mu}$ the mean-field potential. Let $G\in{\mathcal{G}}$
and $t\in[0,1]$. One has $G_{t}:=G_{*}+t(G-G_{*})\in{\mathcal{G}}$ and
${\mathcal{E}}(G_{t})\geq{\mathcal{E}}(G_{*})$. Besides,

$\displaystyle{\mathcal{D}}(\rho_{G_{t}}-\mu)$
$\displaystyle={\mathcal{D}}(\rho_{*}-\mu)+2t{\mathcal{D}}(\rho_{*}-\mu,\rho_{G}-\rho_{*})+o(t)$
$\displaystyle={\mathcal{D}}(\rho_{*}-\mu)+2\,t\,{\rm
Tr}(V_{*}(G-G_{*}))+o(t).$

Thus,

$t\left({\rm
Tr}\left(-\frac{1}{2}\Delta_{2}(G-G_{*})\right)+\frac{\pi^{2}}{2}{\rm
Tr}\left(G_{*}^{2}(G-G_{*})\right)+{\rm Tr}(V_{*}(G-G_{*}))\right)+o(t)\geq
0.$

Hence, dividing by $t$ and letting $t\to 0$, one obtains

${\rm
Tr}\left(\left(-\frac{1}{2}\Delta_{2}+\frac{\pi^{2}}{2}G_{*}^{2}+V_{*}\right)(G-G_{*})\right)\geq
0.$

Proceeding as in [7, Sect. 4.3.2] and [8, Proposition 3.8], set
$H_{*}:=-\frac{1}{2}\Delta_{2}+V_{*}$ and
$L_{*}=H_{*}+\frac{\pi^{2}}{2}G_{*}^{2}$. Therefore, ${\rm Tr}(L_{*}G)\geq{\rm
Tr}(L_{*}G_{*})$, for all $G\in{\mathcal{G}}$. This implies that $L_{*}$ is
bounded from below and

${\rm Tr}(L_{*}G_{*})=\inf\\{{\rm Tr}(L_{*}G),\;G\in{\mathcal{G}}\\}.$

Set $\lambda_{*}=\inf\sigma(L_{*})$ the bottom of the spectrum of $L_{*}$.
Then, $(L_{*}-\lambda_{*})G_{*}=0$. Next, let us expand $G_{*}$ as follows
$G_{*}=\sum_{j\in\mathbb{N}}g_{j}^{*}|\phi_{j}^{*}\rangle\langle\varphi_{j}^{*}|$,
where $g_{j}\geq 0$ for any $j\in\mathbb{N}$ with $\sum_{j}g_{j}=Z$, and
$\overline{\\{\varphi_{j},\;j\in\mathbb{N}\\}}=L^{2}(\mathbb{R}^{2})$. One has
$(L_{*}-\lambda_{*})\varphi_{j}^{*}=0$, for any $j\in\mathbb{N}$ such that
$g_{j}>0$. That is

$H_{*}\varphi_{j}^{*}=\left(-\frac{1}{2}\Delta_{2}+V_{*}\right)\varphi_{j}^{*}=\left(\lambda_{*}-\frac{1}{2}(\pi
g_{j}^{*})^{2}\right)\varphi_{j}^{*},\quad\text{once}\quad g_{j}>0.$

Hence, $s_{j}:=\lambda_{*}-\frac{1}{2}(\pi g_{j}^{*})^{2}$ belongs to the
point spectrum of $H_{*}$. Conversely, if $s\in\sigma(H_{*})$ such that
$s<\lambda_{*}$, then $g=\frac{\sqrt{2}}{\pi}\sqrt{\lambda_{*}-s}$ is an
eigenvalue of $G_{*}$. This yields
$G_{*}=\frac{\sqrt{2}}{\pi}\left(\lambda_{*}-H_{*}\right)_{+}^{1/2}$. Note
that $H_{*}$ is a compact perturbation of $L_{*}$; thus they share the same
essential spectrum. It follows that $H_{*}$ cannot have essential spectrum
below $\lambda_{*}$. ∎

##  5\. Numerical simulation

###  5.1. Thomas-Fermi model

We perform numerical illustrations on the simple Thomas-Fermi model, which can
be seen as a semi–classical limit of the rHF model. The kinetic energy in this
model reads

${\mathcal{E}}^{\rm TF}_{\rm kin}(\rho)=c_{\rm
TF}\int_{\mathbb{R}^{2}}\rho^{5/3}(x){\mathrm{d}}x,\quad\text{with}\quad
c_{{\rm TF}}=\frac{3}{10}(3\pi^{2})^{2/3}.$

Therefore, for a given $0\leq\mu\in L^{1}(\mathbb{R}^{2})$, we define the
Thomas-Fermi energy per unit length as

${\mathcal{E}}^{\rm TF}(\rho):=c_{\rm
TF}\int_{\mathbb{R}^{2}}\rho^{\frac{5}{3}}(x){\mathrm{d}}x+\frac{1}{2}{\mathcal{D}}(\rho-\mu),$

for every $\rho\in{\mathcal{C}}^{\rm TF}$, where

${\mathcal{C}}^{\rm TF}:=\\{\rho\in L^{1}(\mathbb{R}^{2})\cap
L^{5/3}(\mathbb{R}^{2}),\;\rho-\mu\in{\mathcal{C}}\\}.$

We point out that Thomas–Fermi model and its derivatives are widely considered
in the literature, see [17, 12, 2]. The functional ${\mathcal{E}}^{\rm TF}$ is
stricltly convex and has a unique minimizer $\rho_{\rm TF}$ that satisifes the
following self–consistent equation

(5.1)
$\begin{cases}\rho_{TF}=\left(\frac{3}{5c_{TF}}(\lambda-{\mathrm{V}}_{{\rm
TF}})_{+}\right)^{3/2}\\\ {\mathrm{V}}_{\rm TF}={\mathrm{V}}_{\rho_{\rm
TF}-\mu}\end{cases},$

with $\lambda$ the Fermi level of the system is such that
$\int_{\mathbb{R}^{2}}\rho_{\rm TF}=Z:=\int_{\mathbb{R}^{2}}\mu$.
We emphasize that we still use our expression for the potential energy
${\mathcal{D}}$ and the mean–field potential in (2.7) which is involved in the
self–consistent equation (5.1). We next show a numerical solution to (5.1). We
also contrast the outcomes from this equation with the non-regularized reduced
energy in which the Coulomb energy and potential are

(5.2)
${{\mathcal{D}}}_{2}(f,g)=-2\int_{\mathbb{R}^{2}}\int_{\mathbb{R}^{2}}\log(|x-y|)f(x)g(y)\,{\mathrm{d}}x{\mathrm{d}}y$

and

(5.3)
$\widetilde{{\mathrm{V}}}_{f}(x)=-2\int_{\mathbb{R}^{2}}\log(|x-y|)f(y)\,{\mathrm{d}}y.$

###  5.2. Numerical schemes

Python is used to carry out the numerical simulation. The two-dimensional
functions are evaluated on a two-dimensional grid representing
$[-a,a]\times[-b,b]$, with $N_{a}\times N_{b}$ discretization points (we took
$a=b=8$ and $N_{a}=N_{b}=41$). The various integrals inolved in the total
energy are then approximated using quadrature techniques; we used the two-
dimensional trapezoid rule.

We use two different methods in the regularized and non-regularized case. In
the latter, the total energy is approximated using a numerical integration
method. The constraints of the model are simply the positivity of the density
$\rho$ along with the neutrality condition $\int_{\mathbb{R}^{2}}\rho=Z$. The
resulting discretized problem is thus a finite-dimensional optimization
problem on $\rho_{ij}$ (the approximated values of the density on the grid
points) subject to linear constraints. This problem is solved numerically in
Python using the predefined function minimize from Scipy library. In the
regularized case, we apply an iterative process to solve equation (5.1) as
follow: at each iteration $n$, we evaluate the potential
${\mathrm{V}}_{n}:={\mathrm{V}}_{\rho_{n}-\mu}$ using (2.7) and we compute
$\lambda_{n}$ solution to

(5.4) $Z_{n}(\lambda_{n})=0,\quad\text{where}\quad
Z_{n}(\lambda)=\int_{\mathbb{R}^{2}}\left(\frac{3}{5c_{TF}}(\lambda-\Phi_{n})_{+}\right)^{3/2}-Z.$

Note that, since $Z_{n}$ is a non-decreasing function of $\lambda$, a simple
dichotomy can be used to compute efficiently $\lambda_{n}$. We then update the
density by setting

(5.5)
$\widetilde{\rho_{n+1}}=\left(\frac{3}{5c_{TF}}(\lambda_{n}-\Phi_{n})_{+}\right)^{3/2}\quad\text{and}\quad\rho_{n+1}=t\widetilde{\rho_{n+1}}+(1-t)\rho_{n},$

where $t\in[0,1]$ is optimized to minimize the Thomas-Fermi energy (we use
again the Python predefined function minimize). The algorithm is terminated
when a tolerance $\varepsilon$ is reached between two successively evaluated
energies. We used $\varepsilon=10^{-5}$.

For the approximation of the regularized energy $\mathcal{D}(\rho-\mu)$ by
trapezoid rule, we use the formula (2.8) involving the potential, which is
valid since $\rho-\mu$ is neutral, instead of the one involving the Fourier
transform.

###  5.3. Numerical results

We present here the numerical results obtained for our test case, which
corresponds to a homogeneous charged nanowire localized in a given square

(5.6) $\mu_{1}(x)=\mathds{1}(|x_{1}|<2)\mathds{1}(|x_{2}|<2).$

We note that we take $c_{TF}=\frac{2}{3}\pi^{2}$. We start by plotting
$\mu_{1}$ in Figure 1. Then, Figures 2 and 3 show the results for the non-
regularized and regularized models, respectively. The two models’ combined
results are fairly close. Figure 4 shows the error between the two models. We
are able to say the following. The difference between the two densities
appears to be almost zero far from the nanowire. The relative error is of
order $13\%$ between potentials and of order $11\%$ between densities in the
material’s vicinity. Ultimately, there is a $7\%$ error between the two ground
states energies in the two models.

(a) 2D representation.

(b) 3D representation.

Figure 1. Nuclear density $\mu_{1}$ of a simple nanowire.

(a) Potential $\tilde{{\mathrm{V}}}$

(b) Density $\tilde{\rho}$

Figure 2. Results of the non-regularized model for $\mu_{1}$.

(a) Regularized potential ${\mathrm{V}}_{{\rm TF}}-\lambda$

(b) Regularized density $\rho_{{\rm TF}}$

Figure 3. Results of the regularized model for $\mu_{1}$.

(a) $|{\mathrm{V}}_{TF}-\lambda-\tilde{{\mathrm{V}}}|$

(b) $|\rho_{TF}-\tilde{\rho}|$

Figure 4. Errors between the regularized and non-regularized models for
$\mu_{1}$.

From a computational view point, the regularized model ends in 1902s with 13
iterations versus 2025s for the non-regularized model. Moreover, the energy
obtained by the regularized model is lower than that of the non-regularized
model (83,11 versus 89,46). More in-depth study remains to be done with finer
grids and other test cases, but this first simple test allows to see the
interest of the regularized model with respect to the non-regularized model,
where there are no theoretical results for the existence and uniqueness of the
minimizer.

## References

  * [1] A. A. Balandin, F. Kargar, T. T. Salguero, and R. K. Lake, One-dimensional van der waals quantum materials, Materials Today, 55 (2022), pp. 74–91.
  * [2] X. Blanc and C. Bris, Thomas-Fermi type theories for polymers and thin films, Advances in Differential Equations, 5 (2000).
  * [3] V. Chaudhary, P. Neugebauer, O. Mounkachi, S. Lahbabi, and A. EL Fatimy, Phosphorene - an emerging two-dimensional material: recent advances in synthesis, functionalization, and applications, 2D Materials, 9 (2022).
  * [4] Y. Dai, W. Wei, Y. Ma, and C. Niu, Calculations and Simulations of Low-Dimensional Materials: Tailoring Properties for Applications, John Wiley & Sons, 2022.
  * [5] L. Frerick, J. Müller, and T. Thomaser, A Fourier integral formula for logarithmic energy, ArXiv preprint, https://arxiv.org/abs/2209.05439, (2022).
  * [6] A. Geim and I. Grigorieva, Van der Waals heterostructures, Nature, 499 (2013), pp. 419–425.
  * [7] D. Gontier, S. Lahbabi, and A. Maichine, Density functional theory for two-dimensional homogeneous materials, Communication in Mathematical Physics, 388 (2021), pp. 1475–1505.
  * [8] D. Gontier, S. Lahbabi, and A. Maichine, Density functional theory for two-dimensional homogeneous materials with magnetic fields, Journal of Functional Analysis, 285 (2023), p. 110100.
  * [9] D. Gupta, V. Chauhan, and R. Kumar, A comprehensive review on synthesis and applications of molybdenum disulfide (mos2) material: Past and recent developments, Inorganic Chemistry Communications, 121 (2020), p. 108200.
  * [10] M. Hoffmann-Ostenhof and T. Hoffmann-Ostenhof, "Schrödinger inequalities" and asymptotic behavior of the electron density of atoms and molecules, Phys. Rev. A, 16 (1977), pp. 1782–1785.
  * [11] F. Hussain, D. Muhammad Imran, and H. U. Janjua, Density Functional Theory (DFT) Study of Novel 2D and 3D Materials, 07 2017, pp. 269–284.
  * [12] P.-L. L. Isabelle Catto, Claude Le Bris, The Mathematical Theory of Thermodynamic Limits: Thomas–Fermi Type Models, Oxford Mathematical Monographs, Oxford University Press, USA, 1998.
  * [13] T. Kurokawa, Riesz potentials, higher Riesz transforms and Beppo Levi spaces, Hiroshima Mathematical Journal, 18 (1988), pp. 541 – 597.
  * [14] E. Lieb and Loss, Analysis, vol. 14 of Graduate studies in mathematics, American Mathematical Society, 2nd ed., 2001.
  * [15] E. Lieb and W. Thirring, Bound on kinetic energy of fermions which proves stability of matter, Phys. Rev. Lett., 35 (1975), pp. 687–689.
  * [16] , Inequalities for the moments of the eigenvalues of the Schrödinger hamiltonian and their relation to Sobolev inequalities, Studies in Mathematical Physics, Princeton University Press, 1976, pp. 269–303.
  * [17] E. H. Lieb and B. Simon, The Thomas-Fermi theory of atoms, molecules and solids, Advances in Mathematics, 23 (1977), pp. 22–116.
  * [18] H. Moustafa, Computational properties of one-dimensional materials, (2023).
  * [19] Paras, K. Yadav, P. Kumar, D. R. Teja, S. Chakraborty, M. Chakraborty, S. S. Mohapatra, A. Sahoo, M. M. Chou, C.-T. Liang, et al., A review on low-dimensional nanomaterials: Nanofabrication, characterization and applications, Nanomaterials, 13 (2022), p. 160.
  * [20] A. Patra, S. Jana, P. Samal, F. Tran, L. Kalantari, J. Doumont, and P. Blaha, Efficient band structure calculation of two-dimensional materials from semilocal density functionals, The Journal of Physical Chemistry C, 125,20 (2021), pp. 11206–11215.
  * [21] M. Reed and B. Simon, Methods of Modern Mathematical Physics. Analysis of Operators, vol. IV, Academic Press, 1978.
  * [22] J. Saha and A. Dutta, A review of graphene: Material synthesis from biomass sources, Waste and Biomass Valorization, 13 (2022).
  * [23] E. M. Stein, Singular Integrals and Differentiability Properties of Functions, vol. 30 of Princeton Mathematical Series, Princeton University Press, 1970.
  * [24] Q. Tang, Z. Zhou, and Z. Chen, Innovation and discovery of graphene-like materials via density-functional theory computations, WIREs Comput Mol Sci, 5 (2015).
  * [25] H. Toffoli, S. Erkoç, and D. Toffoli, Modeling of nanostructures, (2012).
  * [26] Z. Wang, T. Hu, R. Liang, and M. Wei, Application of zero-dimensional nanomaterials in biosensing, Frontiers in chemistry, 8 (2020), p. 320.


# Machine Learning Techniques for Source Localisation in Elastic Media

Bansi Mandalia1, Steve Greenwald1, Simon Shaw2, Gregory Slabaugh1
1Queen Mary University of London
2Brunel University London

###### Abstract

Coronary Artery Disease (CAD) results from plaque deposit in a coronary
artery. Early diagnosis is imperative, so a non-invasive detection method is
being developed to identify acoustic signals caused by partial occlusions in
the artery. The blood flow in the artery is disturbed and imposes oscillatory
stresses on the artery wall. The deformations caused by the stresses can be
detected at the chest surface. Therefore, by using data simulating these
surface signals, which arise from randomly assigned source positions, machine
learning (ML) can be utilised to predict the source of the occlusion. Seven ML
algorithms were investigated, and the results from this study found that an
ensemble model combining k-Nearest Neighbours and Random Forest had the best
performance. The metrics used to evaluate this was the mean squared error and
Euclidean distance.

###### Index Terms:

machine learning, regression, non-invasive diagnosis, coronary artery disease,
source localisation

##  I Introduction

Coronary Artery Disease (CAD) results from plaque deposit in a coronary
artery. This reduces the amount of oxygenated blood from reaching the heart,
thus early diagnosis is imperative in treating the disease and preventing
death. Hence, non-invasive, accurate detection methods are of great interest.
One of these methods identifies acoustic signals caused by partial occlusions,
as the blood flow in the artery is disturbed and causes oscillatory stresses
on the artery wall. The deformations caused by the stresses can be detected at
the chest surface. Figure 1 is a diagram showing this concept.

Figure 1: Basic diagram representing a person with a partial occlusion in
their heart (CAD) and acoustic signals emanating from this to the detector
surface (Clipart Library 2019).

By using data gathered from simulating these surface signals, machine learning
(ML) can be utilised to predict the source of the occlusion. This can be used
by health professionals for diagnosis and more effective treatment.

Preliminary experiments were conducted and it was found that the method used
was time consuming. Hence, the investigation was limited and this motivates
the use of ML in this report. ML is a set of tools for solving problems by
understanding complex datasets and can learn and improve automatically.

The type of problem this report focuses on is a regression problem. ML
algorithms can be used to predict the value of a continuous attribute. In this
case, we wish to predict the value of the coordinates of the source
disturbance ($x_{c},y_{c},z_{c}$) using the sensor readings
($u_{1},u_{2},u_{3}$). Thus, several ML algorithms were explored: Linear
Regression (LR), XGBoost, Decision Tree (DT), Neural Network (NN), k-Nearest
Neighbours (kNN), Random Forest (RF) and an ensemble model.

The research challenge that this project aims to solve is: given some
generated surface sensor readings can a ML model predict the source of this
signal, and therefore detect the position of the occlusion? The overall aim of
this project is to investigate whether ML can be used to solve this problem,
as there is currently no existing system that does this. Furthermore, this
project aims to find the best performing ML model for this task through
experimenting with different algorithms and tuning methods.

##  II Related Work

Several studies were reviewed and analysed in order to better understand
current methods and limitations of CAD diagnosis using ML. The earliest study
from Hesser et al. (2020) developed ML algorithms in order to predict the
location of an active source using the response from piezoelectric sensors.
The data used for this was gathered by dropping a steel ball on an aluminium
plate at different positions, and measuring the resulting elastic waves
through sensors. Furthermore, numerical data was also gathered by performing
500 simulations of this experiment.

Both the numerical and experimental data was used to train an Artificial
Neural Network (ANN) and a Support Vector Machine (SVM) in order to predict
the impact location of the steel ball based on the wave response of each
sensor. The results showed that ML algorithms can precisely identify the
location of a source. One evaluation method they used was the distance, in the
$x$ and $y$ direction, between the source and the prediction of the ANN and
SVM trained with both numerical and experimental data.

The limitations of this study are that only two types of ML algorithms were
investigated and a small dataset was used to train these models. The
difference between the paper from Hesser et al. (2020) and this study is that
the experimental setup only generated 2D numerical data. Specifically, the
study reviewed was only able to locate the $x$ and $y$ coordinates of source
whereas, in this study the $z$ coordinate will also be located. Furthermore,
the study reviewed used a steel ball and a thin aluminium plate. However, in
this study a generated occlusion embedded in a rectangular hexahedral domain
will be used. Overall, this paper shows that it is possible to use ML
algorithms to locate the coordinates of a source using data from sensors.

The second study from Nasarian et al. (2020) used a dataset containing
workplace and environmental features to train a ML algorithm. The main focus
of this paper is the proposal of a new feature selection algorithm called
heterogeneous hybrid feature selection (2HFS). The dataset used was the
Nasarian CAD dataset and some workplace features included are office location
of patients, noise exposure and pollutants. They also tested their method on
other CAD datasets such as the Z-Alizadeh Sani dataset.

The ML algorithms which were investigated are a Decision Tree (DT), Gaussian
Naïve Bayes, Random Forest (RF), XGBoost, k-Nearest Neighbours (kNN) and
Bernoulli Naïve Bayes. The results showed that that these workplace and
environmental features have an impact on the detection of CAD. Specifically,
the best ML algorithm in this study, XGBoost, achieved an accuracy of 81.23%.

The limitations of this study is that no neural networks (NN) were
investigated and all the subjects in the Nasarian dataset were male, thus the
results are not fully representative of the population. The difference between
the paper from Nasarian et al. (2020) and this study is the type of data used.
The study reviewed used categorical data whereas, in this paper numerical data
in terms of signals will be used to train ML algorithms. Furthermore, the
models in the reviewed study predicted whether a patient was likely to have
CAD whereas, the model in this paper will predict the source of the occlusion
which causes CAD. Overall, this study showed that feature selection during
model training is important and has an effect on accuracy.

The third study from Akella et al. (2021) investigated different ML algorithms
in order to predict the likelihood of a person having CAD. Unlike previous
studies beforehand, this study made their code used to build these ML models
available publicly, so it can be used to assist healthcare workers in
diagnosing CAD. The dataset utilised in this study is ‘the Cleveland dataset’
which consists of 14 variables obtained from 303 patients. Examples of
variables in the dataset include age, blood cholesterol and chest pain.

The ML algorithms they investigated were LR, regression tree, RF, SVM, kNN and
NN. The results shows that NNs achieved the highest accuracy, 93%, in
predicting the presence of CAD. The limitation of this study is that the
dataset used is quite small and quite old, it is from 1988, therefore the data
may not be representative of the population today. Furthermore, some of the ML
algorithms utilised are very similar to each other, e.g. DT and RF, therefore
it would have been interesting to explore more varied algorithms and compare
their results. Similarly to the previous paper, the difference between the
paper from Akella et al. (2021) and this study is the type of data used and
what the ML model predicts.

The fourth paper from Doppala et al. (2022) aimed to improve the accuracy of
ML methods for cardiovascular disease diagnosis by utilising an ensemble
model. The data was obtained from the Mendeley Data Centre, IEEE DataPort and
Cleveland dataset. This paper is different from previous studies as it
combines four algorithms in order to create an ensemble model which achieves a
higher accuracy than previous ML models. The ensemble model was a combination
of naïve Bayes, RF, SVM and XGBoost. Each classifier made a prediction then, a
‘Voting Classifier’ selected the most popular outcome.

The results showed that the proposed ensemble model achieves a greater
accuracy than the benchmark classifiers (DT, RF, Naïve Bayes, Logistic
Regression, SVM, Gradient Boosting, XGBoost). The limitations of this study
are that NNs were not investigated. The difference between the paper from
Doppala et al. (2022) and this study is that the dataset of the reviewed study
consisted of categorical data such as age, gender, resting ECG etc. whereas,
in this paper numerical data in terms of signals will be used to train a ML
model. Overall, this paper showed how a voting classifier can be used in order
to combine several different models together, which improves the accuracy of
CAD diagnosis compared to benchmark classifiers.

##  III Methodology

A computational model that simulates the sensor readings as a result of a
source disturbance in the coronary artery will be utilised. The source
disturbance is centred at a point $\boldsymbol{x}_{c}=(x_{c},y_{c},z_{c})$
inside a rectangular hexahedral domain. It is of the form:

$\boldsymbol{F}(\boldsymbol{x}):=A(1,1,1)^{T}\exp\left(-\frac{\|\boldsymbol{x}-\boldsymbol{x}_{c}\|_{\mathbb{E}}^{2}}{\epsilon}\right)$
(1)

where $A$ is a constant, and $\epsilon$ is small and positive to control the
‘local-ness’ (Shaw and Greenwald 2022).

The body has a dimensions of $[0.0,0.3]\times[-0.05,0.05]\times[0.0,0.05]$, in
metres, along the $x$, $y$ and $z$ axes. The sensor readings measure the
static displacement field as a response to this disturbance. The displacements
$\boldsymbol{u}=(u_{1},u_{2},u_{3})^{T}$ and other quantities, such as
derivatives, are measured at sites on the surface using either five
microphones or four accelerometers. The coordinates of the five microphones
are: $(x,y)$ coordinates: $(0.12,0.01)$, $(0.18,0.01)$, $(0.15,0.00)$,
$(0.12,-0.01)$, $(0.18,-0.01)$. The coordinates of the four accelerometers
are: $(x,y)$ coordinates: $(0.15,0.01)$, $(0.12,0.00)$, $(0.18,0.00)$ and
$(0.15,-0.01)$. These sensors are located on the top face at $z=0.05$. A
schematic diagram of the experimental setup can be seen in Figure 2. In the
simulation there is no intrinsic difference between the simulated microphone
and accelerometer data, other than the site at which it is collected.

Figure 2: Basic diagrams of the experimental set up. The markers in the top
diagram represent the locations of the microphones. The markers in the bottom
diagram represent the locations of the accelerometers.

These sensor readings will then be used to train various ML algorithms to
predict the source of the disturbance. The system should simulate data which
is similar to real-life results, as the purpose of this research is to later
use the ML model on patient data readings. Another requirement of the system
is to generate many readings. This is because, ML algorithms require a large
amount of data for training in order to generate good results when testing.
Furthermore, the ML algorithms should perform just as well during deployment
on unseen data as during training.

The data obtained contains $u_{1}$, $u_{2}$ and $u_{3}$ which represent the
displacements in the $x$, $y$ and $z$ directions. Furthermore for each of
these displacements, $\boldsymbol{u}_{x}$, $\boldsymbol{u}_{y}$ and
$\boldsymbol{u}_{z}$ represent the derivative with respect to the $x$, $y$ and
$z$ direction. These measurements are obtained for each of the 5000 randomly
generated points $\boldsymbol{x}_{c}$. Also, these sets are replicated for two
different narrowness of the source $\epsilon$ = 0.01 and $\epsilon$ = 0.001,
for two different refinements of the finite element mesh 10 $\times$ 5
$\times$ 4 mesh and 20 $\times$ 10 $\times$ 8 mesh, and lastly for two
different approximating polynomial degrees $r$ = 1 and $r$ = 2. Hence there
are eight datasets with 5000 data points in each (Shaw and Greenwald 2022). In
addition, for each run a randomly chosen source location,
$\boldsymbol{x}_{c}=(x_{c},y_{c},z_{c})$, was also obtained.

Some pre-processing steps were taken before applying the ML algorithms.
Firstly, for this investigation the three datasets $u_{1}$, $u_{2}$ and
$u_{3}$ were combined with their derivatives into one singular dataset. This
data was then split into a training and validation dataset, where 70% of the
data was used for training and 30% was used for validation in order to tune
the models. However, when implementing the NN the validation dataset is split
again to create a testing dataset. So, the split for the NN is 70% for
training, 20% for validation and 10% for testing. This is because, NNs use the
validation dataset while training therefore, an unseen dataset is needed when
making predictions to avoid overfitting. After this, the singular dataset
containing the $u_{1}$, $u_{2}$ and $u_{3}$ displacements and derivatives were
normalised using the following equation:

$z=\frac{x-\mu}{\sigma}$ (2)

where $\mu$ is the mean and $\sigma$ is the standard deviation. This is done
to ensure that inputs will be treated equally and have a similar range of
amplitudes.

The algorithms which will be investigated are: LR, DT, RF, XGBoost, kNN, NN
and an ensemble model. Linear regression (LR) is one of the most simple ML
algorithms. This is because it assumes that the underlying relationship
between the independent and dependent variables is linear. Specifically for
this study, multiple linear regression will be used as there are three
predictors (the $x$, $y$ and $z$ coordinates of the source disturbance). The
‘LinearRegression’ module from Scikit-Learn will be utilised (Pedregosa et al.
2011) and there are no parameters to be tuned. Overall, this algorithm will be
implemented as a baseline in order to compare to more complex models.

There are two key attributes of a Decision Trees (DT) which are nodes and
branches. A node represents a subdivision of the DT and there are three types:
the root node, internal node and leaf nodes. The root node is at depth 0 and
represents a choice the algorithm must make which results in two mutually
exclusive subsets of the DT (Song et al. 2015). An internal node represents
splits in the tree between different choices. In a DT for regression these
choices are based on values. Lastly, the leaf nodes represent the final result
of the DT. Branches represent the pathway between different types of nodes
based on the outcomes of the choices made by the nodes.

The ‘DecisionTreeRegressor’ module from Scikit-Learn will be utilised
(Pedregosa et al. 2011). This algorithm is trained using the structure of a
tree, as described before, to predict a continuous output of the coordinates
of the source disturbance. The feature which will be tuned in this model is
‘max_depth’ which is the number of layers in the DT.

The Random Forest (RF) model is an ensemble of decision trees. The combination
of these DTs for regression results in a very powerful model. This algorithm
randomises the training samples in each instance of the decision tree, in
order to reduce the risk of overfitting. The output of this model is the
average prediction of each individual tree. The ‘RandomForestRegressor’ module
from Scikit-Learn will be utilised (Pedregosa et al. 2011). The features which
will be tuned in this model are ‘n_estimators’ which is the maximum number of
estimators (trees) used to train and ‘max_depth’.

XGBoost stands for Extreme Gradient Boosting. The reason for using a “Boosted”
decision tree for regression is because, a DT is known as a weak learner which
is a classifier that performs only slightly better than random guessing, and
they can be easily over-trained (Géron 2017). However, ‘boosting’ makes the DT
become a strong learner. A stronger learner is a classifier which achieves
high accuracy (Géron 2017). Boosting is an example of a re-weighting technique
where mis-classified events are re-weighted from one iteration of training to
the next. Therefore, subsequent training iterations focus on correcting mis-
classified events. The ‘XGBRegressor’ class will be utilised (Chen et al.
2016). The features which will be tuned in this model are ‘n_estimators’,
‘max_depth’ and ‘learning_rate’ which is the step size for the gradient
boosting algorithm.

The k-Nearest Neighbours (kNN) model is a non-parametric approach that does
not have any explicit decision boundaries, rather a sample is assigned the
label of the most similar training sample. It does this by calculating the
Euclidean distance to all the training samples, then selects the ‘k’ closest
samples (neighbours) and assigns them a label. The ‘KNeighborsRegressor’
module from Scikit-Learn will be utilised (Pedregosa et al. 2011). The
features which will be tuned in this model are ‘n_neighbors’, which represents
the number of neighbours, ‘k’. Also, the parameter ‘weights’, which affects
the influence samples have in each neighbourhood.

The neural network (NN) in this report consists of single layers connected
sequentially. The data passes through the input layer, hidden layers and
lastly an output layer. Each perceptron within a layer is defined by weights,
activation functions and biases. NNs also have a loss function which measures
the performance of the NN by comparing the desired output to the actual
output, and calculates the output error (Géron 2017). For this investigation,
the ‘Mean Squared Error‘ loss function will be used. TensorFlow specifically
with Keras will also be used to create a NN in Python.

An ensemble model is a combination of different ML algorithms with a ‘Voting
Regressor’ and is the approach the study from Doppala et al. (2022) also
investigated. They found that that accuracy of their model improved using this
method. A combination of all models was tested however, the ensemble model
performed better when using just RF and kNN to train. This is because,
ensemble learning works best when training with a diverse range of models, as
this varies the types of errors and thus improves the ensemble’s accuracy
(Géron 2017). The purpose of the ‘Voting Regressor’ is to average the
individual predictions to output a final prediction. The ‘VotingRegressor‘
module from Scikit-Learn will be utilised (Pedregosa et al. 2011). The
parameters of the RF and kNN will be tuned separately, and the best values for
each parameter will be used in the ensemble model.

To tune these models a module from Scikit-Learn will be utilised called
‘GridSearchCV’. This technique searches for the optimal value to use for each
specified parameter. Therefore, instead of manually tuning each model,
‘GridSearchCV’ does this automatically. The input for this function are the
parameters that require tuning for each model. This is typically
‘n_estimators’, ‘max_depth’ and ‘learning_rate’. The algorithm then computes
all the possible combinations of these values, and comes to a conclusion on
which is the best combination using cross-validation (Géron 2017).

In order to implement all of these different ML algorithms, Python will be
used with Jupyter Notebooks. The evaluation techniques used will be the Mean
Squared Error (MSE):

$E_{MSE}=\frac{1}{N}\sum^{N}_{i=1}(y_{i}-\hat{y}_{i})^{2}$ (3)

where, $y_{i}$ is the continuous true label of the sample $i$, $\hat{y}_{i}$
is the predicted label for the sample and $N$ is the number of samples in the
dataset. Furthermore, the average distance between the prediction on the
validation dataset and the ground truth of the source disturbance will also be
calculated using the Euclidean distance:

$d(y_{i},\hat{y}_{i})=\sqrt{\sum^{N}_{i=1}(y_{i}-\hat{y}_{i})^{2}}$ (4)

where, $y_{i}$ and $\hat{y}_{i}$ represent the ground truth and prediction
respectively and $N$ represents the number of samples. This is a similar
evaluation metric as used in Hesser et al. (2020). Lastly, 3D plots were
created with the same geometry as the hexahedral domain to visualise the
difference between the ground truth and prediction of the source disturbance.

These evaluation measures will also be used to assess whether the model is
being overfitted or underfitted. The model will display overfitting if it has
a very low training MSE but a high validation MSE. This occurs when the model
is memorising irrelevant details of the training samples, and thus does not
perform well during validation. A model will display underfitting if the
training and validation MSEs are high. In this case, the model is unable to
identify the underlying pattern within the data.

##  IV Results

Several tests were conducted in order to investigate which ML technique
achieved the best results. This included using different subsets of the data
depending on the mesh resolution, $\epsilon$ value, polynomial degree and
whether it contained microphone data, accelerometer data or both. Firstly,
only one single file corresponding to $\epsilon=0.01$, polynomial degree $r=1$
and 10 $\times$ 5 $\times$ 4 mesh and only microphone data was used to train
each model.

Table I shows the validation MSEs after training and tuning the models as well
as the respective MSEs of the $x$, $y$ and $z$ coordinates. The validation MSE
was calculated using equation 3, where $y_{i}$ is from the validation dataset
and $\hat{y}_{i}$ is the models’ prediction on the validation dataset. The
MSEs of the $x$, $y$ and $z$ coordinates was calculated using a similar
methodology however, the prediction for each individual coordinate was used so
three different MSEs for each coordinate was computed.

Model |  Validation MSE ($\times 10^{-5}$) |  MSE ($\times 10^{-5}$)
---|---|---
|  | x | y | z
Linear Regression | 27.1 | 66.1 | 9.66 | 5.53
XGBoost | 6.77 | 8.02 | 3.73 | 8.56
Decision Tree | 5.86 | 8.10 | 2.49 | 6.99
Neural Network | 3.14 | 1450 | 172 | 46.6
kNN | 1.62 | 1.01 | 0.57 | 3.27
Random Forest | 2.53 | 3.72 | 0.72 | 3.15
Ensemble | 1.69 | 1.85 | 0.37 | 2.84
TABLE I: The validation MSEs and the MSEs of the $x$, $y$ and $z$ coordinates
for the different models using only one single file corresponding to
$\epsilon=0.01$, polynomial degree $r=1$ and 10 $\times$ 5 $\times$ 4 mesh and
only microphone data.

These results show that the three best performing models are kNN, RF and an
ensemble model, as both their validation MSE and the MSEs of each coordinate
are lower than the other models. To further investigate this, the average
distance between the prediction and the ground truth was plotted for all
coordinates and models. This can be seen in Figure 3.

Figure 3: The average distance between the prediction on the validation
dataset and the ground truth for all models and coordinates. Only one single
file corresponding to $\epsilon=0.01$, polynomial degree $r=1$ and 10 $\times$
5 $\times$ 4 mesh and only microphone data was used.

This figure shows that the ensemble model was the most accurate in predicting
the coordinates of the source disturbance. In addition, the results show that
the $y$-coordinate of this disturbance was the best predicted coordinate for
all models except LR. Whereas, the $z$-coordinate was the most difficult
coordinate for the kNN, RF and ensemble model to predict.

Figure 4: A 3D plot comparing the ground truth to the prediction for the
ensemble model for a single example. Only one file corresponding to
$\epsilon=0.01$, polynomial degree $r=1$ and 10 $\times$ 5 $\times$ 4 mesh and
only microphone data was used. Figure is not to scale.

A 3D plot showing the truth in comparison to the predicted values for the
ensemble model can be seen in Figure 4. This plot shows that the prediction
made by the ensemble model is near to the ground truth, with only a small
deviation. Specifically, the coordinate with the largest deviation from the
truth and prediction is the $z$-coordinate, which is averaged to be 3.7mm. In
comparison, the $y$-coordinate has the smallest average difference between the
truth and prediction of 1.4mm.

Similarly to before, one single file corresponding to $\epsilon=0.01$
polynomial degree $r=1$ and 10 $\times$ 5 $\times$ 4 mesh was used however,
only accelerometer data was used to train the models. There was no major
difference in the results when comparing microphone and accelerometer data. A
comparison between the microphone and accelerometer data for this
configuration, as well as as for $\epsilon=0.001$ polynomial degree $r=1$ and
10 $\times$ 5 $\times$ 4 mesh can be seen in Figure 5.

Figure 5: The average distance between the prediction using the ensemble model
and the ground truth, for each coordinate. Using only one file corresponding
to polynomial degree $r=1$ and 10 $\times$ 5 $\times$ 4 mesh for
$\epsilon=0.01$ and $\epsilon=0.001$, with microphone and accelerometer data.

It shows that the $x$-coordinate prediction was slightly better, and the
$y$-coordinate prediction was slightly worse for $\epsilon=0.001$.
Furthermore, it shows that for both these configurations, the $y$-coordinate
error was the smallest and the $z$-coordinate error was the largest. This
shows that the narrowness of the source ($\epsilon$) is almost irrelevant in
predicting the location of the source disturbance.

Secondly, two models were created using all data either corresponding to
$\epsilon=0.01$ or $\epsilon=0.001$. This contained data using polynomial
degree $r=1$ and $r=2$, as well as 10 $\times$ 5 $\times$ 4 mesh and 20
$\times$ 10 $\times$ 8 mesh. Furthermore, two subsets using this data were
carried out, first using only microphone data and then using only
accelerometer data.

The parameters used for the RF model were, ‘n estimators’: 800 and ‘max
depth’: 25. The parameters used for the kNN model were, ‘n estimators’: 4 and
‘weights’: ‘distance’. Table II shows the MSEs of the ensemble model trained
using data corresponding to $\epsilon=0.01$ and $\epsilon=0.001$ with
microphone or accelerometer data. Comparing these results to the first test,
where only one file corresponding to $\epsilon=0.01$, polynomial degree $r=1$
and 10 $\times$ 5 $\times$ 4 mesh was used, the $z$-coordinate error has
increased. However, the $x$ and $y$ coordinate errors of the ensemble model
decreased when using these models.

$\epsilon$ Value | Type of Data |  Validation MSE ($\times 10^{-5}$) |  MSE ($\times 10^{-5}$)
---|---|---|---
|  |  | x | y | z
0.01 | Microphone | 1.44 | 0.49 | 0.29 | 3.55
Accelerometer | 1.33 | 0.45 | 0.28 | 3.26
0.001 | Microphone | 1.53 | 0.59 | 0.29 | 3.71
Accelerometer | 1.36 | 0.43 | 0.26 | 3.40
TABLE II: The validation MSEs and the MSEs of the x, y and z coordinates for
the ensemble model using all data corresponding to $\epsilon=0.01$ or
$\epsilon=0.001$ with microphone or accelerometer data.

To investigate this further, Figure 6 shows the average distance between the
ground truth and prediction for the two models using all data corresponding to
$\epsilon=0.01$ or $\epsilon=0.001$, with either microphone or accelerometer
data. These results show that regardless of the type of data, the ensemble
model’s prediction is very similar. This was also seen in Figure 5. This is to
be expected as the simulated microphone and accelerometer data are
fundamentally the same.

Figure 6: The average distance between the prediction using the ensemble model
and the ground truth, for each coordinate. Using all data corresponding to
either $\epsilon=0.01$ or $\epsilon=0.001$, with either microphone or
accelerometer data.

Further comparing this to the previous test using only one file, the overall
average distance between the truth and predicted $x$-coordinate decreased,
$y$-coordinate decreased and $z$-coordinate increased. This shows that,
alongside the MSEs, the $x$ and $y$ coordinate prediction of the ensemble
model improved when using all data corresponding to either $\epsilon=0.01$ or
$\epsilon=0.001$. This is because, the more data the ML algorithms have, the
more samples can be used for training. Therefore, the models can make more
accurate predictions.

Thirdly, two models were created using all data either corresponding to
$\epsilon=0.01$ or $\epsilon=0.001$. This contained data using polynomial
degree $r=1$ and $r=2$, as well as 10 $\times$ 5 $\times$ 4 mesh and 20
$\times$ 10 $\times$ 8 mesh. However, unlike the previous test, both
microphone and accelerometer data were combined. This resulted in all models,
including the ensemble model, having similar results. They were all unable to
accurately predict the source disturbance location. Therefore, the ML
algorithms underfitted as they were unable to learn the underlying pattern of
the data.

Lastly, two universal models were created using all data corresponding to both
$\epsilon=0.01$ and $\epsilon=0.001$. This contained data using polynomial
degree $r=1$ and $r=2$, as well as 10 $\times$ 5 $\times$ 4 mesh and 20
$\times$ 10 $\times$ 8 mesh. Furthermore, either microphone or accelerometer
data was used. These two types of data were investigated separately due to the
results from the previous test.

The parameters used for the RF model were, ‘n estimators’: 800 and ‘max
depth’: 25. The parameters used for the kNN model were, ‘n estimators’: 3 and
‘weights’: ‘distance’. Table III shows that MSEs of the ensemble model using
this configuration. Comparing these results to the second test, where the two
$\epsilon$ values were separated, the MSEs for all the coordinates have
decreased.

Type of Data |  Validation MSE ($\times 10^{-5}$) |  MSE ($\times 10^{-5}$)
---|---|---
|  | x | y | z
Microphone | 1.00 | 0.31 | 0.18 | 2.53
Accelerometer | 0.92 | 0.27 | 0.17 | 2.31
TABLE III: The validation MSEs and the MSEs of the x, y and z coordinates for
the ensemble model using all data combining $\epsilon=0.01$ and
$\epsilon=0.001$, with microphone or accelerometer data.

In addition, Figure 7 shows the average distance between the ground truth and
prediction for the universal models for microphone or accelerometer data.
Comparing these results to the previous tests, the overall average distance
between the ground truth and predicted coordinates has decreased, particularly
for the $x$ and $y$ coordinates. This shows that, alongside the MSEs, the
source disturbance prediction of the ensemble model improved when using all
data corresponding to a combination of $\epsilon=0.01$ and $\epsilon=0.001$.
This further consolidates that the value of $\epsilon$ does not impede the ML
models’ prediction.

Figure 7: The average distance between the prediction using the ensemble model
and the ground truth, for each coordinate. Using data corresponding to
$\epsilon=0.01$ and $\epsilon=0.001$ with microphone or accelerometer data and
an average of the microphone and accelerometer predictions.

An extension to this test was done by combining the results from each
universal model, either using microphone or accelerometer data, together. This
was done by averaging the two predictions to create a final prediction. Then,
the average distance between the ground truth and this final prediction was
calculated and the results can be seen in Figure 7. It shows that the $x$ and
$y$ coordinate prediction improved however, the $z$-coordinate prediction did
not improve compared to when using just the accelerometer data. This can also
be visualised with a 3D plot. Figure 8 shows the ensemble model’s results
using the average of the microphone and accelerometer prediction. Comparing
this to Figure 4, it shows an even smaller deviation from the ground truth and
prediction.

Figure 8: A 3D plot comparing the ground truth to the prediction on the
validation dataset for the ensemble model for a single example. Using a
combination of data corresponding to $\epsilon=0.01$ and $\epsilon=0.001$ with
the average of the microphone and accelerometer prediction. Figure is not to
scale.

The $z$-coordinate prediction may have not improved as the sensors are less
sensitive to signals along the $z$ axis. This may be due to the locations of
the microphones and accelerometers. This can be seen in Figure 2 of the
experimental setup. It shows that the sensors are more spread out along the
$x$ and $y$ axis compared to the $z$ axis.

A further test was conducted in order to investigate the effect of the number
of sensors used. This was done by using the previous universal model with
microphone data. However, the datasets utilised to train and validate the
model correspond to either using one, two or three sensors. The locations of
these sensors can be seen in Figure 2 and they correspond to the results seen
in Figure 9. These results shows that the more sensors used in the
experimental setup, the more accurate the models’ prediction is. However, this
is only up until a certain point as the results from using three and four
sensors are very similar. Therefore, the positions of the sensors may also
affect their sensitivity to signals and thus, the model’s ability to predict
the source position.

Figure 9: The average distance between the prediction using the ensemble model
and the ground truth, for each coordinate. Using data corresponding to
$\epsilon=0.01$ and $\epsilon=0.001$ with microphone data and a different
number of sensors.

##  V Discussion/Conclusion

In conclusion, simulated sensor readings were used to train ML models to
predict the source location of this signal, and thus detect the position of
the occlusion. Several regression ML algorithms were explored, as well as
experimenting with different configurations of the data. In the end, it was
found that ML models can predict, with high accuracy, the location of the
source as a result of the sensor readings.

Overall, there are several outcomes to this work. Firstly, the best performing
ML algorithm was the ensemble model. This ensemble model was a combination of
the kNN and Random Forest algorithms.

Secondly, the best performing configuration was using all data corresponding
to a combination of $\epsilon=0.01$ and $\epsilon=0.001$ with either
microphone or accelerometer data. Furthermore, the average of these two
predictions, to create a final estimate, resulted in the best performing
model. The average MSE of the $x$, $y$ and $z$ coordinates were 2.9$\micro$m2,
1.8$\micro$m2 and 0.024mm2 respectively as seen in Table III. The average
distance between the ground truth and prediction for the $x$, $y$ and $z$
coordinates were 1.1mm, 0.93mm, 3.2mm respectively as seen in Figure 7.

Thirdly, the results show that the narrowness of the source ($\epsilon$) does
not greatly affect the models’ prediction. This may be useful when building a
non-invasive detection method for future real-life scenarios.

Lastly, the results showed that nearly all models predicted the $y$-coordinate
of the source disturbance the most accurately. While the $z$-coordinate was
the most inaccurate. A possible method to improve these results would be to
spread the sensors along the $z$-axis more. However, this is limited to some
extent as the sensors are confined to the skin surface, so the z-distance of
the sensors cannot vary too much. Also, a further investigation is needed in
order to see if the sites at which the microphones and accelerometers are
located improves the prediction of the $z$-coordinate. This is critical in
determining the optimal sensor positions in the real device.

##  VI Future work

Future work for this investigation would be to add noise to the simulated
data. This is because, real-life readings are not perfect and often have
background noise. Also, using time-dependent data would also be beneficial.
These additions would simulate these situations more accurately and help
improve the ML models’ performance.

However, testing the ensemble model developed on real-life data would provide
the most meaningful results. This is because, the model’s performance can be
assessed and improved using the same type of data it will use during
deployment.

All data used in this investigation was numerically generated. However, the
results from this work are intended, in the future, to be used on real-life
CAD patient data. Therefore, legal, social, ethical and sustainability issues
such as privacy need to be considered when applying these ML models to
personal data.

## References

  * [1] Clipart Library. (2019) “Human Body Silhouette.” Available at: http://clipart-library.com/clipart/yikrLoqjT.htm (Accessed: 31/10/22).
  * [2] Clipart Library. (2019) “Unlabelled Diagram Of The Heart.” Available at: http://clipart-library.com/clipart/5cRrEGKbi.htm (Accessed: 31/10/22).
  * [3] Hesser, D. F., Kocur, G. K., & Markert, B. (2020). “Active source localization in wave guides based on machine learning.” Ultrasonics, 106, 106144. https://doi.org/10.1016/j.ultras.2020.106144
  * [4] Nasarian, Elham et al. (2020) “Association between work-related features and coronary artery disease: A heterogeneous hybrid feature selection integrated with balancing approach.” Pattern Recognit. Lett, 133, 33-40. https://doi.org/10.1016/j.patrec.2020.02.010
  * [5] Akella, A., & Akella, S. (2021). “Machine learning algorithms for predicting coronary artery disease: efforts toward an open source solution.” Future science OA, 7(6), FSO698. https://doi.org/10.2144/fsoa-2020-0206
  * [6] Doppala, B. P., Bhattacharyya, D., Janarthanan, M., & Baik, N. (2022). “A Reliable Machine Intelligence Model for Accurate Identification of Cardiovascular Diseases Using Ensemble Techniques.” Journal of healthcare engineering. https://doi.org/10.1155/2022/2585235
  * [7] Shaw, S., Greenwald, S. (2022). “Machine Learning Techniques for Source Localisation in Elastic Media, MSc Project Proposal.”
  * [8] Géron, A. (2017), “Hands-on Machine Learning with Scikit-Learn, Keras, and TensorFlow concepts, tools, and techniques to build intelligent systems.” O’Reilly Media. 978-1-492-03264-9.
  * [9] Pedregosa, F et al. (2011). “Scikit-learn: Machine Learning in Python.” Journal of Machine Learning Research, 12, 2825-2830. https://doi.org/10.48550/arXiv.1201.0490
  * [10] Song, Y. Y., & Lu, Y. (2015). “Decision tree methods: applications for classification and prediction.” Shanghai archives of psychiatry, 27(2), 130–135. https://doi.org/10.11919/j.issn.1002-0829.215044
  * [11] Chen. T, Guestrin. C. (2016). “XGBoost: A Scalable Tree Boosting System.” In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD ’16). Association for Computing Machinery, 785–794. https://doi.org/10.1145/2939672.2939785