# Modeling Interlayer Interactions and Phonon Thermal Transport in Silicene
Bilayer

Sapta Sindhu Paul Chowdhury  Appalakondaiah Samudrala  Santosh Mogurampelly
<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology
Jodhpur
Karwar, Jodhpur, Rajasthan, India - 342030.

###### Abstract

We develop an accurate interlayer pairwise potential derived from the ab-
initio calculations and investigate the thermal transport of silicene bilayers
within the framework of equilibrium molecular dynamics simulations. The
electronic properties are found to be sensitive to the temperature with the
opening of the band gap in the $\Gamma$$\rightarrow$M direction. The
calculated phonon thermal conductivity of bilayer silicene is surprisingly
higher than that of monolayer silicene, contrary to the trends reported for
other classes of 2D materials like graphene and hBN bilayers. This
counterintuitive behavior of the bilayer silicene is attributed to the
interlayer interaction effects and inherent buckling, which lead to a higher
group velocity in the LA1/LA2 phonon modes. The thermal conductivity of both
the mono- and bilayer silicene decreases with temperature as $\kappa\sim
T^{-0.9}$ because of the strong correlations between the characteristic
timescales of heat current autocorrelation function and temperature ($\tau\sim
T^{-0.75}$). The mechanisms underlying phonon thermal transport in silicene
bilayers are further established by analyzing the temperature induced changes
in acoustic group velocity.

††preprint: APS/123-QED

##  I Introduction

Silicene is an intriguing two dimensional (2D) material with unique in-plane
atomic configuration analogous to graphene and out-of-plane structural
similarities with other classes of 2D materials [1, 2, 3]. The structural
uniqueness inherent to silicene offers several distinguishable physical and
electronic properties [4]. Silicene has garnered tremendous interest in the
last decade due to the existence of remarkable flexural phonon scattering [5],
electro-optic effects [6], quantum anomalous Hall effect [7], spin-orbit
couplings [8], giant magnetoresistance [9], and tunable band structure [10].
The flexural phonon scattering in silicene limits the thermal conductivity to
low values, opening new avenues for technological applications such as
thermoelectric devices. More interestingly, the structural and compositional
superiority makes silicene compatible with silicene based electronics [11] and
state-of-the-art semiconductor fabrication technologies [12].

The bilayer silicene poses interesting challenges for technological
applications [13, 14, 15, 16], and therefore, understanding the underlying
physical mechanisms is essential. A majority of the investigations carried out
so far focused on the electronic properties of bilayer silicene [17, 18, 19,
20]. Despite significant progress in elucidating the electronic properties of
mono- and bilayer silicene 2D materials, the phonon thermal transport and the
underlying fundamental heat transport mechanisms are poorly understood [21].
Specifically, while there exist only a few computational works reporting the
thermal conductivity of monolayer silicene, there is little work [22] on the
thermal conductivity of the bilayer.

Understanding the thermal transport mechanisms in unique 2D materials like
silicene has fundamental importance from both the scientific and the
technological points-of-view [23]. Technological application of these 2D
materials in thermal interfacing applications, thermoelectric devices require
thorough understanding of phonon transport at higher temperatures [24]. As a
result, both theoretical and experimental efforts have been undertaken to
understand the thermal transport in these materials. Thermal conductivity of
monolayer, bilayer and twisted bilayers of graphene, h-BN etc. have been
studied extensively with monolayer graphene showing the highest thermal
conductivity. Similarly, many works have focused on the dependency of thermal
conductivity with the external factors like strain, temperature, twist etc.
The thermal conductivity in silicene and the interplay of interactions is not
well understood due to the challenges involved in modeling these materials.
[25, 24].

The computational studies of monolayer silicene focused on the thermal
conductivity and its dependency on temperature [26], strain [27, 28], vacancy
[29], and substrate [30]. The first principles based calculations estimated
the thermal conductivity of silicene to be 9.4 W/(m K) at room temperature,
mainly due to the in-plane vibrations [31]. Using the reactive force field
(ReaxFF) model, molecular dynamics (MD) study found the silicene monolayer to
be thermally stable up to a temperature of 1500 K [32]. Earlier work based on
the Tersoff bond order potential [33] based MD simulations estimated the
thermal conductivity to be 20 W/(m K) [34] but failed to reproduce the buckled
structure of silicene. Similarly, the Modified Embedded Atom Model (MEAM)
potential incorrectly predicts the buckling height (twice the reported value)
in silicene [35]. Stillinger-Weber (SW) potential parameters modified by Zhang
et al. [36] could reproduce the buckling of silicene monolayer in agreement
with experiments and first principle density functional theory (DFT) studies.
The modified SW potential parameter set produces the thermal conductivity of
the monolayer to be in the range 5-15 W/(m K) at 300 K, employing both
equilibrium and non-equilibrium molecular dynamics (EMD and NEMD) approaches
[36].

Despite the above, there is no consensus on how the thermal conductivity of
monolayer silicene scales with temperature and the corresponding phonon
mechanisms. More importantly, the investigations on phonon thermal properties
of multilayered silicene suffer a serious lack of the availability of a
suitable interlayer potential. In this work, we optimize the Lennard-Jones
(LJ) pairwise interaction potential to accurately capture the interlayer
interactions using the ab-initio DFT calculations. The optimized interlayer LJ
potential model is used for studying the structure and thermal transport
properties of the silicene bilayer. Specifically, we study the effect of
temperature on the structural, electronic, and thermal conducting properties
of the silicene bilayer. Moreover, we compare the thermal transport of the
silicene bilayer to that of the monolayer. The thermal conductivity results of
silicene are explained through the calculation of mode dependent phonon
properties of the bilayer and monolayer systems. The phonon characteristic
timescales associated with the heat current autocorrelation functions and the
phonon dynamics are thoroughly analyzed to examine the origins behind the
power-law scaling relation between thermal conductivity and temperature.

The rest of the paper is organized as follows: Section II presents the
theoretical formalism and computational details followed in the paper. Section
III presents the results obtained in the study, followed by a discussion, and
Section IV summarizes the key results of the paper. We present the structural
information in Section III.1, and the optimization of LJ parameters is
detailed in Section III.2. The effect of temperature on the structural
configurations is described in Section III.3 followed by the calculations of
thermal conductivity in Section III.4. Section III.5 explains the thermal
conductivity by analyzing phonon modes.

##  II Theoretical Background and Computational Details

We use the _ab-initio_ density functional theory framework as implemented in
the quantum espresso [37, 38] to predict the ground-state interlayer
separation between silicene layers. The exchange–correlation (XC) energy was
estimated using the Perdew–Burke–Ernzerhof (PBE) functional [39] of the
generalized gradient approximation (GGA). It is noteworthy mentioning that
GGA-PBE functional commonly leads to an underestimation of the van der Waals
(vdW) interactions in layered materials like graphene and h-BN, where the
$\pi$-bonding network plays a critical in determining their properties.
However, the interactions in bilayer or multi-layered silicene materials are
predominantly governed by strong covalent bonds that disrupt the $\pi$-bonding
network between the layers, and therefore, integrating the empirical vdW
correction schemes into our DFT calculations is unnecessary (more details are
available in the Supplemental Material [40]. See also references [41, 42, 43]
therein). The projector augmented wave (PAW) [44] method for pseudopotential
is used to address the electron core interactions. For our calculations, plane
waves with a kinetic energy cutoff of 110 Ry for wavefunction and 880 Ry for
charge density have been used with a Monkhorst-Pack [45] k-point grid of 32
$\times$ 32 $\times$ 1\. We have minimized the Hellmann-Feynman forces acting
on the atoms to $5\times 10^{-4}$ eV/Å, and the total energy convergence
threshold is $10^{-6}$ eV.

The molecular dynamics simulations have been carried out in the Large-scale
Atomic/Molecular Massively Parallel Simulator (Lammps) [46]. In our approach,
we have employed SW [47] potential to describe the intralayer interactions,
while the interlayer interactions are modeled using the LJ potential.
Specifically, the atoms within a given layer interact exclusively through the
SW potential, and the atoms of different layers interact solely via the LJ
potential. The parameters governing the SW potential have been adopted from
Zhang et al. [36], while the parameters for the LJ model are subject to
refinement within the scope of this study. The Hamiltonian used in our
classical physics calculations is given by:

$H=\text{K.E.}+\Phi_{\text{intra}}+\Phi_{\text{inter}},$ (1)

where the terms on the right-hand side represent the kinetic energy,
intralayer interaction potential, and interlayer interaction potential of the
system, respectively. The term, $\Phi_{\text{intra }}$ takes the following
form:

$\displaystyle\Phi_{\text{intra }}=$
$\displaystyle\sum_{i}\sum_{j>i}\sum_{l}\phi_{2}\left(r_{ijl}\right)$ (2)
$\displaystyle+\sum_{i}\sum_{j\neq
i}\sum_{k>j}\sum_{l}\phi_{3}\left(r_{ijl},r_{ikl},\theta_{ijkl}\right).$

In the above equation, $r_{ijl}$ represents the distance between
$i^{\text{th}}$ and $j^{\text{th}}$ atom of $l^{\text{th}}$ layer and
$\theta_{ijkl}$ represents the angle between $i^{\text{th}}$, $j^{\text{th}}$,
and $k^{\text{th}}$ atom of $l^{\text{th}}$ layer. The terms, $\phi_{2}$ and
$\phi_{3}$ are respectively given by:

$\displaystyle\phi_{2}\left(r_{ijl}\right)=$ $\displaystyle
A_{ijl}\epsilon_{ijl}\left[B_{ijl}\left(\frac{\sigma_{ijl}}{r_{ijl}}\right)^{p_{ijl}}-\left(\frac{\sigma_{ijl}}{r_{ijl}}\right)^{q_{ijl}}\right]$
(3)
$\displaystyle\times\exp\left(\frac{\sigma_{ijl}}{r_{ijl}-a_{ijl}\sigma_{ijl}}\right)$

and

$\displaystyle\phi_{3}\left(r_{ijl},r_{ikl},\theta_{ijkl}\right)=$
$\displaystyle\lambda_{ijkl}\epsilon_{ijkl}\left[\cos\theta_{ijkl}-\cos\theta_{0ijkl}\right]^{2}$
(4)
$\displaystyle\times\exp\left(\frac{\gamma_{ijl}\sigma_{ijl}}{r_{ijl}-a_{ijl}\sigma_{ijl}}\right)$
$\displaystyle\times\exp\left(\frac{\gamma_{ikl}\sigma_{ikl}}{r_{ikl}-a_{ikl}\sigma_{ikl}}\right).$

Here, $a,\lambda,\gamma,\cos\theta_{0},A,B,p,q$ are the force field parameters
defining the interaction potential terms, taken from Zhang et al. [36]. The
interlayer potential, $\Phi_{\text{inter}}$ is of the Lennard-Jones form:

$\Phi_{\text{inter}}=\sum_{i\in L_{1}}\sum_{j\in
L_{2}}\phi\left(r_{ij}\right),$ (5)

where,

$\phi\left(r_{ij}\right)=4\epsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right],r_{ij}<r_{c}.$
(6)

Here, $\sigma$ and $\epsilon$ are the parameters of the potential, $r_{ij}$ is
the distance between $i^{\text{th}}$ and $j^{\text{th}}$ atoms of different
layers. The global cutoff for LJ is chosen to be $r_{c}=20$ Å, which is more
than five times $\sigma$. For optimization of the LJ parameters for bilayer
silicene, we have calculated the binding energy using the equations:
$E^{\mathrm{DFT/MD}}_{\mathrm{b}}$ =
$E^{\mathrm{{DFT/MD}}}_{\mathrm{{bilayer}}}$ \-
$E^{\mathrm{{DFT/MD}}}_{\mathrm{{layer1}}}-E^{\mathrm{{DFT/MD}}}_{\mathrm{{layer2}}}$
where, $E^{\mathrm{{DFT}}}_{\mathrm{{b}}}$ and
$E^{\mathrm{{MD}}}_{\mathrm{{b}}}$ are the binding energy from DFT and MD
calculations respectively. Here, $E^{\mathrm{{DFT}}}_{\mathrm{bilayer}}$,
$E^{\mathrm{{DFT}}}_{\mathrm{layer1}}$ and
$E^{\mathrm{{DFT}}}_{\mathrm{layer2}}$ are the total energies of the bilayer,
and the monolayers from DFT calculations and
$E^{\mathrm{{MD}}}_{\mathrm{bilayer}}$, $E^{\mathrm{{MD}}}_{\mathrm{layer1}}$
and $E^{\mathrm{{MD}}}_{\mathrm{layer2}}$ are the total energies of the
bilayer and the monolayer from MD simulations. We minimize the following
objective function to get the optimized $\sigma$ and $\epsilon$:

$\chi^{2}=\frac{1}{W}\int_{0}^{\infty}|E^{\mathrm{DFT}}_{\mathrm{b}}(r)-E^{\mathrm{MD}}_{\mathrm{b}}(r)|^{2}w(r)dr$
(7)
$\approx\frac{1}{W}\sum_{r_{i}}^{r_{cut}}|E^{\mathrm{DFT}}_{\mathrm{b}}(r_{i})-E^{\mathrm{MD}}_{\mathrm{b}}(r_{i})|^{2}w(r_{i}),$
(8)

where, $r$ specifies the interlayer separation between two silicene layers and
$w(r_{i})$ is the weight function, defined as:

$w(r)=\exp\left[-\left(\frac{r-r_{0}}{\zeta}\right)^{2}\right],$ (9)

where, $r_{0}$ is the minima of LJ potential and $\zeta$ controls the
amplitude of $w(r)$ around the minima. We note that the choice of $w(r)$ is
not unique and a relevant weight function can be chosen so that the
significant part of the LJ potential is represented correctly.

We have applied the periodic boundary condition (PBC) in all three directions.
To create aperiodicity and to prevent the interactions between the bilayer
system and its images in the z-direction, we have created 50 Å of vacuum in
the simulation cell. Energy minimization using the conjugate gradient (CG)
[48] and steepest descent (SG) algorithms have been conducted with energy
tolerance of $10^{-8}$ and force tolerance of $10^{-8}$ eV/Å. The minimized
structures have been equilibrated in the isothermal-isobaric (NPT) ensemble at
different temperatures starting from 300 K to 1000 K and at 0 bar pressure for
300 ps. The Nosé-Hoover thermostat [49, 50] and the Nosé-Hoover barostat [51]
with a temperature coupling constant of 0.1 ps and pressure coupling constant
of 1 ps have been used for maintaining a constant temperature and pressure,
respectively. The systems are then equilibrated in a canonical (NVT) ensemble
for 300 ps. The Verlet algorithm [52] with timestep 1 fs has been used to
integrate the equation of motion.

For the generation of production trajectories for the calculations of thermal
conductivity, we simulated larger trajectories in an NVT ensemble. A total of
30 independent trajectories are generated using uncorrelated initial
conditions. Each independent trajectory of length 2 ns is generated using a
timestep of 0.1 fs.

We have used the Green-Kubo (GK) relation based on the fluctuation-dissipation
theorem to calculate the thermal conductivity of the systems using the
following relation:

$\bm{\kappa}(T)=\frac{1}{Vk_{B}T^{2}}\int_{0}^{\infty}\langle\bm{S}(0).\bm{S}(t)\rangle
dt,$ (10)

where, $V$ is the volume, $k_{B}$ is the Boltzmann constant, and $T$ is the
temperature of the system. The ensemble average
$\langle\bm{S}(0).\bm{S}(t)\rangle$ represents the autocorrelation function of
the heat current operator $\bm{S(t)}$, which is calculated from the simulation
data as:

$\bm{S}(t)=\frac{d}{dt}\sum_{i}\bm{r}_{i}\tilde{E}_{i},$ (11)

where $\bm{r}_{i}$ is the position vector concerning the ith atom,
$\tilde{E}_{i}=E_{i}-\langle E_{i}\rangle$ represents the fluctuation in total
energy with respect to the mean energy. $\bm{S}(t)$ is numerically computed in
the following way:

$\bm{S}(t)=\sum_{i}\tilde{E}_{i}\bm{v}_{i}+\frac{1}{2}\sum_{i<j}(\bm{F}_{ij}.(\bm{v}_{i}+\bm{v}_{j}))\bm{r}_{ij},$
(12)

where $\bm{F}_{ij}$ is the force between atom $i$ and atom $j$, $\bm{v}_{i}$
is the velocity and $\bm{r}_{ij}$ is the displacement of ith and jth atoms.

The phonon density of states (DoS) of the systems are calculated at different
temperatures using the velocity autocorrelation function (VACF) as follows:

$D(\omega)=\frac{1}{3Nk_{B}T}\int_{0}^{\infty}\frac{\langle\bm{v}(0).\bm{v}(t)\rangle}{\langle\bm{v}(0).\bm{v}(0)\rangle}e^{i\omega
t}dt,$ (13)

where $\langle\bm{v}(0).\bm{v}(t)\rangle$ defines the VACF, $\omega$ is the
angular frequency, $N$ is the number of atoms in the system. Notably, a
simulation trajectory of 250 ps with a finer timestep of 0.05 fs and a saving
frequency of 0.2 fs is used for the NVT run to capture higher resolution in
VACF and obtain a smoother DoS. The trajectory is recorded at a time interval
of 0.2 fs.

We use the equilibrium statistical mechanics formalism based on the
fluctuation-dissipation theorem to construct the dynamical matrices [53] of
phonon wave vector using the following equation:

$\bm{D}_{k\alpha,k^{\prime}\beta}(\bm{q})=(m_{k}m_{k^{\prime}})^{-\frac{1}{2}}\bm{\Phi}_{k\alpha,k^{\prime}\beta}(\bm{q}).$
(14)

Here, the force constant coefficient
$\bm{\Phi}_{k\alpha,k^{\prime}\beta}(\bm{q})$ at phonon wavevector $\bm{q}$ is
obtained from:

$\bm{\Phi}_{k\alpha,k^{\prime}\beta}(\bm{q})=k_{B}T\bm{G}^{-1}_{k\alpha,k^{\prime}\beta}(\bm{q}).$
(15)

If the displacement of the kth basis atom in the unit cell in the $\alpha$
direction is:

$\bm{u}_{k\alpha}(\bm{q})=\sum_{l}\bm{u}_{lka}\exp(i\bm{q}.\bm{r}_{l}),$ (16)

then the Green’s function coefficient can be given as:

$\bm{G}_{k\alpha,k^{\prime}\beta}(\bm{q})=\langle\bm{u}_{k\alpha}(\bm{q})\cdot\bm{u}^{\ast}_{k^{\prime}\alpha}(\bm{q})\rangle.$
(17)

To calculate the phonon dispersion, we have used a post-processing tool, phana
[54] along with the fix-phonon package [55] as implemented in lammps. The
group velocity is calculated as follows:

$\begin{split}\bm{v}_{g}(\bm{q}\nu)&=\nabla_{q}\omega(\bm{q}\nu)\\\
&=\frac{\partial\omega(\bm{q}\nu)}{\partial\bm{q}}\\\
&=\frac{1}{2\omega(\bm{q}\nu)}\frac{\partial[\omega(\bm{q}\nu)]^{2}}{\partial\bm{q}}\\\
&=\frac{1}{2\omega(\bm{q}\nu)}\left<\bm{e}(\bm{q}\nu)|\frac{\partial
D(\bm{q})}{\partial\bm{q}}|\bm{e}(\bm{q}\nu)\right>.\end{split}$ (18)

Here, $\bm{e}(\bm{q}\nu)$ is the phonon eigenvector at a frequency $\nu$ at
$\bm{q}$. We have used phonopy [56, 57] interfaced with phono-lammps for
calculating the group velocity at different temperatures.

Figure 1: (a) Crystal structure [36] and the snapshot of equilibrated
monolayer silicene displaying out-of-plane corrugations in vacuum at 300 K (b)
crystal structure and the snapshot of equilibrated bilayer silicene displaying
out-of-plane corrugations at 300 K. Due to the substrate-like effect on the
silicene as a result of interlayer interactions, we observe suppressed spatial
fluctuation on the bilayer system compared to the monolayer

##  III Results and Discussion

###  III.1 Structural Configuration

Figure 2: Binding energy with interlayer separation for bilayer silicene. The
black curve is the reference DFT data used for deriving LJ parameters. The
optimized parameters are used for all subsequent MD simulations. The red-
colored curve is the binding energy curve obtained using the optimized LJ
parameters ($\epsilon$ = 0.0346 eV and $\sigma$= 2.586 Å).

We consider the low-buckled geometry of monolayer silicene, which is the most
stable structure of silicene [21]. The hexagonal lattice with a lattice
constant of 3.81 Å, Si-Si bond length of 2.24 Å with a buckling height of
0.427 Å is taken for our calculation as per previous reports [36] (Fig. 1(a)).
The unit cell follows symmetry operations as followed by the C2/m spacegroup.
We extend this unit cell to create a 32 $\times$ 32 $\times$ 1 supercell with
2048 Si atoms (Fig. 1(b)) for subsequent MD simulations. We use the SW2
potential parameters provided by Zhang et al. [36] to describe interlayer
interactions. The parameter set is based on the acoustic part of the phonon
dispersion curve responsible for the maximum contribution to thermal
conductivity. The optimized geometry obtained using the SD and CG minimization
protocols agrees well with the lattice parameters reported. Our DFT
calculations confirm the electronic band structure with the Dirac cone at 0 K
point and the non-magnetic nature of the monolayer in the ground state [21].
To construct the bilayer, we align one layer on top of the other monolayer so
that an AA stacking is obtained. Since no known suitable and accurate
interlayer interaction potential exists for bilayer silicene, we develop the
LJ potential model. The additional degree of freedom unique to silicene
results in a multitude of stacking possibilities, encompassing various
buckling and lattice constants. These include AA-planer, AA-buckled, AA′, AB,
AB′, A′B′, which were explored theoretically to comprehend electronic
properties [42, 43, 13]. Our forthcoming work will focus on developing
manybody interlayer interaction potentials tailored for silicene, enabling an
understanding of thermal transport in multilayer buckled silicene structures.
In the present study, we employ the pairwise potential for the AA-stacked
bilayer silicene to study thermal transport, acknowledging the complexity of
more comprehensive manybody potentials for diverse stackings.

###  III.2 Optimization of the Force Field Parameters

Figure 3: The variation of structural parameters with temperature: (a)
interlayer separation, (b) buckling height, and (c) percentage deviation of
the lattice constant with respect to the corresponding value at 0 K. (d) The
electronic band structure and (e) the phonon dispersion along the
$\Gamma\rightarrow K\rightarrow M\rightarrow\Gamma$ symmetry directions. The
solid black line represents the Fermi level in the electronic band. The solid
curve represents 0 K, and the dashed curve represents 300 K temperature for
both band structures. Figure 4: (a) The averaged (over 30 independent runs)
normalized HCACF as a function of correlation or lag time at 300 K, (b) the
thermal conductivity calculated using 30 independent trajectories (dashed
black curve), and the average thermal conductivity (solid black curve) with
the upper limit of the Green-Kubo integral (Eq. 10) at 300 K, and (c) the
thermal conductivity of silicene mono- and bilayer as a function of
temperature. The filled-in circles and the solid line represents the
simulation data and the fit to power law $\kappa\sim T^{-0.93}$ for the
monolayer system, respectively. The open circles and the dashed line represent
the simulation data and the fit to power law $\kappa\sim T^{-0.88}$ for the
bilayer system. The power-law fitting shows the decreasing trend of the
thermal conductivity with temperature for both the monolayer and bilayer. The
bilayer silicene system shows a higher magnitude of thermal conductivity than
the monolayer.

Developing an appropriate potential for modeling the interlayer interactions
is a challenging yet crucial task. In the past, many interlayer potentials
such as DRIP, KC, and ILP have been developed for graphene and other similar
2D materials [58, 59, 60]. Though robust and effective in modeling complex
geometry such as twisted bilayers and multilayers, simulations with these
potentials are computationally highly expensive. Here, we aim to explore the
pristine bilayer system without twisting or rotational asymmetry. As such, for
simplicity, we develop the pairwise LJ potential to model the interactions
between two silicene layers. The reference force field parameters $\epsilon$
and $\sigma$ are based on the unified force field (UFF) model of Rappe et al.
[61] and a global cutoff $r_{c}=20$ Å which is more than five times $\sigma$
is taken. Previously, in many studies concerning silicene, Si-Si interactions
are modeled using the UFF parameters. However, the interlayer separation (d)
produced by these parameters is inaccurate. Using the UFF LJ parameters,
molecular mechanics minimization produces an interlayer separation of 3.95 Å,
which we plan to optimize by tuning $\epsilon$ and $\sigma$. However, our DFT
calculations with PBE-GGA density functional show the binding energy is
minimum at an interlayer separation of 2.79 Å compatible with previous studies
[43, 42]. Therefore, we optimize the LJ parameters by taking the UFF as the
starting approximation. We have plotted the binding energy curve as a function
of the interlayer separation in Fig. 2. It is clear that the UFF model
parameters are far off from our DFT calculation and the literature [43, 42].
We obtain the minimized $\epsilon$ and $\sigma$ by optimizing the $\chi^{2}$
function (Eq. (8)) following the simplex algorithm developed by Nelder and
Mead [62]. The developed force field parameters for the LJ potential are
$\epsilon$ = 0.0346 eV and $\sigma$ = 2.586 Å using weight function parameter
$\zeta$ = 0.3 Å. Here, we have introduced the weight function parameter to
model the minimum of the binding energy curve accurately. Based on this
optimized interlayer potential and the intralayer SW potential, the unit cell
of the bilayer considered in this work follows the symmetry operations of the
P-3m1 spacegroup. Using the developed LJ force field parameters, the obtained
lattice constant is 3.7941 Å, the buckling height is 0.460 Å, and the bond
length is 2.238 Å. In Table 1, we summarize the structural information of the
systems used in our study. Clearly, there is a slight increase in the buckling
height while the lattice constant decreases by a bit. We observe that the
optimized LJ parameters obtained by us produce the interlayer separation as
well as lattice parameters and buckling height of bilayer silicene consistent
with the first-principles DFT calculations and in agreement with the
literature [42, 43] at 0 K. Our DFT calculations with both ferromagnetic and
anti-ferromagnetic orderings reveal that the bilayer is most stable in the
non-magnetic state. Further details on the validity of our potential is given
on Supplementary Material [40](see also references [59, 63, 58, 17] therein) .
Despite satisfactorily reproducing the structural features of bilayer silicene
correctly, we note that the optimized LJ potential in this work has limited
capabilities in the sense that it is most accurate only for the AA-stacked
silicene, which it is based on.

Table 1: The crystal structures of the monolayer and bilayer silicene in the ground state (0 K) used in our calculations.  | a (Å) | Δh (Å) | Space Group
---|---|---|---
Monolayer | 3.812 | $0.427$ | C2/m
Bilayer | 3.794 | $0.460$ | P-3m1

###  III.3 Effect of Temperature on Structural Features

Because of the buckled topology and structural complexity, silicene is
expected to show high sensitivity to temperature. In this section, we discuss
the results of structure, electronic and phononic properties as a function of
temperature. We use a molecular dynamics trajectory of 500 ps in the
isothermal-isobaric ensemble to calculate the structural parameters for the
bilayer system at all temperatures. Prior to that, we equilibrate the system
in the canonical ensemble for 500 ps at the respective temperatures. The
results of temperature dependency of interlayer separation, buckling height,
and the change in lattice constant are explored in Fig. 3. We observe a
monotonic increase in the interlayer distance between the layers as the
temperature is increased (Fig. 3(a)). This results from the repulsion between
layers due to the increased thermal motion of the atoms at high temperatures.
However, a sharp decrease in buckling height is observed with temperature till
500 K due to the fact that the lattice constant is increasing (Fig. 3(b)). The
buckling height is apparently a constant in the range of 500 K - 750 K and
slightly increases till 1000 K. Compared to the monolayer, the buckling height
is less sensitive to temperatures beyond 500 K. We find that there is little
change in the lattice constant with the temperature in bilayer compared to the
monolayer. In Fig. 3(c), we plot the % deviation of the lattice constant with
respect to the value at 0 K. The maximum observed deviation in the lattice
constant is only 0.8% at an elevated temperature of 1000 K.
Counterintuitively, we observe that the lattice constant increases with
temperature till 300 K and shows a decreasing trend till 500 K. This is not
the case in the monolayer silicene, indicating interlayer interactions and
changes in buckling height might collectively result in the unusual trends for
lattice constant. In Fig. S2(a) (in the Supplemental Material [40]), we find
that as the temperature increases, there is a monotonic decrease in the
lattice constant from $\sim$3.81 to $\sim$3.71 Å as the temperature is
increased from 0 K to 1000 K. However, it is well known that the lattice
thermal expansion of 2D materials varies widely with the choice of interatomic
potentials [64, 65, 66, 67]. An accurate description of these properties will
require the use of ab-initio molecular dynamics simulations (AIMD). We will
address the behavior thermal expansion of bilayer silicene with accurate AIMD
in future communication. As expected, we observe significant corrugations in
the monolayer and bilayer silicene (Fig. 1) at finite temperatures. The degree
of corrugations is slightly suppressed in bilayer silicene compared to its
monolayer counterpart. In the case of the bilayer, one layer of silicene acts
as the substrate stabilizing the other layer, compatible with previous studies
on silicene based heterostructures [68].

The electronic structure and phonon properties exhibit pronounced sensitivity
to structural changes induced by temperature variations. To investigate the
effect of temperature, we employ DFT to calculate the electronic band
structure along the $\Gamma$$\rightarrow$K$\rightarrow$M$\rightarrow$$\Gamma$
direction (Fig. 3(d)) for the structures predicted by classical simulations at
0 and 300 K. Surprisingly, our calculations indicate the presence of a Dirac
cone at 0 K, contrary to previous findings [69]. This difference highlights a
significant issue with the structures predicted by classical molecular
mechanics. It serves to emphasize that the SW+LJ model, even after refinement
in our study, is inadequate for predicting structures intended to be
subsequently used when performing electronic structure calculations using DFT.

We also calculate the temperature effects on phonon dispersion in the same
high symmetry direction as the electronic band structure in the $\bm{q}$-space
(Fig. 3(e)). We observe that both the optical and acoustic branch of the
phonon spectrum undergoes red-shifting as the temperature is increased with a
higher redshift in the optical phonon branch compared to the acoustic branch.
More details about the phonon modes and the effect of temperature are
discussed in Section III.5. In summary, we see significant and unconventional
changes in the structural, electronic, and phononic behavior of the bilayer
compared to the monolayer counterpart with the increase in temperature. In the
following section, we shed light on how the above results affect thermal
transport in the bilayer system.

###  III.4 Thermal Conductivity of Bilayer Silicene

Figure 5: (a) The characteristic time scales associated with short-range
($\tau_{1}$) and long-range ($\tau_{2}$) correlations as a function of
temperature. The filled-in symbols and solid lines represent the simulation
data and the power-law fit for the monolayer system, respectively. The open
symbols and the dashed lines represent the simulation data, and the power law
fit for the bilayer system. Short-range characteristic time scales are shown
in diamonds, and the long-range characteristic time scales are shown in
circles. (b) Thermal conductivity as a function of the long-range
characteristics time scales. Other legends are the same as that of Fig. 4. It
is evident that the temperature has the same effect on $\tau$ as it has on
$\kappa$, which explains convincingly that $\kappa\sim T^{-0.93}$. Figure 6:
(a) Phonon DoS as a function of the angular frequency of vibration for
monolayer (solid black curve) and bilayer (solid red curve) silicene. The
absolute value is rescaled to better represent the dependency on the
interactions. Flexural contribution is plotted with dashed lines (black and
red for monolayer and bilayer, respectively). (b) Dependency of phonon DoS on
temperature for bilayer silicene. Black, red, and blue curves represent 300 K,
500 K, and 1000 K, respectively. The intensity of phonon density of state
decreases with temperature due to increased phonon scattering rate at high
temperatures.

The thermal transport properties of a material are characterized by thermal
conductivity. Using the EMD simulations, we calculate the thermal conductivity
of the silicene bilayer. The equilibrium approach uses the Greek-Kubo relation
(Eq. 10), which integrates the heat current autocorrelation function (HCACF)
to calculate the thermal conductivity. In Fig. 4(a), we calculate the
normalized HCACF of bilayer silicene as a function of correlation time or
lagtime, $\tau$. We observe that the correlation function decays exponentially
with oscillations till $\sim$ 50 ps, beyond which there is a negligible
contribution to thermal conductivity. For calculating the thermal
conductivity, we simulated 30 independent trajectories with a length of 1 ns
with a trajectory saving frequency of 0.02 fs. We note that the choice of the
number of independent trajectories strictly depends on the convergence of the
system and the desired accuracy. Fig. 4(b) shows the thermal conductivity
calculated from individual runs and the average thermal conductivity based on
30 independent runs as a function of the upper limit of the GK integral. The
thermal conductivity converges after about $\sim$20 ps of correlation time for
300 K. For higher temperatures, $\kappa$ converges faster, which can be
attributed to the lower lifetime of phonons at higher scattering. However, to
reduce errors, we have taken 100 ps as the cutoff correlation time, based on
the decay of the heat current autocorrelation function for estimating the
thermal conductivity.

Our calculations show that the thermal conductivity of the silicene bilayer is
14.8 $\pm$ 0.3 W/(m K) at room temperature, much lower than planer 2D
materials such as graphene and h-BN. We see an unusual increase in the thermal
conductivity of about $\sim 15\%$ in bilayer silicene compared to the
monolayer system. This is contrary to previously reported trends for graphene
and h-BN [70, 71]. A further understanding of this result requires a detailed
analysis of phonons conducted in III.5. We notice a monotonic decrease in the
thermal conductivity of both the mono- and bilayer systems with temperature
(Fig. 4(c)), similar to other 2D materials. We attribute this decrease in
thermal conductivity to phonon-phonon Umklapp processes. At higher
temperatures, the Umklapp process increases significantly, thus reducing the
thermal conductivity. Since the thermal conductivity is ideally found to be
inversely proportional to temperature, we fit $\kappa$ to a power-law and
obtain $\kappa\sim T^{-0.93}$ for monolayer and $\kappa\sim T^{-0.88}$ for
bilayer silicene, with an overall trend to be $\kappa\sim T^{-0.9}$. The
deviation from the ideal scaling of $T^{-1}$ for both systems is only
marginal. This deviation from ideal $T^{-1}$ law might be to the quasi-
harmonic nature of the scattering resulting from higher order effects [72].

To further understand the deviation from ideal $\kappa\sim T^{-1}$ law, we fit
the correlation functions with two exponential functions such that:

$\frac{\langle\bm{S}(0).\bm{S}(t)\rangle}{\langle\bm{S}(0).\bm{S}(0)\rangle}=a_{0}\exp(-t/\tau_{1})+(1-a_{0})\exp(-t/\tau_{2}).$
(19)

The fitting gives two characteristic times, a short-range characteristic time,
$\tau_{1}$, and a long-range characteristic time $\tau_{2}$). We find that
$\tau_{1}$ is less than $\tau_{2}$ signifying a faster decay at a lower
correlation time (below $\sim$5 ps). Further, both the $\tau_{1}$ and
$\tau_{2}$ decrease with the increase in temperature (Fig. 5(a)). This
indicates a lower phonon mean free path with the increase in temperature [72].
The lowering of the phonon mean free path contributes to the decrease in
thermal conductivity at higher temperatures. We analyze the correlation
between characteristic times and temperature in Fig. 5(a), and by fitting to a
power-law, we see that for both the monolayer and bilayer systems scale as
$\tau_{1}\sim T^{-0.75}$ and $\tau_{2}\sim T^{-0.77}$. To investigate the
existence of any scaling relation between the characteristics times and the
thermal conductivity, in Fig. 5(b), we compare $\kappa$ and $\tau_{2}$, the
long-range characteristic time as it contains more intricacies of phonons. It
is found that the thermal conductivity has a linear relation with the
characteristics time ($\kappa\sim\tau_{2}$) qualitatively. Although this
simple analysis provides a significant outlook, to further understand the
behavior of bilayer silicene compared to other 2D materials, we explore the
dependency of phonon modes in systems with temperature in the subsequent
section. We note that the quantum corrections are not incorporated in our
calculations due to the unreliable estimate of the quantum effects using these
techniques. No classical counterpart of zero point energy exists; hence every
quantum correction is deemed unreliable [73]. Since the main goal of this work
is to compare and understand the effect of temperature and interactions on
thermal transport, our classical results are valid for all considerations.

###  III.5 Phonons in Bilayer Silicene

Figure 7: (a) Phonon group velocity with angular frequency in bilayer silicene
at 300 K, (b) variation of the group velocity for the bilayer at 300 K, 500 K,
1000 K temperatures, (c) variation of the group velocity for monolayer (black
circles) and bilayer (red circles) and (d) comparison of phonon band structure
for monolayer (black curve) and bilayer (red curve) at 300 K. The acoustic
branches are represented with solid curves while the optical branches are
represented with dashed curves. We observe that group velocity decreases with
increased temperature, leading to reduced thermal conductivity. However,
LA1/LA2 mode velocities are increased in the bilayer compared to the monolayer
leading to higher thermal conductivity.

To understand the effect of interlayer interactions and temperature on thermal
conductivity, we have systematically analyzed the phonon dynamics in the
silicene systems. The phonon density of states provides insights into how the
phonon modes contribute to the thermal conductivity of the systems. We
calculate the phonon density of states and the phonon partial density of
states to understand how the phonon density is influenced by the in-plane and
out-of-plane motion of the atoms. The phonon density of states, calculated
from the MD simulations, is simply a Fourier transform of the VACFs. We plot
the normalized VACF and the contribution of the out-of-plane component to the
VACF as a function of the correlation or lagtime in Fig. S3(a) (in the
Supplemental Material [40]). We find that the VACF decays to zero after $\sim$
10 ps of correlation time. Hence, the maximum time limit for the integration
in Eq. (13) is taken as 50 ps, which is more than sufficient for all the
fluctuations to die out. Interestingly, the out-of-plane component shows
lesser fluctuations over the whole range of correlation time. We plot the
rescaled phonon DoS for the silicene bilayer and monolayer system in Fig.
6(a). We note that the acoustic mode density blueshifts in the bilayer, which
result in an increased thermal conductivity of the bilayer silicene compared
to the monolayer. The flexural component of the DoS contributes to total DoS
at lower frequencies but decays to zero at the optical frequencies.
Nevertheless, both the acoustic and optical modes are blueshifted for the
bilayer system. For understanding the effect of temperature on phonon DoS, the
total phonon density of states for the bilayer is plotted at 300 K, 500 K, and
1000 K temperatures in Figure 6(b). It is explicit that there is a monotonic
decrease in intensity in the phonon dispersion curve with temperature. The
redshift in the peak on the phonon DoS (Fig. S3(b), in the Supplemental
Material [40]) directly correlates to the decrease in thermal conductivity
with temperature for both the silicene systems [74]. This softening of the
phonon modes with temperatures is, in fact, due to the increase in the
interatomic and interlayer separations.

In Fig. 7(a), we have plotted the group velocity of the bilayer system at 300
K. From the calculations, it is found that the flexural ZA1/ZA2 phonon modes
have less group velocity than the other acoustic modes, i.e., TA1/TA2 and
LA1/LA2. This observation points to the fact that due to the inherent
buckling, the flexural modes are prone to higher scattering than the in-plane
modes and, thus, possess lower group velocity. Naturally, the thermal
conductivity due to flexural modes is expected to be limited, unlike in the
cases of graphene and h-BN [70]. However, in silicene, the LA modes contribute
the most to the thermal conductivity [30]. Among the acoustic modes, the
LA1/LA2 modes have higher group velocity in bilayer silicene (Fig. 7(a)). We
observe a blueshift of phonon branches apparent in the bilayer compared to the
monolayer system (Fig. 7(b)), similar to the blueshift seen in the phonon DoS.
Overall, it is found that the bilayer phonon modes have higher group
velocities compared to the monolayer system. Although bilayer silicene is
found to possess lower flexural group velocity (Fig. 7(b)), the much higher
magnitude of the group velocity of the LA1/LA2 modes explains the unusually
higher thermal conductivity of the bilayer silicene compared to the monolayer
counterpart. With temperature, there is a monotonic decrease in the group
velocity of the bilayer system, shown in Fig. 7(c). As the temperature is
raised, the scattering probability of phonons increases, leading to higher
collisions which decrease the group velocity. The reduced group velocity leads
to a lowered thermal conductivity of the systems, as discussed in the previous
section. Clearly, the redshifting of phonon modes directly correlates with the
decrease in thermal conductivity.

We use the classical molecular dynamics trajectories to calculate the phonon
dispersion of the monolayer and bilayer systems. This approach of calculating
the phonon dispersion agrees well with the phonon dispersion for the monolayer
calculated using the density functional perturbation theory (DFPT) (Fig. S4(a)
in the Supplemental Material [40])) and the phonon spectral energy density
(SED) method (Fig. S4(b) in the Supplemental Material [40])) in the
$\Gamma\rightarrow$ M direction. Since our method with the classical
potentials reproduces the results of quantum calculations within a reasonable
limit for uncertainties (Fig. S4(a) in the Supplemental Material [40])), we
use the classical mechanics based techniques for simplicity. The phonon
dispersion curve calculated for the bilayer (Fig. S4(c) in the Supplemental
Material [40])) and the monolayer systems (Fig. S4(a) in the Supplemental
Material [40])) shows a similar trend. We notice a blueshift in both the
optical and the acoustic branches of the bilayer compared to the monolayer,
which directly corresponds to the increase in the thermal conductivity in the
bilayer (Fig. 7(d)). It is found that the in-plane and out-of-plane modes are
strongly coupled in monolayer silicene due to inherent buckling [36]. The
flexural acoustic modes also experience a dampening effect due to the
increased buckling induced out-of-plane scattering. However, because strong
interlayer interactions exist between the silicene layers, the bonding leads
to the generation of high-velocity LA1/LA2 elastic waves. Generally, this is
observed in the case of substrate layer interactions, where increasing the
interaction strength leads to higher thermal transport [75, 76]. This high
velocity of the LA1/LA2 modes contributes the most to the thermal conductivity
of silicene, thereby increasing the thermal conductivity of the bilayer
compared to the monolayer system. To examine the effect of interaction
strength on the group velocity of silicene bilayers, we have carried out a few
simulation experiments by manually changing essentially the interaction
parameter $\epsilon$ appearing in the LJ interaction potential. To implement
this, we introduced a scaling parameter $\lambda$ in the LJ potential such
that,
$\phi\left(r_{ij}\right)=4\lambda\epsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right]$.
We have calculated the group velocities associated with the LA phonon modes of
the bilayer silicene systems for two different values of $\lambda$ (i.e., 0.5
and 1.5), to change the strength of interaction by 50 % weaker and 50 %
stronger with respect to the value of interaction parameter $\epsilon$
optimized in this work. In short, we observe that the group velocities of LA
phonons have a positive correlation with the strength of interaction in the
studied range of interaction strength (Fig. S5 in the Supplemental Material
[40])). Further, the thermal conductivity also exhibits similar changes with a
reduced thermal conductivity of 12.9 $\pm$ 0.4 W/(m K) for $\lambda$ = 0.5,
consistent with the changes observed for the group velocities. Thus, the
enhancement of LA1/LA2 phonon group velocities induced by stronger interlayer
interactions ultimately contributes to unusually higher thermal conductivity
in the bilayer silicene.

##  IV Conclusions

In summary, we developed accurate parameters of pairwise LJ potential for
modeling interlayer interactions between silicene layers and investigated the
phonon thermal transport of mono- and bilayer silecene with the Green-Kubo
approach. The developed set of parameters correctly predicts the interlayer
separation between the layers and phonon band structure in agreement with the
ab-initio DFT calculations [17].

The interlayer interactions and temperature play a major role in determining
the thermal transport properties of the silicene layered materials. Notably,
the behavior of ZA1/ZA2 modes in bilayer silicene exhibits a reduction in
group velocity when compared to the monolayer silicene systems. However,
intriguingly, the thermal conductivity in bilayer silicene is found to be
higher than that in monolayer silicene at any given temperature, which is
attributed to the enhancement in group velocity of LA1/LA2 modes arising due
to strong covalent interactions. Furthermore, the temperature-dependent
variation of silicene’s thermal conductivity follows a power-law relationship,
consistent with the trends observed in other two-dimensional materials. This
decrease in thermal conductivity follows a scaling trend of $\kappa\sim
T^{-0.9}$, which diverges from the typically reported $\kappa\sim T^{-1}$
scaling, owing to the quasi-harmonic nature of phonon scattering at elevated
temperatures.

Our investigation has revealed novel insights into the intricate phonon
transport mechanisms inherent in bilayer silicene, facilitated by the
development of interlayer pairwise interaction potential parameters.
Specifically, our analysys establishes that the phonon thermal conductivity
correlates excellently with the underlying relaxation timescales associated
with the decay of heat current autocorrelation functions. In addition to
understanding the phonon mechanisms in slicene systems, our study may pave
forward their potential utility in thermal interfacing devices and
thermoelectric applications. We propose further experimental as well as
theoretical studies at first-principles level to further establish the
complete nature of thermal transport in silicene.

###### Acknowledgements.

The authors acknowledge the Computer Center of IIT Jodhpur and the HPC Center
at the Department of Physics, Freie Universität Berlin
(10.17169/refubium-26754), for providing computing resources that have
contributed to the research results reported in this paper. SM acknowledges
support from the SERB International Research Experience Fellowship
SIR/2022/000786 and SERB CRG/2019/000106 provided by the Science and
Engineering Research Board, Department of Science and Technology, India. SSPC
acknowledges the Ministry of Education (MoE), Govt. of India, for the
financial support received as a fellowship. AS acknowledges support from the
SERB through grant number RJF/2021/000147.

## References

  * Vogt  _et al._ [2012] P. Vogt, P. De Padova, C. Quaresima, J. Avila, E. Frantzeskakis, M. C. Asensio, A. Resta, B. Ealet, and G. Le Lay, Silicene: Compelling experimental evidence for graphenelike two-dimensional silicon, Phys. Rev. Lett. 108, 155501 (2012).
  * Lin  _et al._ [2013] C.-L. Lin, R. Arafune, K. Kawahara, M. Kanno, N. Tsukahara, E. Minamitani, Y. Kim, M. Kawai, and N. Takagi, Substrate-induced symmetry breaking in silicene, Phys. Rev. Lett. 110, 076801 (2013).
  * Kharadi  _et al._ [2020] M. A. Kharadi, G. F. A. Malik, F. A. Khanday, K. A. Shah, S. Mittal, and B. K. Kaushik, Review—silicene: From material to device applications, ECS Journal of Solid State Science and Technology 9, 115031 (2020).
  * Fleurence  _et al._ [2012] A. Fleurence, R. Friedlein, T. Ozaki, H. Kawai, Y. Wang, and Y. Yamada-Takamura, Experimental evidence for epitaxial silicene on diboride thin films, Phys. Rev. Lett. 108, 245501 (2012).
  * Tao  _et al._ [2015] L. Tao, E. Cinquanta, D. Chiappe, C. Grazianetti, M. Fanciulli, M. Dubey, A. Molle, and D. Akinwande, Silicene field-effect transistors operating at room temperature, Nature Nanotechnology 10, 227 (2015).
  * Matthes  _et al._ [2013] L. Matthes, O. Pulci, and F. Bechstedt, Massive dirac quasiparticles in the optical absorbance of graphene, silicene, germanene, and tinene, Journal of Physics: Condensed Matter 25, 395305 (2013).
  * Wu  _et al._ [2014] S.-C. Wu, G. Shan, and B. Yan, Prediction of near-room-temperature quantum anomalous hall effect on honeycomb materials, Phys. Rev. Lett. 113, 256401 (2014).
  * Liu  _et al._ [2011] C.-C. Liu, W. Feng, and Y. Yao, Quantum spin hall effect in silicene and two-dimensional germanium, Phys. Rev. Lett. 107, 076802 (2011).
  * Xu  _et al._ [2012] C. Xu, G. Luo, Q. Liu, J. Zheng, Z. Zhang, S. Nagase, Z. Gao, and J. Lu, Giant magnetoresistance in silicene nanoribbons, Nanoscale 4, 3111 (2012).
  * Drummond  _et al._ [2012] N. D. Drummond, V. Zólyomi, and V. I. Fal’ko, Electrically tunable band gap in silicene, Phys. Rev. B 85, 075423 (2012).
  * Shan  _et al._ [2023] G. Shan, H. Tan, R. Ma, H. Zhao, and W. Huang, Recent progress in emergent two-dimensional silicene, Nanoscale 15, 2982 (2023).
  * Novoselov  _et al._ [2019] K. S. Novoselov, D. V. Andreeva, W. Ren, and G. Shan, Graphene and other two-dimensional materials, Frontiers of Physics 14, 13301 (2019).
  * Pflugradt  _et al._ [2014] P. Pflugradt, L. Matthes, and F. Bechstedt, Unexpected symmetry and aa stacking of bilayer silicene on $\mathrm{Ag}(111)$, Phys. Rev. B 89, 205428 (2014).
  * Sakai and Oshiyama [2015] Y. Sakai and A. Oshiyama, Structural stability and energy-gap modulation through atomic protrusion in freestanding bilayer silicene, Phys. Rev. B 91, 201405(R) (2015).
  * Ouyang  _et al._ [2018] X.-F. Ouyang, Z.-Y. Song, and Y.-Z. Zhang, Fully spin-polarized current in gated bilayer silicene, Phys. Rev. B 98, 075435 (2018).
  * Shih  _et al._ [2019] P.-H. Shih, T.-N. Do, G. Gumbs, H. D. Pham, and M.-F. Lin, Electric-field-diversified optical properties of bilayer silicene, Opt. Lett. 44, 4721 (2019).
  * Liu  _et al._ [2013] F. Liu, C.-C. Liu, K. Wu, F. Yang, and Y. Yao, $d+i{d}^{{}^{\prime}}$ chiral superconductivity in bilayer silicene, Phys. Rev. Lett. 111, 066804 (2013).
  * Zhang  _et al._ [2015] M.-M. Zhang, L. Xu, and J. Zhang, Topological insulator state in gated bilayer silicene, Journal of Physics: Condensed Matter 27, 445301 (2015).
  * Do  _et al._ [2018] T.-N. Do, P.-H. Shih, G. Gumbs, D. Huang, C.-W. Chiu, and M.-F. Lin, Diverse magnetic quantization in bilayer silicene, Phys. Rev. B 97, 125416 (2018).
  * Do  _et al._ [2019] T.-N. Do, G. Gumbs, P.-H. Shih, D. Huang, and M.-F. Lin, Valley- and spin-dependent quantum hall states in bilayer silicene, Phys. Rev. B 100, 155403 (2019).
  * spe [2018] _Silicene: Structure, properties and applications_ , 1st ed., Springer Series in Materials Science (SPRINGER, 2018).
  * Zhou  _et al._ [2018] Y. Zhou, Z.-X. Guo, S.-Y. Chen, H.-J. Xiang, and X.-G. Gong, Anisotropic in-plane thermal conductivity in multilayer silicene, Physics Letters A 382, 1499 (2018).
  * Gu  _et al._ [2018] X. Gu, Y. Wei, X. Yin, B. Li, and R. Yang, Colloquium: Phononic thermal properties of two-dimensional materials, Rev. Mod. Phys. 90, 041002 (2018).
  * Xiong  _et al._ [2010] D. Xiong, J. Wang, Y. Zhang, and H. Zhao, Heat conduction in two-dimensional disk models, Phys. Rev. E 82, 030101(R) (2010).
  * Wang  _et al._ [2012] L. Wang, B. Hu, and B. Li, Logarithmic divergent thermal conductivity in two-dimensional nonlinear lattices, Phys. Rev. E 86, 040101(R) (2012).
  * Pei  _et al._ [2014] Q.-X. Pei, Z.-D. Sha, Y.-Y. Zhang, and Y.-W. Zhang, Effects of temperature and strain rate on the mechanical properties of silicene, Journal of Applied Physics 115, 10.1063/1.4861736 (2014), 023519.
  * Xie  _et al._ [2016] H. Xie, T. Ouyang, E. Germaneau, G. Qin, M. Hu, and H. Bao, Large tunability of lattice thermal conductivity of monolayer silicene via mechanical strain, Phys. Rev. B 93, 075404 (2016).
  * Kuang  _et al._ [2016] Y. D. Kuang, L. Lindsay, S. Q. Shi, and G. P. Zheng, Tensile strains give rise to strong size effects for thermal conductivities of silicene, germanene and stanene, Nanoscale 8, 3760 (2016).
  * peng Li and qin Zhang [2012] H. peng Li and R. qin Zhang, Vacancy-defect–induced diminution of thermal conductivity in silicene, Europhysics Letters 99, 36001 (2012).
  * Wang  _et al._ [2015] Z. Wang, T. Feng, and X. Ruan, Thermal conductivity and spectral phonon properties of freestanding and supported silicene, Journal of Applied Physics 117, 084317 (2015), https://doi.org/10.1063/1.4913600 .
  * Xie  _et al._ [2014] H. Xie, M. Hu, and H. Bao, Thermal conductivity of silicene from first-principles, Applied Physics Letters 104, 10.1063/1.4870586 (2014), 131906.
  * Berdiyorov  _et al._ [2014] G. R. Berdiyorov, M. Neek-Amal, F. M. Peeters, and A. C. T. van Duin, Stabilized silicene within bilayer graphene: A proposal based on molecular dynamics and density-functional tight-binding calculations, Phys. Rev. B 89, 024107 (2014).
  * Tersoff [1989] J. Tersoff, Modeling solid-state chemistry: Interatomic potentials for multicomponent systems, Phys. Rev. B 39, 5566 (1989).
  * Hu  _et al._ [2013] M. Hu, X. Zhang, and D. Poulikakos, Anomalous thermal response of silicene to uniaxial stretching, Phys. Rev. B 87, 195417 (2013).
  * Pei  _et al._ [2013] Q.-X. Pei, Y.-W. Zhang, Z.-D. Sha, and V. B. Shenoy, Tuning the thermal conductivity of silicene with tensile strain and isotopic doping: A molecular dynamics study, Journal of Applied Physics 114, 10.1063/1.4815960 (2013), 033526.
  * Zhang  _et al._ [2014] X. Zhang, H. Xie, M. Hu, H. Bao, S. Yue, G. Qin, and G. Su, Thermal conductivity of silicene calculated using an optimized stillinger-weber potential, Phys. Rev. B 89, 054310 (2014).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. Fabris, G. Fratesi, S. de Gironcoli, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, J. Phys. Condens. Matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. D. Jr, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Kucukbenli, M. .Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Ponce, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with QUANTUM ESPRESSO : a modular and open-source software project for quantum simulations of materials, J. Phys. Condens. Matter 29, 465901 (2017).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * [40] See supplemental material at [url] for (s1) comments on the suitability of the refined lj model for phonon properties (s2) the effect of temperature on structural configuration, (s3) the vacf and the variation of the peak of phonon dos with temperature, and (s4) phonon dispersion and thermodynamic properties, (s5) phonon dispersion and group velocity, and (s6) note on vdw corrections.
  * Wang and Wu [2017] X. Wang and Z. Wu, Intrinsic magnetism and spontaneous band gap opening in bilayer silicene and germanene, Phys. Chem. Chem. Phys. 19, 2148 (2017).
  * Fu  _et al._ [2014] H. Fu, J. Zhang, Z. Ding, H. Li, and S. Meng, Stacking-dependent electronic structure of bilayer silicene, Applied Physics Letters 104, 10.1063/1.4870534 (2014), 131904.
  * Padilha and Pontes [2015] J. E. Padilha and R. B. Pontes, Free-standing bilayer silicene: The effect of stacking order on the structural, electronic, and transport properties, The Journal of Physical Chemistry C 119, 3818 (2015), https://doi.org/10.1021/jp512489m .
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Special points for brillouin-zone integrations, Phys. Rev. B 13, 5188 (1976).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comp. Phys. Comm. 271, 108171 (2022).
  * Stillinger and Weber [1985] F. H. Stillinger and T. A. Weber, Computer simulation of local order in condensed phases of silicon, Phys. Rev. B 31, 5262 (1985).
  * Hestenes and Stiefel [1952] M. R. Hestenes and E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Res. Natl. Inst. Stand. Technol 49, 409 (1952).
  * Nosé [1984] S. Nosé, A molecular dynamics method for simulations in the canonical ensemble, Molecular Physics 52, 255 (1984), https://doi.org/10.1080/00268978400101201 .
  * Hoover [1985] W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A 31, 1695 (1985).
  * Nosé and Klein [1983] S. Nosé and M. Klein, Constant pressure molecular dynamics for molecular systems, Molecular Physics 50, 1055 (1983), https://doi.org/10.1080/00268978300102851 .
  * Verlet [1967] L. Verlet, Computer ”experiments” on classical fluids. i. thermodynamical properties of lennard-jones molecules, Phys. Rev. 159, 98 (1967).
  * Kong  _et al._ [2011] L. T. Kong, C. Denniston, and M. H. Müser, An improved version of the green’s function molecular dynamics method, Computer Physics Communications 182, 540 (2011).
  * Kong [2011] L. T. Kong, Phonon dispersion measured directly from molecular dynamics simulations, Computer Physics Communications 182, 2201 (2011).
  * Kong  _et al._ [2009a] L. T. Kong, G. Bartels, C. Campañá, C. Denniston, and M. H. Müser, Implementation of green’s function molecular dynamics: An extension to lammps, Computer Physics Communications 180, 1004 (2009a).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scr. Mater. 108, 1 (2015).
  * Togo [2023] A. Togo, First-principles phonon calculations with phonopy and phono3py, J. Phys. Soc. Jpn. 92, 012001 (2023).
  * Wen  _et al._ [2018] M. Wen, S. Carr, S. Fang, E. Kaxiras, and E. B. Tadmor, Dihedral-angle-corrected registry-dependent interlayer potential for multilayer graphene structures, Phys. Rev. B 98, 235404 (2018).
  * Naik  _et al._ [2019] M. H. Naik, I. Maity, P. K. Maiti, and M. Jain, Kolmogorov–crespi potential for multilayer transition-metal dichalcogenides: Capturing structural transformations in moiré superlattices, The Journal of Physical Chemistry C 123, 9770 (2019), https://doi.org/10.1021/acs.jpcc.8b10392 .
  * Leven  _et al._ [2016] I. Leven, T. Maaravi, I. Azuri, L. Kronik, and O. Hod, Interlayer potential for graphene/h-bn heterostructures, Journal of Chemical Theory and Computation 12, 2896 (2016), pMID: 27168429, https://doi.org/10.1021/acs.jctc.6b00147 .
  * Rappe  _et al._ [1992] A. K. Rappe, C. J. Casewit, K. S. Colwell, W. A. Goddard, and W. M. Skiff, Uff, a full periodic table force field for molecular mechanics and molecular dynamics simulations, Journal of the American Chemical Society 114, 10024 (1992), https://doi.org/10.1021/ja00051a040 .
  * Nelder and Mead [1965] J. A. Nelder and R. Mead, A Simplex Method for Function Minimization, The Computer Journal 7, 308 (1965), https://academic.oup.com/comjnl/article-pdf/7/4/308/1013182/7-4-308.pdf .
  * Silva  _et al._ [2020] A. Silva, V. E. P. Claerbout, T. Polcar, D. Kramer, and P. Nicolini, Exploring the stability of twisted van der waals heterostructures, ACS Applied Materials & Interfaces 12, 45214 (2020), pMID: 32894936, https://doi.org/10.1021/acsami.0c13971 .
  * Herrero and Ramírez [2020] C. P. Herrero and R. Ramírez, Thermodynamic properties of graphene bilayers, Phys. Rev. B 101, 035405 (2020).
  * Pozzo  _et al._ [2011] M. Pozzo, D. Alfè, P. Lacovig, P. Hofmann, S. Lizzit, and A. Baraldi, Thermal expansion of supported and freestanding graphene: Lattice constant versus interatomic distance, Phys. Rev. Lett. 106, 135501 (2011).
  * Magnin  _et al._ [2014] Y. Magnin, G. D. Förster, F. Rabilloud, F. Calvo, A. Zappelli, and C. Bichara, Thermal expansion of free-standing graphene: benchmarking semi-empirical potentials, Journal of Physics: Condensed Matter 26, 185401 (2014).
  * Lindsay and Broido [2010] L. Lindsay and D. A. Broido, Optimized tersoff and brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene, Phys. Rev. B 81, 205441 (2010).
  * Li  _et al._ [2018] G. Li, L. Zhang, W. Xu, J. Pan, S. Song, Y. Zhang, H. Zhou, Y. Wang, L. Bao, Y.-Y. Zhang, S. Du, M. Ouyang, S. T. Pantelides, and H.-J. Gao, Stable silicene in graphene/silicene van der waals heterostructures, Advanced Materials 30, 1804650 (2018).
  * Liu  _et al._ [2014] H. Liu, N. Han, and J. Zhao, Band gap opening in bilayer silicene by alkali metal intercalation, Journal of Physics: Condensed Matter 26, 475303 (2014).
  * Kong  _et al._ [2009b] B. D. Kong, S. Paul, M. B. Nardelli, and K. W. Kim, First-principles analysis of lattice thermal conductivity in monolayer and bilayer graphene, Phys. Rev. B 80, 033406 (2009b).
  * Ouyang  _et al._ [2020] W. Ouyang, H. Qin, M. Urbakh, and O. Hod, Controllable thermal conductivity in twisted homogeneous interfaces of graphene and hexagonal boron nitride, Nano Letters 20, 7513 (2020), pMID: 32898421, https://doi.org/10.1021/acs.nanolett.0c02983 .
  * Tretiakov and Scandolo [2004] K. V. Tretiakov and S. Scandolo, Thermal conductivity of solid argon from molecular dynamics simulations, The Journal of Chemical Physics 120, 3765 (2004), https://pubs.aip.org/aip/jcp/article-pdf/120/8/3765/12243463/3765_1_online.pdf .
  * Fan  _et al._ [2017] Z. Fan, L. F. C. Pereira, P. Hirvonen, M. M. Ervasti, K. R. Elder, D. Donadio, T. Ala-Nissila, and A. Harju, Thermal conductivity decomposition in two-dimensional materials: Application to graphene, Phys. Rev. B 95, 144309 (2017).
  * Chakraborty  _et al._ [2018] H. Chakraborty, S. Mogurampelly, V. K. Yadav, U. V. Waghmare, and M. L. Klein, Phonons and thermal conducting properties of borocarbonitride (bcn) nanosheets, Nanoscale 10, 22148 (2018).
  * Ong and Pop [2011] Z.-Y. Ong and E. Pop, Effect of substrate modes on thermal transport in supported graphene, Phys. Rev. B 84, 075471 (2011).
  * Guo  _et al._ [2011] Z.-X. Guo, D. Zhang, and X.-G. Gong, Manipulating thermal conductivity through substrate coupling, Phys. Rev. B 84, 075470 (2011).


# Novel Fine-Tuned Attribute Weighted Naïve Bayes NLoS Classifier for UWB
Positioning

Fuhu Che⋆, Qasim Zeeshan Ahmed⋆, Fahd Ahmed Khan†, and Faheem A. Khan⋆ ⋆F.
Che, Q. Z. Ahmed, and F. A. Khan is with the School of Computing and
Engineering at the University of Huddersfield, United Kingdom. Email:
{fuhu.che, q.ahmed<EMAIL_ADDRESS>†F. A. Khan is with the School of
Electrical and Computer Engineering, University of Oklahoma, Tulsa, OK, USA.
Email<EMAIL_ADDRESS>Part of this work was presented at the IEEE
International Black Sea Conference on Communications and Networking 2022,
Sofia, Bulgaria.

###### Abstract

In this paper, we propose a novel Fine-Tuned attribute Weighted Naïve Bayes
(FT-WNB) classifier to identify the Line-of-Sight (LoS) and Non-Line-of-Sight
(NLoS) for UltraWide Bandwidth (UWB) signals in an Indoor Positioning System
(IPS). The FT-WNB classifier assigns each signal feature a specific weight and
fine-tunes its probabilities to address the mismatch between the predicted and
actual class. The performance of the FT-WNB classifier is compared with the
state-of-the-art Machine Learning (ML) classifiers such as minimum Redundancy
Maximum Relevance (mRMR)- $k$-Nearest Neighbour (KNN), Support Vector Machine
(SVM), Decision Tree (DT), Naïve Bayes (NB), and Neural Network (NN). It is
demonstrated that the proposed classifier outperforms other algorithms by
achieving a high NLoS classification accuracy of $99.7\%$ with imbalanced data
and $99.8\%$ with balanced data. The experimental results indicate that our
proposed FT-WNB classifier significantly outperforms the existing state-of-
the-art ML methods for LoS and NLoS signals in IPS in the considered scenario.

###### Index Terms:

UWB, ML, WNB, mRMR, and IPS.

##  I Introduction

Mobile target tracking is extensively used in commercial, civil, industrial,
and military applications, such as object tracking, autonomous vehicles,
forklift control, military reconnaissance, etc [1, 2, 3]. Ultra-wideband (UWB)
technology, utilizes exploratory signals with an extremely large spectrum,
that provides significant advantages in tracking and positioning [4, 5, 6, 7].
Beside this various alternative technologies have been extensively researched
and employed in indoor scenarios such as Radio Frequency Identification
(RFID), Bluetooth, WiFi fingerprinting, etc. Recently, [9, 3] has shown that
WiFi with channel state information can be utilized for Indoor Positioning
System (IPS) and have the same positioning cost and its accuracy is not
inferior to UWB. Non-line-of-sight (NLoS) conditions are a major hindrance in
achieving high positioning accuracy [10, 11]. In NLoS, the direct signal path
between the transmitters and the receivers is blocked by the obstacles causing
estimation error [8, 7]. To deal with the impact of NLoS, several signal-
feature-based Machine Learning (ML) classification methods have been proposed
[12, 15, 16, 11].

Non-parametric regression technique based on Support Vector Machine (SVM) and
Gaussian process to estimate the ranging error was proposed in [4]. It was
shown that for signals under NLoS conditions, $28\%$ had a ranging error of
less than one meter. Deep learning methods such as Convolutional Neural
Network (CNN) and Long Short-Term Memory (LSTM) were employed in [12, 15]. The
highest classification accuracy of $82.14\%$ was achieved by employing CNN and
stacked LSTM together [12]. Other ML techniques such as SVM, Random Forest
(RF), and Multi-Layer Perceptron (MLP) were applied in [16] and an accuracy of
$82.80\%$, $91.9\%$, and $91.20\%$ was reported for these techniques,
respectively. A semi-supervised learning approach has also been adopted for
NLoS identification. It was shown that NLoS identification accuracy could be
improved to $94.7\%$ with the aid of semi-supervised SVM via leveraging
unlabeled data [11]. An unsupervised ML approach based on Gaussian mixture
models (EM-GMM) to identify the NLoS links from the unlabelled data was
proposed in [17]. [18] employed the ranging and received signal strength (RSS)
to identify the NLoS condition using LoS, NLoS soft, and NLoS hard
classification. Furthermore, in [19] the power delay profile was exploited to
accelerate the training of NN-based classifiers. However, the dataset was
perfectly balanced. Finally, in [20], a different scenario was considered for
testing with respect to the training for generalization of IPS.

Feature-based ML methods were well-studied in the aforementioned works,
however, none of these works considered the correlation between different
features of the signals. Moreover, in an imbalanced dataset, where a limited
amount of NLoS samples are present, it is difficult for existing algorithms to
train a robust model [12]. To address these shortcomings, we propose a novel
feature-based Fine-Tuned attribute Weighted Naïve Bayes (FT-WNB) classifier.
To the best of the authors’ knowledge, no work has proposed signal-feature
weighting to enhance NLoS classification with an imbalanced dataset. The main
contributions of this work are as follows:

  * •

A novel signal-feature-based FT-WNB algorithm for NLoS classification
addresses class imbalance by assigning different weights to each feature and
fine-tunes the classifier probabilities to improve classification accuracy.

  * •

We conduct extensive experimentation to validate the superiority of the
proposed FT-WNB classifier in comparison to the state-of-art algorithms. The
performance is compared in terms of the confusion matrix, Receiver Operating
Characteristics (ROC) curve, Area Under Curve (AUC), precision, recall, and
classification accuracy.

  * •

We show that the FT-WNB classification algorithm outperforms existing
algorithms for different ratios of NLoS and LoS samples, even in the presence
of imbalanced dataset.

##  II Fine-Tuned Weighted Naïve Bayes (FT-WNB) Classification Algorithm

###  II-A UWB Range Theory

UWB transmission and detection algorithms can be found in more details in [13,
14]. Time-of-Arrival (ToA) is used to calculate the estimated distance $d$
between the tag and the anchor [15],

$d=\begin{cases}c\times\mathrm{ToA}+\epsilon,&\ \ \mathrm{LoS}\\\
c\times\mathrm{ToA}+\epsilon+b,&\ \ \mathrm{NLoS}\\\ \end{cases}$ (1)

where, $c$ is the speed of light in $\mathrm{m/s}$, $\mathrm{ToA}$ is the
propagation time taken from the tag to the anchor, $\epsilon$ represents the
measurement noise and $b$ is the NLoS measurement error caused by the
blockages and obstacles.

###  II-B Weighted Naïve Bayes

Naïve Bayes (NB) classifies the data with a corresponding maximum posterior
probability. As features are assumed to be conditionally independent, the
predicted class $\hat{l}^{k}$ at the $k$-th instance of data can be determined
as

$\displaystyle\hat{l}^{k}=\mathop{\arg\max}_{l}P(l)\prod_{i=1}^{I}P(X_{i}^{k}|l),$
(2)

where $P(X_{i})$ is the probability of the $i$-th attribute, $P(l)$ is the
apriori probability of class $l$ and $I$ is the number of attributes. The
underlying principle of WNB algorithm is that some attributes are of more
importance than others [21]. This importance is specified by weighting the
attributes and the predicted class, $\hat{l}^{k}$, in this case, is determined
as

$\displaystyle\hat{l}^{k}=\mathop{\arg\max}_{l}P(l)\prod_{i=1}^{I}P({X_{i}^{k}|l})^{w(i)},$
(3)

where $w(i)$ is the weight associated with the $i$-th attribute.

.

Figure 1: Block diagram of Fine-tuned WNB algorithm. Figure 2: Experimental
data collected in the studio.

###  II-C FT-WNB Algorithm

Figure 3: Signal-features correlation matrix.

The overall framework of the FT-WNB algorithm is described in Fig. 1. We first
initialize the conditional probabilities as defined in (3) by applying the
attribute weights to predict the class label. In the training phase, if the
training instance is classified incorrectly, the predicted class $l_{P}$ has a
higher posterior probability than its true class $l_{T}$. Therefore, the
conditional probabilities of signal features, $P(X_{i}^{k}|l_{T})$ and
$P(X_{i}^{k}|l_{P})$ can be fine tuned in the next iteration as

$\displaystyle P_{t+1}(X_{i}^{k}|l_{T})$ $\displaystyle=$ $\displaystyle
P_{t}(X_{i}^{k}|l_{T})+\xi_{t+1}(X_{i}^{k},l_{T}),$ (4) $\displaystyle
P_{t+1}(X_{i}^{k}|l_{P})$ $\displaystyle=$ $\displaystyle
P_{t}(X_{i}^{k}|l_{P})-\xi_{t+1}(X_{i}^{k},l_{P}),$ (5)

where $t$ represents the number of iteration and $\xi_{t+1}(\cdot)$ are the
updated step size calculated as

$\displaystyle\\!\\!\\!\\!\xi_{t+1}(X_{i}^{k},l_{T})\\!\\!\\!\\!\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\\!\beta\\!\\!\left(\alpha\\!\\!\left(\max
P_{t}(X_{i}^{k}|l_{T})\right)\\!\\!-\\!\\!P(X_{i}^{k}|l_{T})\right)\\!\\!e(l_{T},l_{P}),$
(6) $\displaystyle\\!\\!\\!\\!\xi_{t+1}(X_{i}^{k},l_{P})\\!\\!\\!\\!\\!\\!$
$\displaystyle=$ $\displaystyle\\!\\!\\!\\!\\!\\!\beta\\!\\!\left(\alpha
P(X_{i}^{k}|l_{P})\\!\\!-\\!\\!\min\left(P_{t}(X_{i}^{k}|l_{P})\right)\right)\\!\\!e(l_{T},l_{P}),$
(7)

where $\alpha,\beta\in(0,1)$ controls the size of the updated step for
$P_{t}(X_{i}^{k}|l_{T})$ and $P_{t}(X_{i}^{k}|l_{P})$, respectively.
$\max(P_{t}(X_{i}^{k}|l_{T}))$ and $\min(P_{t}(X_{i}^{k}|l_{P}))$ are the
maximum and the minimum value for $P_{t}(X_{i}^{k}|l_{T})$ and
$P_{t}(X_{i}^{k}|l_{P})$, respectively. The error term, $e(l_{T},l_{P})$, is
calculated as

$\displaystyle e(l_{T},l_{P})=|P(l_{T})-P(l_{P})|.$ (8)

Finally, the steps involved in the training and testing phase are summarised
in Algorithms $1$ and $2$, respectively.

Algorithm 1 : Training Phase

Training dataset consisting of signal features.

The conditional probability vector

Phase 1: Initializing conditional probabilities

$\square\quad\textbf{{Step 1:}}$

for each $X_{i}^{k}$ do

calculate $\hat{l}^{k}$ by (2)

end for

for each weighted $X_{i}^{k}$ do

calculate $\hat{l}^{k}$ by (3)

end for Estimate the conditional probability vector

Phase 2: Fine tuning conditional probabilities

$\square\quad\textbf{{Step 2:}}$

if $l_{T}\neq l_{P}$ then

Calculate the error by (8)

Calculate the $\max(P_{t}(X_{i}^{k}|l_{T}))$ and
$\min(P_{t}(X_{i}^{k}|l_{P}))$

for each weighted figure $X_{i}^{k}$ do

calculate $\xi_{t+1}(X_{i}^{k},l_{T})$ by (6)

update $P_{t+1}(X_{i}^{k}|l_{T})$ by (4)

calculate $\xi_{t+1}(X_{i}^{k},l_{P})$ by (7)

update $P_{t+1}(X_{i}^{k}|l_{P})$ by (5)

end for

end if

if $l_{T}=l_{P}$ then

end if

return Classification model

Algorithm 2 : Testing Phase

Test dataset with a mixture of LoS and NLoS signals

Predicted class label of $X_{i}^{k}$

Algorithm:

  1. 1.

Estimate the class membership probability $P(X_{i}^{k}|l)$.

  2. 2.

Predict the class label $\hat{l}^{k}$ of each class $X_{i}^{k}$.

  3. 3.

Return $\hat{l}^{k}$.

##  III Experimental Setup and Data Collection

The required dataset was collected in a studio environment ($4.8m\times
3.5m$). The floorplan of the infrastructure is shown in Fig. 2. DecawaveⓇ
MDEK-$1001$ UWB kits, based on DW-1000 chip, with configuration as mentioned
in [22], are used to generate the dataset. During data collection for the LoS
scenario, there is no obstacle between the anchor and tag. For the NLoS
scenario, an iron sheet is placed between the anchors and the tag. The plates
are placed very close to the tag, generally less than a meter away. In this
case, the first path of the signal transmission is either completely blocked
or distorted through multiple reflections. As a result of this, the estimated
range will be biased due to the signal propagation delay. From the collected
dataset, we randomly select $100$ NLoS and $1000$ LoS signals to generate data
imbalance. We consider two more features in addition to the $10$ features of
[23]: i) the index of the detected first-path ($FP\\_index$) and ii) the power
ratio between the estimated received power and first-path power. The features
and the correlation between them are illustrated in Fig. 3. As expected, the
features are clearly not equally important nor independent. It can be observed
that the features amplitude of the first path ($F1$), the amplitude of the
second path ($F2$), and the amplitude of the third path ($F3$), as well as
first path power level ($FPPL$), have a relatively higher correlation than
others.

Figure 4: Samples of Error vs. Power Difference

##  IV Performance Evaluation

Visualization of our samples is shown in Fig. 4 where the power difference
versus the distance error is plotted. There are $1000$ LoS samples and $100$
NLoS signals, respectively. In this figure, the green samples represent the
LoS signals and the blue color represents the NLoS signals. From the figure,
it can be observed that few NLoS samples can be easily differentiated from LoS
samples. While a large number of samples for NLoS cannot be easily
differentiated from LoS samples. Therefore, for accurate classification, an
appropriate algorithm is required.

Figure 5: Performance comparison for the proposed FT-WNB algorithm with state-
of-the-art ML algorithms. TABLE I: Running time and Confusion matrix of the
proposed FT-WNB algorithm and state-of-the-art ML algorithms

| Running | LoS | NLoS |  |  |  |  |
---|---|---|---|---|---|---|---|---
Algorithms | Time (s) | CR | CR | TP | FN | FP | TN | Accuracy($\%$)
mRMR-KNN | 0.049 | 98.6% | 87% | 986 | 14 | 13 | 87 | 97.5
mRMR-SVM | 0.092 | 98.2% | 90% | 982 | 18 | 10 | 90 | 97.5
mRMR - DT | 0.104 | 98.4% | 89% | 984 | 16 | 11 | 89 | 97.5
mRMR - NB | 0.048 | 98.9% | 90% | 989 | 11 | 10 | 90 | 98.1
mRMR - NN | 0.061 | 98.8% | 92% | 988 | 12 | 8 | 92 | 98.2
FT-WNB | 0.041 | 99.5% | 98% | 995 | 5 | 2 | 98 | 99.4

The performance of the proposed FT-WNB classifier is compared with ML
classifiers in Fig. 5. For these classifier, a subset of features is selected
based on minimum redundancy maximum relevance (mRMR) feature selection. mRMR
is a feature selection approach that tends to select features with a high
correlation with the class and a low correlation between themselves to
maximize relevance toward the features and simultaneously to minimize the
redundancy among features [24]. The running time, confusion matrix and the
correct rate (CR) is summarised in Table I. The confusion matrix depicts
metrics, True Positive (TP), False Positive (FP), False Negative (FN) and True
Negative (TN), respectively. In order to compare the classification
performance, we plot precision, recall and classification accuracy in Fig. 5.
It can be observed that the best performing algorithm among the existing ML
algorithms is the mRMR-NN, which has a precision of approximately $99.2\%$,
recall around $98.8\%$ and overall classification accuracy is $98.2\%$. mRMR-
NB follows closely and has similar performance. However, it is evident from
Fig. 5, that the proposed FT-WNB algorithm outperforms all the state-of-the-
art ML algorithms, i.e., mRMR-KNN, -SVM, -DT, -NB and -NN by achieving a
precision of $99.7\%$, recall of $99.5\%$ and an overall classification
accuracy of $99.4\%$. From Table I it can be observed that the proposed FT-WNB
algorithm, TN$=98$ which indicates that $98\%$ NLoS samples are correctly
classified and only $2\%$ NLoS samples are incorrectly classified. The
classification performance is significantly better than the second best case
which is $92\%$ achieved by mRMR-NN algorithm. TP$=995$, refers to LoS
resulting in a correct rate of $99.5\%$ which means $5$ samples out of $1000$
samples were inaccurately classified. Again, it is $1.3\%$ higher than the
second-best case. The average running time of FT-WNB is $0.041$s which is
better than other considered algorithms where the process of mRMR takes
significant time. To sum up, in comparison to mRMR-KNN, mRMR-SVM, mRMR-DT,
mRMR-NB, and mRMR-NN ML algorithms, it can be observed that the FT-WNB
algorithm performs better for all performance metrics of the confusion matrix.

Figure 6: Classification accuracy of the FT-WNB algorithm under imbalanced
datasets having varying ratios of data. Figure 7: Comparisons for the existing
UWB position methods.

Fig. 6 shows the overall classification accuracy in the case of an imbalanced
dataset. The ratio (NLoS: LoS) was set to $0.1$, $0.5$, and $1.0$,
respectively. It can be observed from Fig. 6 that as the ratio increases, the
results of all algorithms improve because the model is trained better with
equal number of LoS and NLoS samples. However, it is noticeable that the FT-
WNB algorithm always maintains a higher classification accuracy for all
considered ratios. This shows that the proposed FT-WNB classifier is highly
robust and can assign appropriate weights to the desired features,
guaranteeing better NLoS identification as compared to existing ML algorithms.

We also compares the proposed algorithm with the existing UWB NLoS
identification methods as mentioned in the introduction section. The results
in terms of accuracy are shown in Fig. 7. It can be clearly observed that the
proposed FT-WNB algorithm performs better than MLP, EM-GMM, MLP, CNN-LSTM, and
CNN methods which are some of the existing positioning methods in UWB IPS.

Figure 8: Impact of number of samples evaluation

For the fine-tuning step, the impact of a number of selected samples is shown
in Fig. 8. It can be observed that the accuracy of the algorithm improves as
the number of samples increase and the proposed model benefits from the fine-
tuning step. The classification accuracy improves from $98.2\%$ to $98.6\%$ by
collecting the first $10$ samples and then further increases to $99.5\%$ by
just collecting the first $40$ samples. However, there is no significant gain
in accuracy beyond that point. Therefore, in our simulations, the number of
samples required for fine-tuning is fixed at $40$.

TABLE II: Comparison with [18]. Algorithms | Processing |  $\%$ LoS |  $\%$ NLoS |  $\%$ Total
---|---|---|---|---
| Time (s) | Accuracy | Accuracy | Accuracy
Range & RSS | 0.0371 | 97.4 | 91 | 96.8
All features | 0.041 | 99.5 | 98 | 99.4
Difference | 3.9 ms | 2.1 ($\%$) | 7 ($\%$) | 2.6 ($\%$)

There is a tradeoff between the collected features and accuracy. In this
paper, our focus was to collect a maximum number of features that yield the
highest accuracy. In Table II, we present the comparison of the processing
time considering only two features (Range and RSS [18]) and all the $12$
gathered features. It can be observed that a difference of $3.9$ms leads to an
improved accuracy of $99.4\%$ as compared to $96.8\%$. However, if the
application does not require high accuracy, we could reduce the number of
collected features and improve the processing time of the algorithm.

Our algorithm has been tested and trained in various scenarios as adopted in
[20]. It can be observed from Table III that the proposed WNB has a
classification accuracy of over $98\%$, even after being trained in a studio
and tested in a room measuring $(4.8m\times 5.4m)$. This demonstrates the
versatility of FT-WNB and its ability to perform in different environments,
regardless of where it was trained.

TABLE III: Performance of FT-WNB in a different environment. Training | Testing | TP | FP | FN | TN | Accuracy
---|---|---|---|---|---|---
| Scenario | Scenario |  |  |  | ($\%$)
Studio | Studio | 995 | 5 | 2 | 98 | 99.4
Studio | Room | 990 | 10 | 5 | 95 | 98.6

##  V Conclusion

In this work, we propose the FT-WNB algorithm which is a feature-based NLoS
classification method to improve classification accuracy by automatically
weighting different signal features, for imbalanced LoS and NLoS samples. The
performance of the FT-WNB classifier is compared with state-of-the-art machine
learning algorithms using simulation. It is shown that the proposed FT-WNB
algorithm performed better in terms of confusion matrix, ROC, and AUC. In
addition, the proposed FT-WNB algorithm adapts better to the different ratios
of LoS and NLoS signal data. Our findings demonstrate that the proposed FT-WNB
classifier efficiently alleviates the influence of the signal features in an
imbalanced dataset in a harsh mixed LoS and NLoS scenario.

## References

  * [1] F. Zafari, A. Gkelias and K. K. Leung, “A Survey of Indoor Localization Systems and Technologies,” IEEE Communications Surveys & Tutorials, vol. 21, no. 3, pp. 2568-2599, third quarter 2019.
  * [2] Q. Z. Ahmed, K. -H. Park and M. -S. Alouini, “Ultrawide Bandwidth Receiver Based on a Multivariate Generalized Gaussian Distribution,” in IEEE Transactions on Wireless Communications, vol. 14, no. 4, pp. 1800-1810, April 2015.
  * [3] P. Roy and C. Chowdhury, “A survey on ubiquitous WiFi-based indoor localization system for smartphone users from implementation perspectives,” CCF Trans. on Pervasive Computing and Interaction, 1-21, 2022.
  * [4] H. Wymeersch et.al., “A Machine Learning Approach to Ranging Error Mitigation for UWB Localization,” IEEE Transactions on Communications, vol. 60, no. 6, pp. 1719-1728, June 2012.
  * [5] Q. Z. Ahmed and L. -L. Yang, “Reduced-rank adaptive multiuser detection in hybrid direct-sequence time-hopping ultrawide bandwidth systems,” in IEEE Transactions on Wireless Communications, vol. 9, no. 1, pp. 156-167, January 2010.
  * [6] Q. Z. Ahmed, L. -L. Yang and S. Chen, “Reduced-Rank Adaptive Least Bit-Error-Rate Detection in Hybrid Direct-Sequence Time-Hopping Ultrawide Bandwidth Systems,” in IEEE Transactions on Vehicular Technology, vol. 60, no. 3, pp. 849-857, March 2011.
  * [7] F. Che et al., “Feature-Based Generalized Gaussian Distribution Method for NLoS Detection in Ultra-Wideband (UWB) Indoor Positioning System,” IEEE Sensors Journal, vol. 22, no. 19, pp. 18726-18739, 2022.
  * [8] F. Che, A. Ahmed, Q. Z. Ahmed, S. A. R. Zaidi and M. Z. Shakir, “Machine Learning Based Approach for Indoor Localization Using Ultra-Wide Bandwidth (UWB) System for Industrial Internet of Things (IIoT),” 2020 International Conference on UK-China Emerging Technologies (UCET), Glasgow, UK, pp. 1-4, 2020.
  * [9] Z. Zhang, H. Du, S. Choi and S. H. Cho, “TIPS: Transformer Based Indoor Positioning System Using Both CSI and DoA of WiFi Signal,” IEEE Access, vol. 10, pp. 111363-111376, 2022.
  * [10] W. Wang, D. Marelli and M. Fu, “Multiple-Vehicle Localization Using Maximum Likelihood Kalman Filtering and Ultra-Wideband Signals,” IEEE Sensors Journal, vol. 21, no. 4, pp. 4949-4956, 2021.
  * [11] T. Wang, K. Hu, Z. Li, K. Lin, J. Wang and Y. Shen, ”A Semi-Supervised Learning Approach for UWB Ranging Error Mitigation,” IEEE Wireless Communications Letters, vol. 10, no. 3, pp. 688-691, March 2021.
  * [12] C. Jiang, J. Shen, S. Chen, Y. Chen, D. Liu and Y. Bo, ”UWB NLOS/LOS Classification Using Deep Learning Method,” IEEE Communications Letters, vol. 24, no. 10, pp. 2226-2230, Oct. 2020.
  * [13] Q. Z. Ahmed, W. Liu and L. -L. Yang, ”Least Mean Square Aided Adaptive Detection in Hybrid Direct-Sequence Time-Hopping Ultrawide Bandwidth Systems,” VTC Spring 2008 - IEEE Vehicular Technology Conference, Marina Bay, Singapore, pp. 1062-1066, 2008.
  * [14] Q. Z. Ahmed and L. -L. Yang, ”Normalised Least Mean-Square Aided Decision-Directed Adaptive Detection in Hybrid Direct-Sequence Time-Hopping UWB Systems,” 2008 IEEE 68th Vehicular Technology Conference, Calgary, AB, Canada, pp. 1-5, 2008.
  * [15] M. Stahlke, S. Kram, F. Ott, T. Feigl and C. Mutschler, “Estimating TOA Reliability With Variational Autoencoders,” IEEE Sensors Journal, vol. 22, no. 6, pp. 5133-5140, 2022.
  * [16] C.L. Sang et al., “Identification of NLOS and Multi-Path Conditions in UWB Localization Using Machine Learning Methods,” Applied Sciences, vol. 10, no. 11, pp. 3980-4005, June 2020.
  * [17] J. Fan and A. S. Awan, ”Non-Line-of-Sight Identification Based on Unsupervised Machine Learning in Ultra Wideband Systems,” IEEE Access, vol. 7, pp. 32464-32471, 2019.
  * [18] V.Barralet al., ”NLOS identification and mitigation using low-cost UWB devices.” Sensors, 19(16), 3464, 2019.
  * [19] V. B. Vales, T. Domínguez-Bolaño, C. J. Escudero and J. A. García-Naya, “Using the Power Delay Profile to Accelerate the Training of Neural Network-Based Classifiers for the Identification of LOS and NLOS UWB Propagation Conditions,” IEEE Access, vol. 8, pp. 220205-220214, 2020.
  * [20] V. Barral, C. J. Escudero, J. A. García-Naya and P. Suárez-Casal, “Environmental Cross-Validation of NLOS Machine Learning Classification/Mitigation with Low-Cost UWB Positioning Systems”, in Sensors 2019, 19, 5438.
  * [21] S. Ruan, H. Li, C. Li and K. Song, “Class-Specific Deep Feature Weighting for Naïve Bayes Text Classifiers,” IEEE Access, vol. 8, pp. 20151-20159, 2020.
  * [22] DW1000 User Manual. Available online: https://www.decawave.com/sites/default/files/resources
  * [23] F. Che et al., ”Weighted Naive Bayes Approach for Imbalanced Indoor Positioning System Using UWB,” IEEE International Black Sea Conference on Communications and Networking (BlackSeaCom), 2022, pp. 72-76.
  * [24] A. Samat, E. Li, P. Du, S. Liu and J. Xia, “GPU-Accelerated CatBoost-Forest for Hyperspectral Image Classification Via Parallelized mRMR Ensemble Subspace Feature Selection,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 14, pp. 3200-3214, 2021.