# Resource Constrained Neural Networks for
5G Direction-of-Arrival Estimation in Micro-controllers

Piyush Sahoo1, Romesh Rajoria1, Shivam Chandhok2, S. J. Darak1, Danilo Pau3,
Hem-Dutt Dabral4
1Algorithms to Architecture Lab, ECE, IIIT Delhi, India
2Mohamed bin Zayed University of AI, UAE,
3STMicroelectronics, Milan, Italy,
4STMicroelectronics, Greater Noida, India.

###### Abstract.

With the introduction of shared spectrum sensing and beam-forming based multi-
antenna transceivers, 5G networks demand spectrum sensing to identify
opportunities in time, frequency, and spatial domains. Narrow beam-forming
makes it difficult to have spatial sensing (direction-of-arrival, DoA,
estimation) in a centralized manner, and with the evolution of paradigms such
as artificial intelligence of Things (AIOT), ultra-reliable low latency
communication (URLLC) services and distributed networks, intelligence for edge
devices (Edge-AI) is highly desirable. It helps to reduce the data-
communication overhead compared to cloud-AI-centric networks and is more
secure and free from scalability limitations. However, achieving desired
functional accuracy is a challenge on edge devices such as microcontroller
units (MCU) due to area, memory, and power constraints. In this work, we
propose low complexity neural network-based algorithm for accurate DoA
estimation and its efficient mapping on the off-the-self MCUs. An ad-hoc
graphical-user interface (GUI) is developed to configure the STM32
NUCLEO-H743ZI2 MCU with the proposed algorithm and to validate its
functionality. The performance of the proposed algorithm is analyzed for
different signal-to-noise ratios (SNR), word-length, the number of antennas,
and DoA resolution. In-depth experimental results show that it outperforms the
conventional statistical spatial sensing approach.

DoA estimation, edge-AI, 5G, neural networks, X-CUBE-AI, micro-controller,
STM32

††conference: First International Conference on AI-ML Systems; Oct. 21–23,
2021; Bangalore, India

##  1\. Introduction

A never-ending expansion of network connectivity has led to the evolution of
the Internet of things (IoT) capable to route heterogeneous information
between people and devices on a planetary scale (Chang et al., 2021; Shi et
al., 2016; Chettri and Bera, 2020). From simple temperature and pressure
sensors to smart home and smart industry processing units relying on faster
networks, IoT has resulted in exponential growth, offering increased
productivity, better security, as well as the availability of a huge amount of
business-worth data (Chang et al., 2021; Shi et al., 2016; Chettri and Bera,
2020). The analysis of such data to get useful trends and patterns is itself a
challenging task and recent progress in Artificial Intelligence (AI) has led
to promising solutions for exploiting the hidden information and correlations
embodied into the data (Li et al., 2018; Al-Garadi et al., 2020; Zolanvari et
al., 2019). Conventional cloud-AI IoT networks perform data processing and
inferences at resource-rich data centers. Unfortunately, such a centralized
approach suffers from data-communication overhead, high unpredictable latency,
low security, privacy, and poor scalability (Verhelst and Moons, 2017; Lin et
al., 2020; Li et al., 2020; Cass, 2019). Furthermore, data centers are not
environmentally friendly, energy parsimonious, and green (Ferreira et al.,
2018).

To overcome some of these challenges, the Edge-AI distributed approach is
being considered as a solution in which edge devices are featuring
intelligence to perform inference either independently or in collaboration
with other processing units (Verhelst and Moons, 2017; Lin et al., 2020; Li et
al., 2020; Cass, 2019). Edge-AI is being actively explored for speech, video,
environmental, industry 4.0, and autonomous driving applications. For
instance, edge-AI includes assisting devices such as cognitive hearing aids
(Wang, 2017), optical character recognition (Memon et al., 2020; Amin and
Hossain, 2021), and smart home appliances such as google home, Alexa, etc.
With the introduction of shared spectrum sensing due to limited spectrum
availability, 5G networks demand spectrum sensing to identify spectrum
opportunities in time, frequency and spatial domains (Sharma et al., 2018;
Ahmad et al., 2020; Lien et al., 2017). The beam-forming based multi-antenna
transceivers allow spatially separated users to communicate simultaneously
over the common spectrum leading to significant improvement in network
throughput (Gao et al., 2016). However, narrow beam-forming makes it difficult
to have spatial sensing (direction-of-arrival, DoA, estimation) in a
centralized manner, and hence, intelligence at user devices is highly
desirable. Edge-AI is also important to bring next-generation computing
paradigms such as artificial intelligence of Things (AIOT), ultra-reliable low
latency communication (URLLC) services, smart industry 4.0, and distributed
networks to reality (Dohler et al., 2017; Wang et al., 2019).

In the last decade, machine and deep learning algorithms have significantly
outperformed hand-crafted rule-based and classic machine learning approaches
in various application domains. In wireless communications, deep learning
algorithms are widely used for spectrum sensing and characterization, channel
estimation, and resource allocations (Wang et al., 2020, 2017; Mao et al.,
2018). However, most of their deployment is limited to cloud-AI networks due
to the availability of virtually unlimited resources (memory, computational
power, etc.). The deployment of these algorithms in Edge-AI networks needs
careful design and low power mapping on hardware since achieving the desired
target accuracy is a challenge due to area, memory and power-constrained
micro-controller units (MCU) usually feature. The work presented in this paper
attempts to address these challenges for spatial sensing applications in
wireless networks by targeting the highest accuracy with the lowest
implementation costs and using the MCU resources parsimoniously.

Figure 1. Overall architecture for the proposed DoA estimation framework.

###  1.1. Problem Formulation

The goal of this work is to develop a lightweight DoA estimation framework
using artificial neural networks that can be easily deployed on off-the-shelf
resource-constrained MCUs. Various works have explored DoA estimation for
wireless signals and its implementation on hardware such as FPGA/ASICs (Joshi
et al., 2021; Liu et al., 2018; Hussain et al., 2018). To the best of our
knowledge, DoA estimation for wireless signals using neural networks for MCUs
has not been explored yet in the literature and is the focus of the proposed
work. The DoA estimation-based spatial sensing is also useful in radar, sonar,
and navigation to track an object, biomedical to detect tumors, and artery
wall movements.

We consider the analog-front-end with multiple antennas and assume that the
MCU receives a baseband signal consisting of one active transmission with
unknown DoA. The problem of DOA estimation can be formulated as a $K$ class
classification problem, where each class corresponds to a possible angle or
direction of signal arrival, depending on the antenna-array placement and
desired resolution. Given an angle range $A$ of antenna array and a specific
DoA resolution $R$, we define a set of DoA classes
$C=\\{c_{1},c_{2},c_{3}....c_{K}\\}$, where the number of classes, $K=A/R$. In
most of the cases, $A=360^{\circ}$. Figure 1 shows the overall architecture of
the proposed DoA estimation framework using a deep-learning-based supervised
learning paradigm. Formally, we first generate the training and testing data
following the wideband signal generation model described in Section 2.1. We
then devise a lightweight neural network architecture and train a deep
learning-based DoA classification model as described in Section 2.2. Finally,
we develop a GUI interface for deployment and performance validation as
discussed in Section 3. Experimental results are discussed in Section 4 and
Section 5 concludes the paper.

##  2\. Proposed Work

In this section, we discuss the wireless signal model used to generate the
dataset and, the design details of the neural network used to estimate the
DoA.

Figure 2. Wideband spectrum digitization approach for dataset generation.

###  2.1. Wireless Signal Model for Dataset

We consider a wideband spectrum consisting of several disjoint, uncorrelated,
and far-field signals. These signals are imposed on the $L$-antenna receiver
with either Uniform Linear Array (ULA). The signal received at the
$l^{\mathrm{th}}$ antenna is defined as follows:

(1) $\displaystyle x_{l}(t)$
$\displaystyle=\sum_{m=1}^{M}a_{m}\left(t-\tau_{l}\left(\theta_{m}\right)\right)e^{j2\pi
f_{m}\left(t-\tau_{l}\left(\theta_{m}\right)\right)}+n_{l}(t)$ (2)
$\displaystyle\approx\sum_{m=1}^{M}a_{m}(t)e^{j2\pi
f_{m}\left(t-\tau_{1}\left(\theta_{m}\right)\right)}+n_{l}(t)$

where M is the number of narrow band signals, $a_{m}(t)$ is the amplitude of
$m^{\mathrm{th}}$ narrowband signal of carrier frequency, $f_{m}$ and DoA,
$\theta_{m}$. Furthermore, $\tau_{l}(\theta_{m})$ is the time delay observed
by the $m^{\mathrm{th}}$ signal at the $l^{\mathrm{th}}$ antenna and
$n_{l}(t)$ is the additive white Gaussian noise (AWGN) at the
$l^{\mathrm{th}}$ antenna based on the SNR of the communication channel. The
output of the $l^{\mathrm{th}}$ antenna is digitized by means of Nyquist
sampling where the received signal, $x_{l}(t)$ is mixed with a mixing
function, $m(t)$ = $\sum_{b\in\beta}\alpha_{l,b}e^{-j2\pi(b-1)Bt}$ as shown in
Fig. 2. Digitization is performed on the set of non-contiguous frequency bands
i.e. $\beta\in\\{1,2,\cdots,N\\}$. Here, N is the total number of non-
overlapping frequency bands in the wideband spectrum, $\alpha_{l,b}$ is a
mixing coefficient of $b^{th}$ frequency band at $l^{\text{th }}$ antenna
carefully chosen from the Gaussian distribution, and B is the bandwidth of a
frequency band. Next we define the Fourier transform (FT) of the mixed signal,
$\tilde{x}_{l}(t)$, as follows

(3) $\displaystyle\tilde{X}_{l}(f)=\sum_{m=1}^{M}e^{j2\pi
f_{m}\tau_{l}\left(\theta_{m}\right)}\sum_{b\in\beta}\alpha_{l,b}A_{m}\left(f-\left(f_{m}-(b-1)B\right)\right)$

where $A_{m}(f)$ is the FT of mth signal. The mixed signal is passed through
the low pass filter (LPF) of cut-off frequency B/2 as $\tilde{X}_{l}(f)$
covers images over the entire frequency. From the LPF, we get the filtered
signal which is digitized at a rate of B Hz. The output of the ADC can be
expressed as

(4) $\displaystyle{Y}_{l}(e^{j2\pi
fNT})=\sum_{b\in\beta_{busy}}\alpha_{l,b}A_{b}(f-(f_{b}-(b-1)B))e^{j2\pi
f_{b}\tau_{l}(\theta_{b})}$

where $\beta_{busy}\in\beta$ is a selected set of occupied frequency bands,
$A_{b}(f)$ is the FT of the signal in $b^{th}$ frequency band with frequency
and DoA of $f_{b}$ $Hz$ and $\theta_{b}$, respectively. The FT of the output
of all L antennas in the matrix form is given as

(5) $\displaystyle\textbf{Y}(f)=\textbf{S}\hskip 2.84526pt\textbf{Z}(f)$

where S is a $L\hskip 2.84526pt\times\hskip 2.84526ptM$ steering matrix where
$\textbf{S}_{l,m}=e^{j2\pi f_{m}\tau{l}(\theta_{m})}$. Note that the steering
matrix, S contains two variables, carrier frequency, $f_{m}$ and DoA,
$\theta_{m}$ for each of the narrowband signal. The signal received at the
antennas is added with AWGN noise based on the SNR value of the channel to
generate the final signal corresponding to the combination of DoA and SNR
values. In this work, we focus on DoA estimation with known carrier frequency
and hence, one frequency band is passed to neural network at a time. The
details of dataset generation are given in Section 4.

###  2.2. Neural Network Architecture

In this section, the proposed neural network architecture for DoA estimation
is presented. Specifically, we use a convolution-based neural network (CNN)
due to their capability in capturing spatial correlations within input signals
at different scales (or receptive field) as demonstrated in various computer
vision (Krizhevsky et al., 2012; Redmon et al., 2016; He et al., 2017, 2016;
Ren et al., 2015; Ronneberger et al., 2015; Shelhamer et al., 2017) and signal
processing (O’Shea et al., 2018, 2016; Chandhok et al., 2019; West and O’Shea,
2017) applications. They are also data-efficient due to the inherent inductive
bias of convolution operations (Cohen and Shashua, 2017). Furthermore, CNNs
use parameter sharing across inputs, employ operations that are easily both
parallelizable and pipelined and entail $O(1)$ sequential operations. This is
in contrast to recurrent networks (RNN) which employ sequential processing
leading to higher inference time or fully connected networks (FCN), which
require more parameters and thus more computational efforts and storage
resources when compared with the CNN counterparts (Zhang et al., 2021; Raja,
2021). Furthermore, CNNs are designed to seamlessly work with spatial signals
which are in contrast to FCN or RNNs which require additional processing steps
to tackle spatial input data. Thus, we choose CNN for building a lightweight
memory efficient, therefore deployable on MCU, framework.

The proposed CNN model consists of a total of 7 layers. The first two layers
are 2D convolutional layers with kernel size 2x2 and Rectified Linear Unit
(ReLU) (Agarap, 2018) activation function. Over multiple experiments, we
observed that using ReLU allows faster convergence, to learn and perform
better. The proposed architecture also includes two 2D max-pooling layers in
between which enable progressive reduction of the activations spatial size,
reduce the parameters and computational complexity of the network. Finally, we
have two Fully Connected Layers (FC). The first FC layer uses ReLU and the
last layer has a softmax function which maps the outputs to probability values
for classification. We train our model end-to-end with categorical cross-
entropy loss. For training it, we use the Adam optimizer (Kingma and Ba, 2015)
due to its training stability and fast convergence properties. Specifically,
we use a learning rate of $0.0001$ and $\beta_{1}=0.5$, $\beta_{2}=0.999$
parameters for Adam the optimizer. Our network is implemented using Keras
library (Chollet et al., 2015) with Tensorflow backend (Abadi et al., 2015).
Due to limited space constraints, we have not reported various CNN
configurations explored during this project. The discussion is limited to a
model which offered the best performance fitting also memory constraints of
the MCU.

The next step is to pre-process the model for deployment on the chosen micro-
controller and corresponding steps are given in Fig. 3. Depending on the
application requirements, the model is quantized with appropriate word length
(WL). Currently, three WLs, 64-bit double-precision floating-point (DPFL),
32-bit single-precision floating-point (SPFL), and 8-bit integer (INT8) are
supported. The quantized model is then exported to as a frozen .pb graph and
the file footprint is reduced by discarding metadata and gradient information
that were needed during the training procedure. The frozen graph is then
converted to the tflite file using the Tensorflow lite conversion procedure.
The resulting file is then used to deploy it on the STM32 MCU taking advantage
of the X-CUBE-AI automatic deployment procedure integrated into STM32CubeMX
(STMicroelectronics, [n.d.]). Please refer to video tutorials (demo link)
explaining the detailing function of the proposed DoA estimation framework.

Figure 3. Block diagram depicting the tensorflow steps used to convert
proposed CNN model in tflite format that can be deployed on the STM32 micro-
controller.

##  3\. GUI for Model Deployment on MCU

The deployment of the neural network model on a micro-controller and its
performance analysis needs an appropriate user-friendly Graphical User
Interface (GUI). Such GUI is application-dependent and needs to be developed
from scratch. In this section, we provide details of the proposed GUI which
enabled validation of the CNN model on the NUCLEO-H743ZI2 board and comparison
with x86 software implementation. Fig. 4 shows the proposed GUI in which users
can configure the system parameters such as SNR, number of antennas ($L$), DoA
resolution ($R$), and DoA of the active signal along with an appropriate model
for a given word-length (WL). The GUI then generates appropriate test data,
validates the desired model on the Nucleo-STM32 board, and displays the
results given by the software (x86) and hardware (STM32) implementations along
with the inference time.

Figure 4. Proposed GUI to validate the neural-network models in both x86
software and NUCLEO-STM32 hardware modes.

To perform the validation on Nucleo-STM32H743ZI2, we develop seamless
interaction between the GUI interface and Nucleo-STM32 board to communicate
model files and input/output data. In this work, we used X-CUBE-AI-based
python libraries to establish the communication between the PC and the NUCLEO-
STM32 MCU. There are two python modules composing the HOST API package which
play an important role in validating generated neural networks running on a
Nucleo-STM32 board connected to the host PC using a standard serial interface
(ST-Link V2/V3 virtual com port).

The first module is _stm32nn_ which allows discovering the neural networks
which are available in the Nucleo-STM32 firmware. It enables GUI application
to retrieve important characteristics of the deployed neural network (tensor
input/output format, multiply-accumulate operation complexity, RAM/ROM size
used by CNN, etc.) and perform inference with the client data (default mode).
Furthermore, it also helps to establish the UART communication with the
Nucleo-STM32 board. The second module is the _stm32pb_ which implements the
low-level input/output functions required to communicate with the firmware
running on the STM32 device. This firmware communicates with the GUI through
the UART interface. For doing Nucleo-STM32 validation we deployed the X-CUBE-
AI automatically generated aiValidaiton-application on the edge device
(Nucleo-STM32 board) with the help of STM32CubeIDE which is an Integrated
Development Environment for STM32 Microcontrollers. This is a built-in
application of X-CUBE-AI that doesn’t require any handcrafting to be deployed
on Nucleo-STM32.

We used the Nucleo-STM32H743ZI2 board to validate our DoA estimation
application. This board features a STM32H743ZI2 MCU (ARM Cortex-M7 core based)
which allows IoT practitioners to choose among different combinations of
performance and power consumption by selecting among four different CPU run-
modes namely Range-0 to Range-3. Run mode Range-0 provides the highest
performance at the cost of higher power consumption. Successive run-modes
(Range-1, Range-2, and Range-3) can be used to lower the system’s power
consumption and performance by lowering the CPU operating frequency.
STM32H743ZI2 MCU can run at clock frequencies of up to 480MHz and embeds 1
MByte RAM and 2 MBytes of non-volatile flash memory. Due to memory
constraints, we need to optimize the size of the CNN model which can be
embedded into it. X-CUBE-AI expansion (of STM32CubeMX) package was used to
automatically generate STM32-optimized C-code and run time library
corresponding to application-specific, pre-trained CNN model. Using X-CUBE-AI,
multiple neural network models can be embedded into the STM32H743ZI2
microcontroller provided that the combined size of all the models is within
the memory limit of the STM32H743ZI2. The GUI application also validates the
performance of the model on the host processor (x86), which is another feature
of X-CUBE-AI. It uses the TensorFlow Lite model file (.tflite) to validate the
accuracy of the neural network library generated by X-CUBE-AI. For GUI design,
the PyQt5 python library and Qt Designer platform is used for adding different
features to make GUI user-friendly.

##  4\. Experiments and Performance Analyses

In this section, we describe the experimental setup and characterize the
performances of the proposed framework in comparison with baseline approaches
for different values of $R$, $L$, SNR, and WL. We begin with the details of
the dataset generated in this project.

Proposed Dataset: The proposed dataset is generated by using the wideband
signal model (described in section 2.1) and MATLAB software. The data
generation code requires various input parameters such as $R$, SNR, $L$,
sampling frequency, and the number of observations $N$. For the experiments
presented in this section, we consider five different DoA resolutions, 5
different size antenna arrays, and 20 different SNR values as shown in Table
1. Thus, the total number of DoA classes ranges from 8 in $R=45$ to 40 in
$R=9$. We generate _N_ $=2500$ observations for each DoA class. Each
observation or data point comprises of complex-valued signal of dimension
$l\times d$, where $d$ denotes the dimension of each signal and $l$ is the
_antenna count_. The complex-valued data is then converted into the real and
imaginary parts. Unfortunately to date, Since Keras does not handle natively
complex numbers, they are then concatenated with each other to form a combined
matrix of size $(l*2)\times d$, which acts as the input to the neural network.
We keep $80\%$ samples as a training set and $20\%$ for the test set.
Furthermore, out of the training set, $20\%$ samples are kept aside as a
validation set. All datasets are free and available online (here).

Antenna Count ($L$) | 2-6
---|---
DoA Resolution ($R$) | 45, 36, 24, 18, 9
SNR | 0-20dB
DoA classes | $K=360/R$
Observations ($N$) | 2500 per class
Table 1. Summary of signal parameters used for dataset generation Figure 5.
Comparison of DoA accuracy in % for different values of $L$ and SNRs with a
fixed $R=45$.

First, we begin with analyzing the effect of $L$ on the DoA estimation
accuracy for different SNRs and $R$. In Fig. 5, we compare the DoA accuracy
for dataset with fixed DoA resolution of 45. Here, we vary the number of
antennas from 2-6 and the SNR range is 0-20 dB. It can be observed that the
DoA accuracy increases significantly with the increase in the number of
antennas due to additional information gained via a signal received from the
spatially separated antenna elements. Furthermore, DoA accuracy increases with
the increase in SNR due to the reduced impact of noise. Note that the increase
in DoA estimation accuracy with the increase in the number of antennas is
substantial at low SNR. Furthermore, for the model with only 2 antennas, the
average validation accuracy is $96\%$ and the average test accuracy is
$95.6\%$ over the chosen SNR range. In Fig. 6, we repeat the same experiments
for the dataset with the DoA resolution of 36 i.e. more number of candidate
DoAs. The results are similar to Fig. 5 except for slight degradation in
estimation accuracy. From Fig. 5 and Fig. 6, we can conclude that the higher
the number of antennas, higher is the accuracy. However, the gain in
performance is not significant for the higher number of antennas which means
the number of antennas must be selected carefully. This is because an increase
in the number of antennas results in an increase in the area and power
complexity of the analog front-end due to multiple radio-frequency chains and
higher computational complexity of digital-front-end and baseband algorithms
due to a higher number of digitized samples.

Figure 6. Comparison of DoA accuracy in % for different values of $L$ and SNRs
with a fixed $R=36$.

Next, we study the difference in the DoA estimation accuracy for different $R$
and subsequent improvement after increasing $L$. In Fig. 7, it can be observed
that the DoA estimation accuracy decreases significantly with the decrease in
the DoA resolution. This is due to the increase in the number of classes with
the decrease in $R$ making a classification problem more challenging. To
further increase the accuracy, we need to increase the number of antennas,
especially at low SNRs. This also demands a reconfigurable antenna and DoA
estimation architecture that can adapt based on the desired accuracy.
Similarly, we study the difference in the DoA estimation accuracy for
different $R$ and subsequent improvement after increasing SNR. Since SNR
depends on the wireless environment, it can not be controlled and hence, the
wireless channel, with the desired SNR, must be chosen via appropriate
resource allocation strategies.

Figure 7. Comparison of the effect of $R$ and $L$ on the DoA accuracy in % for
different SNR. Figure 8. Comparison of the effect of $R$ and $SNR$ on the DoA
accuracy in % for different $L$.

Next, we study the effect of word length on the performance of the proposed
algorithm. As shown in Fig. 9(a) and Fig. 9(b), WL of 8-bit integer results in
slight degradation of the performance compared to SPFL WL. We also observed
that the SPFL and DPFL WLs offer nearly identical performance and the DPFL
results are excluded to maintain the clarity of the plots. The WL optimization
is critical to meet the memory constraints of the micro-controller. For
instance, the proposed algorithm with $R=45$, and $L=2$ needs 0.8 Mebibyte
(MiB) and 0.3 kibibytes (KiB) flash and on-chip RAM, respectively for SPFL WL.
This reduces to 0.2 MiB (-25%) and 0.09 KiB (-33%) of flash and on-chip RAM,
respectively for INT8 WL. Furthermore, the number of multiply-accumulate
computations is reduced from 431560 SPFL operations to 424568 INT8 operations.
This has resulted in the reduction in execution time from 6.4 ms to 2 ms
(-32%). Thus, careful selection of WL can offer significant savings in memory,
computational complexity and execution time without compromising on the DoA
estimation accuracy significantly.

(a)

(b)

Figure 9. Comparison of the effect of WL on the DoA accuracy in % for (a)
$R=45$ and (b) $R=36$ with $L=\\{2,4\\}$.

In the end, we present the comparison of the proposed algorithm with the well-
known spatial spectrum estimation approach known as the MUltiple SIgnal
Classification (MUSIC) algorithm (Gupta and Kar, 2015). The basic idea of the
MUSIC algorithm is to perform characteristic decomposition of the covariance
matrix of an array output data, resulting in a signal subspace orthogonal with
a noise subspace corresponding to the signal components. The two orthogonal
subspaces are then used to constitute a spectrum function that can be used for
spectral peak search and detect DoA signals (Gupta and Kar, 2015). In Fig.
10(a), we compare the performance of proposed and MUSIC algorithms for
$R=\\{45,36\\}$ with $L=2$. It can be observed that the proposed algorithm
outperforms MUSIC at all SNRs and offers an average of 5% higher accuracy. As
shown in Fig. 10(b), the performance of both algorithms improves with the
increase in $L$ from 2 to 4. However, the gain in performance is significant
in the proposed algorithm compared to MUSIC due to carefully designed CNN
architecture exploiting the spatial information received from additional
antennas.

(a)

(b)

Figure 10. Comparison of DoA accuracy in % for MUSIC and proposed approach
for: (a) Different $R$, and (b) Different $L$ .

##  5\. Conclusion

In this work, we developed a lightweight, accurate deep learning-based
framework for the task of direction-of-arrival (DoA) estimation in wireless
networks. We developed a user-friendly graphical user interface (GUI) to
deploy and validate the performance of the proposed algorithm on any micro-
controller such as the STM32H743ZI2 development board. In-depth performance
analysis for different signal-to-noise ratios (SNR), word-length, number of
antennas, and DoA resolution validated the superiority of the proposed
algorithm over the state-of-the-art approach. Future works include the
extension of the proposed framework for multi-band DoA classification and
real-time tracking of DoA for upcoming vehicular networks. Furthermore, we
would like to incorporate the support for additional word lengths such as
INT16 and INT32 which are not available in existing tools.

## References

  * (1)
  * Abadi et al. (2015) Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng. 2015\.  TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems.  https://www.tensorflow.org/ Software available from tensorflow.org.
  * Agarap (2018) Abien Fred Agarap. 2018\.  Deep Learning using Rectified Linear Units (ReLU).  _ArXiv_ abs/1803.08375 (2018).
  * Ahmad et al. (2020) W. S. H. M. W. Ahmad, N. A. M. Radzi, F. S. Samidi, A. Ismail, F. Abdullah, M. Z. Jamaludin, and M. N. Zakaria. 2020\.  5G Technology: Towards Dynamic Spectrum Sharing Using Cognitive Radio Networks.  _IEEE Access_ 8 (2020), 14460–14488.  https://doi.org/10.1109/ACCESS.2020.2966271
  * Al-Garadi et al. (2020) Mohammed Ali Al-Garadi, Amr Mohamed, Abdulla Khalid Al-Ali, Xiaojiang Du, Ihsan Ali, and Mohsen Guizani. 2020.  A Survey of Machine and Deep Learning Methods for Internet of Things (IoT) Security.  _IEEE Communications Surveys Tutorials_ 22, 3 (2020), 1646–1685.  https://doi.org/10.1109/COMST.2020.2988293
  * Amin and Hossain (2021) Syed Umar Amin and M. Shamim Hossain. 2021.  Edge Intelligence and Internet of Things in Healthcare: A Survey.  _IEEE Access_ 9 (2021), 45–59.  https://doi.org/10.1109/ACCESS.2020.3045115
  * Cass (2019) Stephen Cass. 2019\.  Taking AI to the edge: Google’s TPU now comes in a maker-friendly package.  _IEEE Spectrum_ 56, 5 (2019), 16–17.  https://doi.org/10.1109/MSPEC.2019.8701189
  * Chandhok et al. (2019) Shivam Chandhok, H. Joshi, A. Subramanyam, and S. Darak. 2019\.  Novel Deep Learning Framework for Wideband Spectrum Characterization at Sub-Nyquist Rate.  _arXiv: Signal Processing_ (2019).
  * Chang et al. (2021) Zhuoqing Chang, Shubo Liu, Xingxing Xiong, Zhaohui Cai, and Guoqing Tu. 2021.  A Survey of Recent Advances in Edge-Computing-Powered Artificial Intelligence of Things.  _IEEE Internet of Things Journal_ (2021), 1–1.  https://doi.org/10.1109/JIOT.2021.3088875
  * Chettri and Bera (2020) Lalit Chettri and Rabindranath Bera. 2020.  A Comprehensive Survey on Internet of Things (IoT) Toward 5G Wireless Systems.  _IEEE Internet of Things Journal_ 7, 1 (2020), 16–32.  https://doi.org/10.1109/JIOT.2019.2948888
  * Chollet et al. (2015) François Chollet et al. 2015\.  Keras.  https://keras.io.
  * Cohen and Shashua (2017) Nadav Cohen and A. Shashua. 2017.  Inductive Bias of Deep Convolutional Networks through Pooling Geometry.  _ArXiv_ abs/1605.06743 (2017).
  * Dohler et al. (2017) Mischa Dohler, Toktam Mahmoodi, Maria A. Lema, Massimo Condoluci, Fragkiskos Sardis, Konstantinos Antonakoglou, and Hamid Aghvami. 2017.  Internet of skills, where robotics meets AI, 5G and the Tactile Internet. In _2017 European Conference on Networks and Communications (EuCNC)_. 1–5.  https://doi.org/10.1109/EuCNC.2017.7980645
  * Ferreira et al. (2018) Joao Ferreira, Gustavo Callou, Albert Josua, and Paulo Maciel. 2018.  Estimating the Environmental Impact of Data Centers. In _2018 IEEE 17th International Symposium on Network Computing and Applications (NCA)_. 1–4.  https://doi.org/10.1109/NCA.2018.8548326
  * Gao et al. (2016) Zhen Gao, Linglong Dai, Wei Dai, Byonghyo Shim, and Zhaocheng Wang. 2016.  Structured Compressive Sensing-Based Spatio-Temporal Joint Channel Estimation for FDD Massive MIMO.  _IEEE Transactions on Communications_ 64, 2 (2016), 601–617.  https://doi.org/10.1109/TCOMM.2015.2508809
  * Gupta and Kar (2015) Pooja Gupta and S.P. Kar. 2015.  MUSIC and improved MUSIC algorithm to estimate direction of arrival. In _2015 International Conference on Communications and Signal Processing (ICCSP)_. 0757–0761.  https://doi.org/10.1109/ICCSP.2015.7322593
  * He et al. (2017) Kaiming He, Georgia Gkioxari, Piotr Dollár, and Ross B. Girshick. 2017.  Mask R-CNN.  _2017 IEEE International Conference on Computer Vision (ICCV)_ (2017), 2980–2988.
  * He et al. (2016) Kaiming He, X. Zhang, Shaoqing Ren, and Jian Sun. 2016\.  Deep Residual Learning for Image Recognition.  _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ (2016), 770–778.
  * Hussain et al. (2018) Ahmed A. Hussain, Nizar Tayem, Muhammad Omair Butt, Abdel-Hamid Soliman, Abdulrahman Alhamed, and Saleh Alshebeili. 2018.  FPGA Hardware Implementation of DOA Estimation Algorithm Employing LU Decomposition.  _IEEE Access_ 6 (2018), 17666–17680.  https://doi.org/10.1109/ACCESS.2018.2820122
  * Joshi et al. (2021) Himani Joshi, Sumit J. Darak, Mohammad Alaee-Kerahroodi, and Bhavani Shankar Mysore Rama Rao. 2021\.  Reconfigurable and Intelligent Ultrawideband Angular Sensing: Prototype Design and Validation.  _IEEE Transactions on Instrumentation and Measurement_ 70 (2021), 1–15.  https://doi.org/10.1109/TIM.2021.3058406
  * Kingma and Ba (2015) Diederik P. Kingma and Jimmy Ba. 2015.  Adam: A Method for Stochastic Optimization.  _CoRR_ abs/1412.6980 (2015).
  * Krizhevsky et al. (2012) A. Krizhevsky, Ilya Sutskever, and Geoffrey E. Hinton. 2012\.  ImageNet classification with deep convolutional neural networks.  _Commun. ACM_ 60 (2012), 84 – 90.
  * Li et al. (2020) En Li, Liekang Zeng, Zhi Zhou, and Xu Chen. 2020\.  Edge AI: On-Demand Accelerating Deep Neural Network Inference via Edge Computing.  _IEEE Transactions on Wireless Communications_ 19, 1 (2020), 447–457.  https://doi.org/10.1109/TWC.2019.2946140
  * Li et al. (2018) He Li, Kaoru Ota, and Mianxiong Dong. 2018.  Learning IoT in Edge: Deep Learning for the Internet of Things with Edge Computing.  _IEEE Network_ 32, 1 (2018), 96–101.  https://doi.org/10.1109/MNET.2018.1700202
  * Lien et al. (2017) S. Y. Lien, S. L. Shieh, Y. Huang, B. Su, Y. L. Hsu, and H. Y. Wei. June 2017\.  5G New Radio: Waveform, Frame Structure, Multiple Access, and Initial Access.  _IEEE Communication Magazine_ 55 (June 2017), 64–71.
  * Lin et al. (2020) Sen Lin, Zhi Zhou, Zhaofeng Zhang, Xu Chen, and Junshan Zhang. 2020.  Edge Intelligence in the Making: Optimization, Deep Learning, and Applications.
  * Liu et al. (2018) Zhang-Meng Liu, Chenwei Zhang, and Philip S. Yu. 2018\.  Direction-of-Arrival Estimation Based on Deep Neural Networks With Robustness to Array Imperfections.  _IEEE Transactions on Antennas and Propagation_ 66, 12 (2018), 7315–7327.  https://doi.org/10.1109/TAP.2018.2874430
  * Mao et al. (2018) Qian Mao, Fei Hu, and Qi Hao. 2018.  Deep Learning for Intelligent Wireless Networks: A Comprehensive Survey.  _IEEE Communications Surveys Tutorials_ 20, 4 (2018), 2595–2621.  https://doi.org/10.1109/COMST.2018.2846401
  * Memon et al. (2020) Jamshed Memon, Maira Sami, Rizwan Ahmed Khan, and Mueen Uddin. 2020.  Handwritten Optical Character Recognition (OCR): A Comprehensive Systematic Literature Review (SLR).  _IEEE Access_ 8 (2020), 142642–142668.  https://doi.org/10.1109/ACCESS.2020.3012542
  * O’Shea et al. (2016) T. O’Shea, Johnathan Corgan, and T. Clancy. 2016.  Convolutional Radio Modulation Recognition Networks.  _ArXiv_ abs/1602.04105 (2016).
  * O’Shea et al. (2018) T. O’Shea, Tamoghna Roy, and T. Clancy. 2018.  Over-the-Air Deep Learning Based Radio Signal Classification.  _IEEE Journal of Selected Topics in Signal Processing_ 12 (2018), 168–179.
  * Raja (2021) Sudeep Raja. 2021\.  FNNs, RNNs, LSTM and BLSTM.  http://cse.iitkgp.ac.in/~psraja/FNNs%20,RNNs%20,LSTM%20and%20BLSTM.pdf
  * Redmon et al. (2016) Joseph Redmon, S. Divvala, Ross B. Girshick, and A. Farhadi. 2016.  You Only Look Once: Unified, Real-Time Object Detection.  _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ (2016), 779–788.
  * Ren et al. (2015) Shaoqing Ren, Kaiming He, Ross B. Girshick, and J. Sun. 2015\.  Faster R-CNN: Towards Real-Time Object Detection with Region Proposal Networks.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 39 (2015), 1137–1149.
  * Ronneberger et al. (2015) O. Ronneberger, P. Fischer, and T. Brox. 2015.  U-Net: Convolutional Networks for Biomedical Image Segmentation. In _MICCAI_.
  * Sharma et al. (2018) Shree Krishna Sharma, Tadilo Endeshaw Bogale, Long Bao Le, Symeon Chatzinotas, Xianbin Wang, and Björn Ottersten. 2018.  Dynamic Spectrum Sharing in 5G Wireless Networks With Full-Duplex Technology: Recent Advances and Research Challenges.  _IEEE Communications Surveys Tutorials_ 20, 1 (2018), 674–707.  https://doi.org/10.1109/COMST.2017.2773628
  * Shelhamer et al. (2017) Evan Shelhamer, Jonathan Long, and Trevor Darrell. 2017\.  Fully Convolutional Networks for Semantic Segmentation.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 39 (2017), 640–651.
  * Shi et al. (2016) Weisong Shi, Jie Cao, Quan Zhang, Youhuizi Li, and Lanyu Xu. 2016.  Edge Computing: Vision and Challenges.  _IEEE Internet of Things Journal_ 3, 5 (2016), 637–646.  https://doi.org/10.1109/JIOT.2016.2579198
  * STMicroelectronics ([n.d.]) STMicroelectronics. [n.d.].  STM32CubeMX.  https://www.st.com/en/development-tools/stm32cubemx.html
  * Verhelst and Moons (2017) Marian Verhelst and Bert Moons. 2017.  Embedded Deep Neural Network Processing: Algorithmic and Processor Techniques Bring Deep Learning to IoT and Edge Devices.  _IEEE Solid-State Circuits Magazine_ 9, 4 (2017), 55–65.  https://doi.org/10.1109/MSSC.2017.2745818
  * Wang (2017) DeLiang Wang. 2017\.  Deep Learning Reinvents the Hearing Aid.  _IEEE Spectr._ 54, 3 (March 2017), 32–37.  https://doi.org/10.1109/MSPEC.2017.7864754
  * Wang et al. (2017) Tianqi Wang, Chao-Kai Wen, Hanqing Wang, Feifei Gao, Tao Jiang, and Shi Jin. 2017\.  Deep learning for wireless physical layer: Opportunities and challenges.  _China Communications_ 14, 11 (2017), 92–111.  https://doi.org/10.1109/CC.2017.8233654
  * Wang et al. (2019) Xiaofei Wang, Yiwen Han, Chenyang Wang, Qiyang Zhao, Xu Chen, and Min Chen. 2019\.  In-Edge AI: Intelligentizing Mobile Edge Computing, Caching and Communication by Federated Learning.  _IEEE Network_ 33, 5 (2019), 156–165.  https://doi.org/10.1109/MNET.2019.1800286
  * Wang et al. (2020) Yu Wang, Jie Yang, Miao Liu, and Guan Gui. 2020\.  LightAMC: Lightweight Automatic Modulation Classification via Deep Learning and Compressive Sensing.  _IEEE Transactions on Vehicular Technology_ 69, 3 (2020), 3491–3495.  https://doi.org/10.1109/TVT.2020.2971001
  * West and O’Shea (2017) Nathan E. West and T. O’Shea. 2017.  Deep architectures for modulation recognition.  _2017 IEEE International Symposium on Dynamic Spectrum Access Networks (DySPAN)_ (2017), 1–6.
  * Zhang et al. (2021) Aston Zhang, Zachary C. Lipton, Mu Li, and Alexander J. Smola. 2021.  Dive into Deep Learning.  _arXiv preprint arXiv:2106.11342_ (2021).
  * Zolanvari et al. (2019) Maede Zolanvari, Marcio A. Teixeira, Lav Gupta, Khaled M. Khan, and Raj Jain. 2019.  Machine Learning-Based Network Vulnerability Analysis of Industrial Internet of Things.  _IEEE Internet of Things Journal_ 6, 4 (2019), 6822–6834.  https://doi.org/10.1109/JIOT.2019.2912022


# Physics and modeling of multi-domain FeFET with domain wall induced Negative
Capacitance

Nilesh Pandey  Graduate Student Member, IEEE  and Yogesh Singh Chauhan
Fellow, IEEE Nilesh Pandey, and Yogesh Singh Chauhan are with Department of
Electrical Engineering, Indian Institute of Technology - Kanpur, Kanpur,
208016, India. e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>work was
supported in part by the Swarna Jayanti Fellowship under Grant
DST/SJF/ETA02/2017-18 and in part by the Department of Science and Technology
through the FIST Scheme under Grant SR/FST/ETII-072/2016.

###### Abstract

In this paper, we present the dynamics and modeling of multi-domains in the
ferroelectric FET (FeFET). Due to the periodic texture of domains, the
electrostatics of the FeFET exhibit an oscillatory conduction band profile. To
capture such oscillations, we solve coupled 2-D Poisson’s equation with the
net ferroelectric energy density (gradient energy + free energy +
depolarization energy) equation. Multi-domain dynamics are captured by
minimizing the net ferroelectric energy leading to a thermodynamically stable
state. Furthermore, we show that the motion of domain walls originates local
bound charge density in the ferroelectric region, which induces the negative
capacitance (NC) effect. The strength of domain wall-induced NC is determined
by the gradient energy of the ferroelectric material. FeFET exhibits
variability in the drain current with domain period due to the inherent NC
effect. Additionally, the impact of domain wall transition
(soft$\rightleftharpoons$hard) on the device’s electrostatic/transport is also
analyzed. The model also accurately captures both nucleations of a new domain
and the motion of the domain wall. Furthermore, the model is thoroughly
validated against experimental results and phase-field simulations.

Multi-Domain, Domain dynamics, Domain wall induced negative capacitance,
Ferroelectric energy dynamics, Green’s function.

##  1 Introduction

Ferroelectricity discovery in the Hafnium oxide-based material led to the
enormous research interest in the CMOS compatible ferroelectric FET (FeFET),
which possesses the profound possibility to utilize it as an efficient non-
volatile memory (NVM) [1]-[9]. Ferroelectric layer in the FeFET exhibits
multi-domain texture, which leads to an energy minima (thermodynamically
stable) state [10]-[16]. Numerous studies based on the phase-field (numerical)
simulations have demonstrated that the density of domains in the FE region is
determined by the interplay between various ferroelectric energy
components[17]-[23]. Furthermore, the presence of negative capacitance (NC)
effect via domain wall motion is observed in the ferroelectric material
[24]-[26]. The concept of negative capacitance is extensively studied in FETs
to break the fundamental limit of sub-threshold slope [28]-[43]. Nevertheless,
the physics of domain formation in the ferroelectric material is crystalline.
But, the modeling (analytical/compact) and physics of domain dynamics in a
FeFET remains elusive.

It is well established that the interaction between various energy components
determines the state and density of domains in a ferroelectric material
[10]-[12],[17]-[26]. Therefore, we follow the same approach to capture both
nucleation and motion of the domain wall. In this paper, we present an
explicit analytical model of multi-domain FeFET. The 2-D Poisson’s equation is
solved with lateral and vertical gradients in the domain polarization.
Subsequently, obtained electric field profile is used to calculate the
ferroelectric layer’s gradient and depolarization energy density. Net
ferroelectric energy (free energy + gradient energy + depolarization energy)
is minimized to model the dynamics of multi-domains. Furthermore, the model
can capture the secondary effects such as negative capacitance induced by
domain wall motion and domain wall transition from soft $\rightleftharpoons$
hard type caused by ferroelectric thickness scaling. The negative capacitance
effect is a vital function of domain wall width, which causes variability in
the drain current characteristics. Additionally, the experimental results and
phase-field simulations thoroughly validate the developed model and algorithm.
Fig. 1(a) shows the schematic of an double-gate FeFET with the coordinate axis
used in this paper.

Figure 1: (a) Schematic of double gate FeFET (the half structure is shown).
(b) Validations of the developed model with the phase-field simulations [26]
and experimental data [27]. The metal work-function and ferroelectric material
parameters are tuned to calibrate the model. Default ferroelectric material
parameters are taken from [26].  Figure 2: (a) 3-D polarization distribution
in the FE region. (b) 2-D surface plot of the polarization profile.
Polarization is maximum at the metal/FE interface and minimum at FE/DE
interface. (c) Schematic of polarization wave model with an upward and
downward domain. (d) Characteristics of y-directional function f(y) used to
incorporate gradient in the polarization along the vertical direction. Default
parameters: L = 50 nm, $t_{fe}$ = 4 nm, $t_{ox}$ = 1 nm, W = 1 $\mu m$,
$V_{gs}$ = $V_{ds}$ = 0 V, $\epsilon_{x}$ = 22, $\epsilon_{y}$ = 18 and
$\epsilon_{ox}$ = 3.9, $\epsilon_{si}$ = 11.8, $t_{1}=t_{fe}+t_{ox}$, and
$t_{2}=t_{1}+t_{si}$.

##  2 Development of Electrostatics model with Multi-Domains

The electrostatic model of multi-domain FeFET is mainly divided into four
subsections. In the first part, the mathematical formulation of the
polarization profile is obtained. Subsequently, the polarization wave model is
used in the 2-D Poisson’s equation to derive the 2-D potential functions in
the various regions. Afterward, obtained electrostatic model (electric fields
and polarization profile) is used to calculate the net ferroelectric energy.
Finally, the domain period and dynamics of domains in the ferroelectric region
are captured by minimizing net ferroelectric energy.

###  2.1 Polarization wave model

We assume that polarization of multi-domain in the ferroelectric region is
periodic with upward and domain period of $x_{a}$ and $x_{b}$, respectively,
and transition from upward to downward domain occurs in the finite domain wall
width ($d_{w}$). The periodic texture of domains and finite domain wall width
are valid approximations observed in experiments and numerical phase-field
simulations [10]-[13],[17]-[20]. Therefore, the mathematical formulation of
the polarization profile is obtained in Fourier series form.

$\displaystyle P(x,y)=$ $\displaystyle
f(y)\left(\frac{a_{0}}{2}+\frac{2P_{s}}{d_{w}T_{D}}\sum_{j}\frac{A_{j}cos\left(k_{j}x\right)+B_{j}sin\left(k_{j}x\right)}{k_{j}^{2}}\right)$
(1)

where, $T_{D}$ is the total domain period ($x_{a}+x_{b}$), $d_{w}$ is the
domain wall width, $P_{s}$ is the spontaneous polarization and all other
Fourier series coefficients are given in the Appendix section.

The y-directional function $f(y)$ is included in (1) to incorporate the
vertical directional gradient in domain polarization.

$\displaystyle f(y)=\left(\frac{-1}{t_{fe}^{2}}\right)y^{2}+1$ (2)

Due to the gradient in polarization along the y-direction, we formulate the
function $f(y)$ such that: it is maximum at the FE/metal interface ($y=0$),
and it is minimum at the FE/DE interface ($y=t_{fe}$). Therefore, the
polarization of a domain is maximum at the metal/FE interface and approaches
zero at the FE/DE interface (also observed in phase-field simulations, [17],
[18]-[20]). Fig. 2 (a) and (b) show polarization’s 3-D and surface
distribution in the FE region. Fig. 2(c) shows the schematic of polarization
wave with a domain wall width $d_{w}$. Fig. 2(d) shows the characteristic of
$f(y)$.

###  2.2 Electrostatic Potential model

Fig. 1(a) shows the schematic diagram of the double gate multi-domain FeFET.
The 2-D Poisson’s equation in the various regions is expressed as.

1\. Ferroelectric region ($0<y<t_{fe}$)

$\displaystyle\epsilon_{x}\frac{\partial^{2}\phi_{fe}}{\partial
x^{2}}+\epsilon_{y}\frac{\partial^{2}\phi_{fe}}{\partial
y^{2}}=-\frac{1}{\epsilon_{0}}\left(\frac{\partial P(x,y)}{\partial
x}+\frac{\partial P(x,y)}{\partial y}\right)$ (3)

2\. Insulator region ($t_{fe}<y<t_{fe}+t_{ox}$)

$\displaystyle\frac{\partial^{2}\phi_{ox}}{\partial
x^{2}}+\frac{\partial^{2}\phi_{ox}}{\partial y^{2}}=0$ (4)

3\. Channel region ($t_{fe}+t_{ox}<y<t_{fe}+t_{ox}+t_{si}$)

$\displaystyle\frac{\partial^{2}\phi^{si}}{\partial
x^{2}}+\frac{\partial^{2}\phi^{si}}{\partial
y^{2}}=\frac{qN_{a}}{\epsilon_{si}}$ (5)

This work aims to model the conduction band barrier’s oscillations in a FeFET,
which are only significant in the weak inversion region. Therefore, we assume
the depletion approximation in the channel $\left(\ref{poisson_si}\right)$.
The elementary task is to capture domain nucleation and dynamics in the FE
region. Domain dynamics incorporate vertical and lateral gradients in
polarization waves (3). The y-directional gradient induces the bound charges
($\rho=-{\partial P}/{\partial y}$) at the FE-DE interface. On the other hand,
the lateral gradient in the polarization contributes to domain wall energy
density. We will show that various ferroelectric energy components interplay
determines a thermodynamically stable domain state in the FE region.

Eq. (3)-(5) are a system of 2-D non-homogeneous partial differential equation
(PDE). Green’s function approach obtains the analytical solutions of these
equations. Earlier, Green’s function method is used to obtain the 2-D
potential distributions of the various DG-MOSFET structures [44]-[51]. In our
recent work, an analytical model of MFIS-NCFET is reported[51]. However, the
only mono-domain state is considered in the FE region. Hence, the PDE in FE
region reduces to a Laplace equation $\left({\partial P}/{\partial
x}={\partial P}/{\partial y}=0\right)$. Therefore, the developed model in [51]
cannot capture and study the domain dynamics in the ferroelectric region.

Here, for the first time, we report an analytical and explicit model for
multi-domain FeFET with x-y gradients in the polarization (3).

$\displaystyle\phi^{fe}(x,y)=\phi_{G_{x}}(x,y)+\phi_{G_{y}}(x,y)+\phi^{fe}_{\rightleftharpoons}(x,y)+\phi^{fe}_{\pm}(x,y)$
(6)

Potential components $\phi_{G_{x}}(x,y)$, and $\phi_{G_{y}}(x,y)$ are obtained
by following Green’s identity.

$\displaystyle G_{i}=\iint\left(\frac{\partial P(x,y)}{\epsilon_{x}\partial
x}+\frac{\partial P(x,y)}{\epsilon_{y}\partial
y}\right){}G(x,y;x^{\prime},y^{\prime})dx^{\prime}dy^{\prime}$ (7)

$\displaystyle\phi_{G_{x}}(x,y)=\left(\frac{4P_{s}}{t_{fe}\epsilon_{y}d_{w}T_{D}}\right)\sum_{n}\frac{sin\left(k_{n}^{I}y\right)}{sinh\left(k_{n}^{I}L\right)}\left\\{\lambda^{n}\sum_{j}\left(\frac{B_{j}\xi^{j,n}(x)-A_{j}\zeta^{j,n}(x)}{k_{j}}\right)+\Lambda^{n}\sum_{j}\left(\frac{A_{j}\xi^{j,n}(x)+B_{j}\zeta^{j,n}(x)}{k_{j}^{2}}\right)\right\\}$
$\displaystyle+\left(\frac{a_{0}}{t_{fe}\epsilon_{y}}\right)\sum_{n}\frac{sin\left(k_{n}^{I}y\right)\Lambda^{n}}{\left(k_{n}^{I}\right)^{2}}\left\\{1-\left(\frac{sinh\left(k_{n}^{I}x\right)+sinh\left(k_{n}^{I}\left(L-x\right)\right)}{sinh\left(k_{n}^{I}L\right)}\right)\right\\}$
(8)
$\displaystyle\phi_{G_{y}}(x,y)=\left(\frac{4P_{s}}{L\epsilon_{x}d_{w}T_{D}}\right)\sum_{m}\frac{sin\left(k_{m}x\right)}{k_{m}}\left\\{\xi^{m}(y)\sum_{j}\left(\frac{B_{j}\lambda^{m}-A_{j}\Lambda^{m}}{k_{j}}\right)+\zeta^{m}(y)\sum_{j}\left(\frac{A_{j}\lambda^{m}+B_{j}\Lambda^{m}}{k_{j}^{2}}\right)\right\\}$
(9) $\zeta^{j,n}(x)$, and $\xi^{j,n}(x)$ are function of $x$ with two series
indices $n$, and $j$.
$\displaystyle\xi^{j,n}(x)=\left(\frac{1}{k_{j}^{2}+\left(k_{n}^{I}\right)^{2}}\right)\left\\{cos\left(k_{j}x\right)sinh\left(k_{n}^{I}L\right)-sinh\left(k_{n}^{I}\left(L-x\right)\right)-sinh\left(k_{n}^{I}x\right)cos\left(k_{j}L\right)\right\\}$
(10)
$\displaystyle\zeta^{j,n}(x)=\left(\frac{1}{k_{j}^{2}+\left(k_{n}^{I}\right)^{2}}\right)\left(sin\left(k_{j}x\right)sinh\left(k_{n}^{I}L\right)-sin\left(k_{j}L\right)sinh\left(k_{n}^{I}x\right)\right)$
(11) $\zeta^{m}(y)$, and $\xi^{m}(y)$ are function of $y$ with $m$ series
indices
$\displaystyle\xi^{m}(y)=\left\\{\frac{-1}{t_{fe}^{2}}\left[\frac{y^{2}}{k_{m}}+\frac{2}{k_{m}^{3}}\left(1-\frac{cosh\left(k_{m}\left(t_{fe}-y\right)\right)}{cosh\left(k_{m}t_{fe}\right)}\right)-\frac{2t_{fe}}{k_{m}^{2}}\left(\frac{sinh\left(k_{m}y\right)}{cosh\left(k_{m}t_{fe}\right)}\right)\right]+\frac{1}{k_{m}}\left(1-\frac{cosh\left(k_{m}\left(t_{fe}-y\right)\right)}{cosh\left(k_{m}t_{fe}\right)}\right)\right\\}$
(12)
$\displaystyle\zeta^{m}(y)=\left\\{\frac{-1}{t_{fe}^{2}}\left[\frac{y}{k_{m}}-\frac{sinh\left(k_{m}y\right)}{k_{m}^{2}cosh\left(k_{m}t_{fe}\right)}\right]\right\\}$
(13)

Gradients of polarization wave ${\partial P(x,y)}{/\partial x}$, and
${\partial P(x,y)}{/\partial y}$ are evaluated by (1). Subsequently, evaluated
gradients and Green’s functions of the FE region are plugged in (7), which
leads to the 2-D potential distributions given by $\left(\ref{phi_Gx}\right)$,
and $\left(\ref{phi_Gy}\right)$, respectively. The remaining two terms
$\phi^{fe}_{\rightleftharpoons}(x,y),\;\text{and}\;\phi^{fe}_{\pm}(x,y)$
originated due to the FE region’s left/right ($x=0,\;\text{and}\;x=L$) and
top/bottom ($y=0,\;\text{and}\;y=t_{fe}$) boundaries, respectively (given in
the Appendix section).

The 2-D potential equations of oxide and channel regions are obtained by using
Green’s functions of respective regions into the Green’s identity[44].

$\displaystyle\phi^{ox}(x,y)=\phi_{\rightleftharpoons}^{ox}(x,y)+\phi_{\pm}^{ox}(x,y)+\phi^{ox}_{n=0}(x)$
(14)
$\displaystyle\phi^{si}(x,y)=E_{g}/2q+V_{ds}\frac{x}{L}+\frac{2}{L}\sum_{m=1}^{\infty}sin\left(k_{m}x\right)D_{sb}^{m,n}\times$
$\displaystyle\left(\frac{cosh\left(k_{m}\left(t_{2}-y\right)\right)+cosh\left(k_{m}\left(t_{1}-y\right)\right)}{\epsilon_{si}k_{m}sinh(k_{m}t_{si})}\right)$
(15)

The mathematical expressions of $\phi^{ox}(x,y)$ and approach to calculating
lateral Fourier series $\left(D_{sf}^{m,n}\text{, and }D_{sb}^{m,n}\right)$
coefficients are given in the Appendix.

###  2.3 Energy dynamics in the Ferroelectric material

The equilibrium (applied bias = 0 V) and non-equilibrium (applied bias $\neq$
0 V) configuration of the domains are determined by minimizing net
ferroelectric energy leading to a thermodynamically stable state. The net
energy density of the ferroelectric material is given as [11]-[13].

$\displaystyle f_{net}=f_{free}+f_{dep}+f_{grad}$ (16)

where, $f_{free}$ is free energy density and expressed in the Taylor series of
ferroelectric polarization [11]-[12].

$\displaystyle f_{free}=f_{0}+\alpha P^{2}(x,y)+\beta P^{4}(x,y)+\gamma
P^{6}(x,y)$ (17)

Depolarizing energy density is evaluated by the local distributions of
polarization texture and local electric fields [17].

$\displaystyle f_{dep}=$
$\displaystyle\frac{1}{2}\left[\epsilon_{x}\left(E^{fe}_{x}(x,y)\right)^{2}+\epsilon_{y}\left(E^{fe}_{y}(x,y)\right)^{2}\right]+E^{fe}_{y}(x,y)P(x,y)$
(18)

Contributions of local stray fields and polarization in $f_{dep}$ are
incorporated by the first and second terms of (18), respectively.

$\displaystyle E^{fe}_{x}(x,y)=-\frac{\partial\phi_{fe}(x,y)}{\partial x}$
(19) $\displaystyle E^{fe}_{y}(x,y)=-\frac{\partial\phi_{fe}(x,y)}{\partial
y}$ (20)

Domain polarization of ferroelectric material exhibits gradients along the
spatial $x$ and $y$ directions [17]-[20]. These gradients in the polarization
are incorporated by the ${\partial P(x,y)}/{\partial x}$ and ${\partial
P(x,y)}/{\partial y}$ in (3). The gradient energy density of the ferroelectric
material is expressed as.

$\displaystyle f_{grad}=\frac{1}{2}\left[\frac{g_{11}}{2}\left(\frac{\partial
P(x,y)}{\partial x}\right)^{2}+\frac{g_{44}}{2}\left(\frac{\partial
P(x,y)}{\partial y}\right)^{2}\right]$ (21)

where, the material-dependent ferroelectric coefficients are taken from [26].
We consider a stress-free condition in the Landau coefficients, and strain is
incorporated in the Taylor series coefficients of (17). This approximation is
validated by the numerical phase-field simulations [18].

###  2.4 Algorithm to calculate domain period for zero and non-zero applied
bias

Firstly, we evaluate the equilibrium domain period ($x_{e}$) for the zero
applied voltages. At zero bias, upward and downward domain, periods will be
identical [10]-[13], [17]-[20]. Therefore, the total domain period is
calculated as.

$\displaystyle T_{D}=x_{a}+x_{b}+4d_{w}\;(x_{a}=x_{b}=x_{e})\Rightarrow
T_{D}=2x_{e}+4d_{w}$ (22)

The equilibrium domain period is evaluated by minimizing net ferroelectric
energy $\left(F_{net}\right)$.

$\displaystyle
F_{net}=W\int_{0}^{t_{fe}}\int_{0}^{L}\left(f_{free}+f_{dep}+f_{grad}\right)dxdy$
(23)

where, $W$ is the width of device. Typical value of domain wall width is
$\approx$ 0.5 nm$-$1 nm [17]-[20]. Hence, after considering an appropriate
value of $d_{w}$ the only unknown in (23) is $x_{e}$, which is calculated by
the minimization of net ferroelectric energy as.

$\displaystyle\frac{\partial F_{net}}{\partial x_{e}}=0\;(\text{solve for
$x_{e}$})$ (24)

For a given set of physical parameters such as: $t_{fe}$, $t_{ox}$, $L$,
$t_{si}$, and ferroelectric material parameters ($\alpha$, $\beta$, $\gamma$
etc.) the above equation is an only function of $x_{e}$ which is solved
analytically.

The final task is to evaluate the domain period for non-zero applied voltages
(non-equilibrium state). Non-zero applied bias triggers the domain wall motion
[17]-[20]. Hence, the period of the upward and downward domain is no longer
remains identical ($x_{a}\neq x_{b}$). Depending on the polarity of applied
voltage, the expansion/reduction in a domain period is observed. An $\delta x$
variation in an upward domain period leads to the same amount of variation in
the downward domain period [17]-[20]. Therefore, the total domain period
remains the same as in an equilibrium condition period (22).

$\displaystyle T_{D}=x_{e}+\delta x+x_{e}-\delta x+4d_{w}=2x_{e}+4d_{w}$ (25)

However, due to the non-equilibrium condition, $x_{a}\neq x_{b}$, and the
amount of shift is given by $\delta x$, which is calculated by the
minimization of net ferroelectric energy.

$\displaystyle\frac{\partial F_{net}}{\partial\left(\delta
x\right)}=0\;(\text{solve for $\delta x$})$ (26)

Domain dynamics is captured by plugging $\delta x$ (For each applied voltages)
in (6).

##  3 Results and Discussion

To check the robustness and accuracy of the model, we consider both types of
validations: numerical simulations and experimental results. Fig. 1(b) shows
validations of the developed model against phase-field simulations [26], and
experimental results [27]. Fig. 3 shows the electrostatic potential
distribution (at zero bias) of a FeFET with multi-domains. Due to the adjacent
periodic nature of polarization charges in the FE region (see 2(a)), the
electrostatic of the device also exhibits periodicity. The value of
$\Delta\phi$ (difference between maximum and minimum values of the potential)
is $\approx$ 500 mV and 40 mV at the metal/FE ($y=0$) interface and at the
mid-channel $\left(y=t_{fe}+t_{ox}+t_{si}/2\right)$, respectively. Such
significant variations in the channel potential will exponentially alter the
thermionic current transport.

Figure 3: (a) 3-D surface plot of potential band energy. (b) 2-D surface plot
of electrostatic potential energy. Source side Fermi level is taken as
reference. Parameters: if not mentioned, then default values of parameters are
taken.  Figure 4: (a) Bound charge density in the FE region due to multi-
domains. (b) Emergence of local negative permittivity caused by domain wall
motion.

###  3.1 Negative Capacitance via Domain wall motion

The lateral gradient of polarization $\left(\partial P/\partial x\right)$
profile in (3) induces the bound charge density along the lateral direction.
Fig. 4(a) shows the bound charge density ($Q_{b}$) distribution in the FE
region for various domain periods. Charge density $Q_{b}$ exhibits the
periodic nature with positive and negative values. Fig. 4(b) shows the local
permittivity of the FE region (for a domain period = $T_{D}$), which is
calculated by the electric displacement vector. The emergence of negative
permittivity is observed due to the negative slope of polarization vector
w.r.t. electric field vector $\left(\partial P/\partial E<0\right)$. The
signature of negative permittivity is the origin of negative capacitance
induced by domain wall motion. Note that, recently, an experimental study also
demonstrated the presence of NC by domain dynamics in a ferroelectric
capacitor [25]. Therefore, the developed model can capture inherent negative
capacitance by domain wall motion.

Figure 5: (a) Conduction band energy plotted at mid-channel. Oscillations in
the conduction band are originated due to the multi-domains effect. (b) The
transition from mono-domain state to multi-domain state with the variations in
physical parameters.

Fig. 5(a) shows the conduction band profile plotted at the middle of the
channel for various $t_{fe}$ values. Due to enhanced fringing fields, the
negative capacitance effect shoots up as ferroelectric thickness increases
[33]. Therefore, the barrier height in the conduction band rises with an
increasing $t_{fe}$. The number of domains in the FE region is determined by
the interactions in various energy components given by (16). An enhancement in
$t_{fe}$ corresponds to a wider dispersion in domain polarization along the
vertical direction. Therefore, the slope of $f(y)$ decreases (see Fig. 2(d)),
which reduces gradient energy density along the y-direction. This reduction in
the gradient energy relaxes the compulsion to nucleate denser domain patterns
for achieving energy minima. Now, energy can be minimized by the lesser number
of domains. Hence, the domain period increases, which reduces the rate of
oscillations in the conduction band, as shown in Fig. 5(a). On the other hand,
as $t_{ox}$ increases, the depolarizing fields increases [10],[17]-[20], [24].
Therefore, more number domains nucleates to compensate for the increased
depolarizing energy, as shown in Fig. 5(b). Note that similar kinds of
observations were also reported in a pioneer work by Bratkovsky et al. [10],
and the phase-field simulations [17]-[20]. Therefore, Fig. 5(b) shows a
further validation of the developed model and algorithm to capture nucleation
of domains in the FE region.

###  3.2 Electronic transport with multi-domains

The drain to source current is calculated by the Quasi-Fermi potential of the
channel [51].

$\displaystyle I_{ds}=\frac{\mu
WkT\left(1-e^{-qV_{ds}/kT}\right)}{\int_{0}^{L}\frac{dx}{\int_{t_{1}}^{t_{1}+t_{si}}n_{i}e^{q\phi_{si}(x,y)/kT}dy}}$
(27)

Fig. 6(a) shows a short channel device’s conduction band energy plot (L = 15
nm). The number of domain walls varies with the variations in $t_{fe}$ (see
Fig. 5). Hence, the position of local NC shifts with the alteration in
$t_{fe}$. Thus, the maximum barrier height of the conduction band exhibits a
positional variation with ferroelectric thickness. Shifting in the barrier
height signifies the presence of a domain wall at that particular position.
Additionally, enhancement in the barrier height with a larger $t_{fe}$ is due
to the increased NC effect. Fig. 6(b) shows the $I_{ds}-V_{gs}$
characteristics of a FeFET. Solid curves are used for FeFET with multi-domains
which includes the negative capacitance effect (due to the domain wall
motion), and dashed curves are plotted for the conventional high-k FET where
the domain wall does not exist (hence, the NC effect is absent). Enhancement
in the $t_{fe}$ raises the conduction band barrier height (see Fig. 5(a) and
Fig. 6(a)), which leads to a reduction in OFF current (NC effect gets stronger
with the larger value of $t_{fe}$). On the other hand, conventional high-k FET
exhibits the reverse trend with $t_{fe}$, as shown in Fig. 5(a) (dashed
curves). The barrier height reduces with an increment in the FE layer
thickness. Therefore, the OFF current value rises with an increasing value of
$t_{fe}$.

Figure 6: (a) Mid-channel conduction band plot for the short channel device. A
shift in the barrier height occurs due to the nucleation of a domain wall. (b)
Drain current characteristics for short and long channel device. Figure 7: (a)
Transition of domain wall from hard to soft type with increasing value of
domain wall width. (b) The magnitude of local permittivity increases with the
larger $d_{w}$.

###  3.3 Soft $\rightleftharpoons$ Hard domain wall transition

The phase-field simulations demonstrated that scaling in the FE layer triggers
a domain wall transition from hard to soft type [18]-[20], [26]. This second
order effect is incorporated in the model by transforming the polarization
profile in the FE layer. Fig. 7(a) shows the polarization profile with the
variations in domain wall width ($d_{w}$). The hard domain wall texture is
observed at $d_{w}\approx 0$. However, a gradual enhancement in $d_{w}$ leads
to a soft domain wall transition. Therefore, transition of the domain wall
with $t_{fe}$ can be incorporated by defining a parameter $\alpha_{0}$.

$\displaystyle d_{w}=\frac{\alpha_{0}}{t_{fe}}d_{w_{0}}$ (28)

$\alpha_{0}$ can be tuned to capture soft/hard domain wall transitions, and
$d_{w_{0}}$ is the starting domain period (without transition). Impact of
domain wall transition on negative permittivity is shown in Fig. 7(b). The
magnitude of local permittivity increases with the larger $d_{w}$, which can
be understood by analyzing the ${\partial P}/{\partial E}$ ratio.

$\displaystyle\frac{\partial P}{\partial E}=\frac{{\partial P}/{\partial
x}}{{\partial E}/{\partial x}}$ (29)

An increase in $d_{w}$ reduces ${\partial P}/{\partial x}$ ratio (see Fig.
7(a)), however, the ratio ${\partial E}/{\partial x}$ also decreases with the
larger $d_{w}$. A larger permittivity value with the $d_{w}$ states that the
electric field vector is more sensitive than the polarization vector for
spatial variations $\left(\Delta\left({\partial E}/{\partial
x}\right)>\Delta\left({\partial P}/{\partial x}\right)\right)$. An alternative
explanation for increased permittivity with $d_{w}$ can be understood by
observing the conduction band profile, as shown in Fig. 8(a). As $d_{w}$
increases, the barrier height decreases. Thus, the NC’s strength reduces with
the larger $d_{w}$. The strength of negative capacitance is inversely
proportional to the ferroelectric layer capacitance $\left(\text{NC}\propto
1/\left|C_{FE}\right|\right)$, a larger value of $\left|C_{FE}\right|$
deteriorates the capacitance matching leading to a weaker NC effect
[Khan]-[36]. Since a larger magnitude of local permittivity enhances the
$\left|C_{FE}\right|$ value, the barrier height decreases with the increasing
$d_{w}$. Therefore, the transition towards the soft domain wall increases the
sub-threshold current, as observed in Fig. 8(b). Note that variability is
higher in the short channel device (L = 15 nm) than the long channel device (L
= 50 nm). This may be attributed to the fact that the SCEs dominate at smaller
channel lengths, which further reduces the impact of the NC effect, leading to
a higher variability.

.

Figure 8: (a) Mid-channel conduction band energy for the various value of
domain wall widths. Enhancement in $d_{w}$ leads to a weaker NC effect. (b)
Variability in the drain current caused by the domain wall width.

##  4 Conclusion

An analytical physics-based model of multi-domain FeFET is developed by
solving 2-D Poisson’s equation with the ferroelectric energy dynamics state
equation. The process of domain wall nucleation and motion occurs to minimize
the net system energy. The dynamics of multi-domains (nucleation of a new
domain and domain wall motion) are captured by minimizing net system energy.
Ferroelectric energy directly depends on physical parameters such as $t_{fe}$
and $t_{ox}$. Hence, the domain period can be engineered by the alterations in
the FE physical parameters. The motion of domain wall leads to the origination
of local negative permittivity in the FE layer. A wider domain wall enhances
the magnitude of local permittivity leading to a weaker negative capacitance
effect. Furthermore, the transition of domain wall
soft$\rightleftharpoons$hard introduces variability in the $I_{ds}-V_{gs}$
characteristics.

## References

  * [1] T. S. Boscke, J. Muller, D. Brauhaus, U. Schroder, and U. Bottger, “Ferroelectricity in hafnium oxide: CMOS compatible ferroelectric field effect transistors,” in IEDM Tech. Dig., Dec. 2011, pp. 24.5.1–24.5.4, doi: 10.1109/IEDM.2011.6131606.
  * [2] S. Mueller, J. Mueller, A. Singh, S. Riedel, J. Sundqvist, U. Schroeder, and T. Mikolajick, “ Incipient ferroelectricity in Al-doped HfO2 thin films,” Adv. Funct. Mater. 22, 2412–2417 (2012). https://doi.org/10.1002/adfm.201103119.
  * [3] J. Mu¨ller, T. S. Böscke, U. Schröder, S. Mueller, D. Brauhaus, U. Böttger, L. Frey, and T. Mikolajick, “ Ferroelectricity in simple binary ZrO2 and HfO2,” Nano Lett. 12, 4318 (2012). https://doi.org/10.1021/nl302049k.
  * [4] S. Fujii et al., “First demonstration and performance improvement of ferroelectric HfO2-based resistive switch with low operation current and intrinsic diode property,” in Proc. VLSI Tech. Symp., 2016, pp. 1–2.
  * [5] M. Trentzsch et al. A 28 nm HKMG super low power embedded NVM technology based on ferroelectric FETs, in IEEE International Electron Devices Meeting (2016), pp. 294–297.
  * [6] J. Muller, T. S. Boscke, U. Schroder, R. Hoffmann, T. Mikolajick, and L. Frey, “Nanosecond polarization switching and long retention in a novel MFIS-FET based on ferroelectric HfO2,” IEEE Electron Device Lett., vol. 33, no. 2, pp. 185–187, Feb. 2012, doi: 10.1109/ LED.2011.2177435.
  * [7] E. T. Breyer, H. Mulaosmanovic, T. Mikolajick, and S. Slesazeck, “Reconfigurable NAND/NOR logic gates in 28 nm HKMG and 22 nm FD-SOI FeFET technology,” in IEDM Tech. Dig., Dec. 2017, pp. 28.5.1–28.5.4, doi: 10.1109/IEDM.2017.8268471.
  * [8] H. Mulaosmanovic et al., “Evidence of single domain switching in hafnium oxide based FeFETs: Enabler for multi-level FeFET memory cells,” in IEDM Tech. Dig., vol. 3, Dec. 2015, pp. 26.8.1–26.8.3, doi: 10.1109/IEDM.2015.7409777.
  * [9] K. Florent et al., “Vertical ferroelectric HfO2 FET based on 3-D NAND architecture: Towards dense low-power memory,” in IEDM Tech. Dig., Dec. 2018, pp. 2–5, doi: 10.1109/IEDM.2018.8614710.
  * [10] A. M. Bratkovsky, & A. P. Levanyuk, Abrupt appearance of the domain pattern and fatigue of thin ferroelectric films. Phys. Rev. Lett. 84, 3177–3180, https://doi.org/10.1103/PhysRevLett.84.3177 (2000).
  * [11] M. E. Lines, & A. M. Glass, Principles and Applications of Ferroelectrics and Related Materials (Clarendon Press, Oxford, 1977).
  * [12] K. M. Rabe, C. H. Ahn & J. M. Triscone, Physics of Ferroelectrics: a Modern Perspective Vol. 105, Springer Science & Business Media (2007).
  * [13] I. Luk’yanchuk, , A. Sené, & V. M. Vinokur, Electrodynamics of ferroelectric films with negative capacitance. Phys. Rev. B 98, 024107 (2018).
  * [14] Anna N. Morozovska et al. Interaction of a 180° ferroelectric domain wall with a biased scanning probe microscopy tip: Effective wall geometry and thermodynamics in Ginzburg-Landau-Devonshire theory. Phys. Rev. B 78 (2008), 125407, https://doi.org/10.1103/PhysRevB.78.125407.
  * [15] E. A. Eliseev, A. N. Morozovska, S. V. Kalinin, Y. L. Li, Jie Shen, M. D. Glinchuk, L. Q. Chen, and V. Gopalan, Surface Effect on Domain Wall Width in Ferroelectrics, J. Appl. Phys. 106, 084102 (2009).
  * [16] I. A. Luk’yanchuk, L. Lahoche, & A. Sene, Universal Properties of Ferroelectric Domains. Phys. Rev. Lett. 102, 10.1103/PhysRevLett.102.147601 (2009).
  * [17] H. W. Park, J. Roh, Y. B. Lee & C. S. Hwang, Modeling of negative capacitance in ferroelectric thin films, Adv. Mater., vol. 31, Jun. 2019.
  * [18] A.K. Saha, S.K. Gupta, Multi-Domain Negative Capacitance Effects in Metal-Ferroelectric-Insulator-Semiconductor/Metal Stacks: A Phase-field Simulation Based Study. Sci Rep 10, 10207 (2020). https://doi.org/10.1038/s41598-020-66313-1.
  * [19] A. K. Saha, M. Si, K. Ni, S. Datta, P. D. Ye, and S. K. Gupta, “Ferroelectric thickness dependent domain interactions in FEFETs for memory and logic: A phase-field model based analysis,” in IEEE International Electron Devices Meeting (IEDM) (IEEE, San Francisco, CA, 2020), pp. 4.3.1–4.3.4.
  * [20] A. K. Saha & S. K. Gupta, Negative capacitance effects in ferroelectric heterostructures: A theoretical perspective”, J. Appl. Phys., vol. 129, no. 8, Feb. 2021.
  * [21] M. Hoffmann, et al., On the stabilization of ferroelectric negative capacitance in nanoscale devices. Nanoscale 10, 10891–10899 (2018).
  * [22] M. Hoffmann, et al., Unveiling the double-well energy landscape in a ferroelectric layer. Nature 565, 464–467 (2019). https://doi.org/10.1038/s41586-018-0854-z.
  * [23] P. Zubko, et al., Negative capacitance in multidomain ferroelectric superlattices. Nature 534, 524–528 (2016).
  * [24] A. M. Bratkovsky, & A. P. Levanyuk, Very large dielectric response of thin ferroelectric films with the dead layers. Phys. Rev. B 63, 132103, https://doi.org/10.1103/PhysRevB.63.132103 (2001).
  * [25] A. K. Yadav, et al., Spatially resolved steady-state negative capacitance. Nature 565, 468–471 (2019).
  * [26] A. K. Saha and S. K. Gupta, ”Multi-Domain Ferroelectric FETs with Negative and Enhanced Positive Capacitance for Logic Applications,” 2021 International Conference on Simulation of Semiconductor Processes and Devices (SISPAD), 2021, pp. 77-80, doi: 10.1109/SISPAD54002.2021.9592573.
  * [27] Z. Krivokapic et al., ”14nm Ferroelectric FinFET technology with steep subthreshold slope for ultra low power applications,” 2017 IEEE International Electron Devices Meeting (IEDM), 2017, pp. 15.1.1-15.1.4, doi: 10.1109/IEDM.2017.8268393.
  * [28] S. Salahuddin and S. Datta, “Use of negative capacitance to provide voltage amplification for low power nanoscale devices,” Nano Lett., vol. 8, no. 2, pp. 405–410, 2007, doi: 10.1021/nl071804g.
  * [29] G. A. Salvatore, D. Bouvet, and A. M. Ionescu, “Demonstration of subthrehold swing smaller than 60mV/decade in Fe-FET with P(VDF-TrFE)/SiO2 gate stack,” in IEDM Tech. Dig., Dec. 2008, pp. 1–4, doi: 10.1109/IEDM.2008.4796642.
  * [30] M. Kobayashi and T. Hiramoto, “Device design guideline for steep slope ferroelectric FET using negative capacitance in sub-0.2V operation: Operation speed, material requirement and energy efficiency,” in VLSI Symp. Tech. Dig., Jun. 2015, pp. T212–T213, doi: 10.1109/VLSIT.2015.7223678.
  * [31] H. Ota, T. Ikegami, J. Hattori, K. Fukuda, S. Migita, A. Toriumi, ”Fully coupled 3-D device simulation of negative capacitance FinFETs for sub 10 nm integration”, IEDM Tech. Dig., pp. 12.4.1-12.4.4, Dec. 2016, DOI: 10.1109/IEDM.2016.7838403
  * [32] W. Cao, , K Banerjee, Is negative capacitance FET a steepslope logic switch?. Nat Commun 11, 196 (2020). DOI: https://doi.org/10.1038/s41467-019-13797-9
  * [33] G. Pahwa, A. Agarwal, Y. S. Chauhan, ”Numerical investigation of short-channel effects in negative capacitance MFIS and MFMIS transistors: Subthreshold behavior”, IEEE Trans. Electron Devices, vol. 65, no. 11, pp. 5130-5136, Nov. 2018, DOI: 10.1109/TED.2018.2870519.
  * [34] G. Pahwa et al., ”Analysis and Compact Modeling of Negative Capacitance Transistor with High ON-Current and Negative Output Differential Resistance—Part II: Model Validation,” in IEEE Transactions on Electron Devices, vol. 63, no. 12, pp. 4986-4992, Dec. 2016, doi: 10.1109/TED.2016.2614436.
  * [35] G. Pahwa et al., ”Compact Model for Ferroelectric Negative Capacitance Transistor With MFIS Structure,” in IEEE Transactions on Electron Devices, vol. 64, no. 3, pp. 1366-1374, March 2017, doi: 10.1109/TED.2017.2654066.
  * [36] G. Pahwa et al., ”Physical Insights on Negative Capacitance Transistors in Nonhysteresis and Hysteresis Regimes: MFMIS Versus MFIS Structures,” in IEEE Transactions on Electron Devices, vol. 65, no. 3, pp. 867-873, March 2018, doi: 10.1109/TED.2018.2794499.
  * [37] X. Zhang et al., “Analysis on performance of ferroelectric NC-FETs based on real-space Gibbs-free energy with atomic channel structure,” IEEE Trans. Electron Devices, vol. 66, no. 2, pp. 1100-1106, 2019.
  * [38] Asif Islam Khan, et al. ”Negative capacitance in a ferroelectric capacitor.” Nature materials 14.2 (2015): 182.
  * [39] Ku, Hansol, and Changhwan Shin. ”Transient response of negative capacitance in P (VDF 0.75- TrFE 0.25) organic ferroelectric capacitor.” IEEE Journal of the Electron Devices Society 5.3 (2017): 232-236.
  * [40] M. Hoffmann, B. Max, T. Mittmann, U. Schroeder, S. Slesazeck and T. Mikolajick, ”Demonstration of High-speed Hysteresis-free Negative Capacitance in Ferroelectric Hf0.5Zr0.5O2,” 2018 IEEE International Electron Devices Meeting (IEDM), San Francisco, CA, 2018, pp. 31.6.1-31.6.4.
  * [41] Daniel JR Appleby, et al. ”Experimental observation of negative capacitance in ferroelectrics at room temperature.” Nano letters 14.7 (2014): 3864-3868.
  * [42] A. Cano and D. Jiménez, “Multidomain ferroelectricity as a limiting factor for voltage amplification in ferroelectric field-effect transistors,” Appl. Phys. Lett., vol. 97, no. 13, p. 133509, 2010.
  * [43] S. Kasamatsu, S. Watanabe, C. S. Hwang, and S. Han, “Emergence of negative capacitance in multidomain ferroelectric—Paraelectric nanocapacitors at finite bias,” Adv. Mater., vol. 28, no. 2, pp. 335–340, 2016.
  * [44] J. D. Jackson, Classical Electrodynamics3rd edn, chapter-1, (Wiley, (1999)).
  * [45] P.-S. Lin and C.-Y. Wu, “A new approach to analytically solving the two-dimensional Poisson’s equation and its application in short-channel MOSFET modeling,” IEEE Trans. Electron Devices, vol. ED-34, no. 9, pp. 1947–1956, Sep. 1987.
  * [46] J.-Y. Guo and C.-Y. Wu, “A new 2-D analytic threshold-voltage model for fully depleted short-channel SOI MOSFETs,” IEEE Trans. Electron Devices, vol. 40, no. 9, pp. 1653–1661, Sep. 1993, DOI: 10.1109/16.231571.
  * [47] A. Nandi, N. Pandey, S. Dasgupta, ”Analytical modeling of DG-MOSFET in subthreshold regime by green’s function approach”, IEEE Trans. Electron Devices, vol. 64, no. 8, pp. 3056-3062, Aug. 2017, 10.1109/TED.2017.2708603.
  * [48] A. Nandi and N. Pandey, “Accurate analytical modeling of junctionless DG-MOSFET by Green’s function approach,” Superlattices Microstruct., vol. 111, pp. 983–990, Nov. 2017, doi: 10.1016/j. spmi.2017.07.062.
  * [49] A. Nandi, N. Pandey and S. Dasgupta, ”Analytical Modeling of Gate-Stack DG-MOSFET in Subthreshold Regime by Green’s Function Approach,” in IEEE Transactions on Electron Devices, vol. 65, no. 10, pp. 4724-4728, Oct. 2018.
  * [50] N. Pandey, H.-H. Lin, A. Nandi, and Y. Taur, “Modeling of short-channel effects in DG MOSFETs: Green’s function method versus scale length model,” IEEE Trans. Electron Devices, vol. 65, no. 8, pp. 3112–3119, Aug. 2018, DOI: 10.1109/TED.2018.2845875.
  * [51] N. Pandey and Y. S. Chauhan, ”Analytical Modeling of Short-Channel Effects in MFIS Negative-Capacitance FET Including Quantum Confinement Effects,” in IEEE Transactions on Electron Devices, vol. 67, no. 11, pp. 4757-4764, Nov. 2020, doi: 10.1109/TED.2020.3022002.

## Appendix

$D_{sf}^{m,n}$, and $D_{sb}^{m,n}$ are calculated by the potential continuity
condition at FE/DE, and DE/channel interface respectively.

$\displaystyle
D_{sf}^{m,n}=\frac{d_{1}^{m}d_{\alpha}+d_{2}^{m}d_{\beta}}{d_{4}^{m}d_{1}^{m}-\left(d_{2}^{m}\right)^{2}};\;D_{sb}^{m,n}=\frac{d_{2}^{m}d_{\alpha}+d_{4}^{m}d_{\beta}}{d_{4}^{m}d_{1}^{m}-\left(d_{2}^{m}\right)^{2}}$
(30)

Note that, the dynamics of multi-domains introduces $d_{MD}^{m,n,j}$ Fourier
coefficient in the above equations.

$\displaystyle
d_{\alpha}=d_{5}^{m}+d_{MD}^{m,n,j}+d_{6}^{m,n}-d_{8}^{m,n}-d_{n=0}^{m};$ (31)
$\displaystyle d_{\beta}=d_{7}^{m,n}+d_{n=0}^{m}-d_{3}^{m}$ (32)
$\displaystyle d_{MD}^{m,n,j}=d_{Gx}^{m,n,j}+d_{Gy}^{m,n,j}$ (33)

Fourier series coefficients: $d_{i}^{m}\;\&\;d_{h}^{m,n}$ are evaluated by
multiplying sin$\left(k_{m}x\right)$ in the potential functions and
integrating from 0 to L. Where, $i=1,2,3,4,5\;\&\;h=6,7,8$.
$\lambda^{n},\Lambda^{n},\lambda^{m},\text{and}\;\Lambda^{m}$ are the Fourier
series coefficients with $n$ and $m$ indices.

$\displaystyle\lambda^{n}=\frac{2}{t_{fe}^{2}\left(k_{n}^{I}\right)^{3}}-\frac{2(-1)^{n+1}}{t_{fe}\left(k_{n}^{I}\right)^{2}}+\frac{1}{k_{n}^{I}}$
(34)
$\displaystyle\Lambda^{n}=-\frac{2(-1)^{n+1}}{\left(k_{n}^{I}\right)^{2}t_{fe}^{2}}$
(35)
$\displaystyle\lambda^{m}=\left(1+(-1)^{m+1}cos\left(k_{j}L\right)\right)\left(\frac{k_{m}}{k_{m}^{2}-k_{j}^{2}}\right)$
(36)

$\displaystyle\Lambda^{m}=(-1)^{m+1}sin\left(k_{j}L\right)\left(\frac{k_{m}}{k_{m}^{2}-k_{j}^{2}}\right)$
(37) $\displaystyle
d_{Gx}^{m,n,j}=\left(\frac{4P_{s}}{t_{fe}\epsilon_{y}d_{w}T_{D}}\right)\sum_{n}\frac{sin\left(k_{n}^{I}t_{fe}\right)}{sinh\left(k_{n}^{I}L\right)}\sum_{j}\left\\{\lambda^{n}\left(\frac{B_{j}\xi^{m,j,n}-A_{j}\zeta^{m,j,n}}{k_{j}}\right)+\Lambda^{n}\left(\frac{A_{j}\xi^{m,j,n}+B_{j}\zeta^{m,j,n}}{k_{j}^{2}}\right)\right\\}$
(38)
$\displaystyle+\left(\frac{a_{0}}{t_{fe}\epsilon_{y}}\right)\sum_{n}\frac{sin\left(k_{n}^{I}t_{fe}\right)\Lambda^{n}}{\left(k_{n}^{I}\right)^{2}}\left\\{\frac{(-1)^{m+1}}{k_{m}}-\frac{k_{m}}{k_{m}+\left(k_{n}^{I}\right)^{2}}\left(1+(-1)^{m+1}\right)\right\\}$
(39) $\displaystyle
d_{Gy}^{m,n,j}=\left(\frac{2P_{s}}{\epsilon_{x}d_{w}T_{D}}\right)\sum_{m}\frac{1}{k_{m}}\left\\{\xi^{m}(t_{fe})\sum_{j}\left(\frac{B_{j}\lambda^{m}-A_{j}\Lambda^{m}}{k_{j}}\right)+\zeta^{m}(t_{fe})\sum_{j}\left(\frac{A_{j}\lambda^{m}+B_{j}\Lambda^{m}}{k_{j}^{2}}\right)\right\\}$
(40)
$\displaystyle\xi^{m,j,n}=\frac{\left(1+(-1)^{m+1}cos\left(k_{j}L\right)\right)}{k_{j}^{2}+\left(k_{n}^{I}\right)^{2}}\left\\{\frac{k_{m}}{k_{m}^{2}-k_{j}^{2}}-\frac{k_{m}}{k_{m}^{2}+\left(k_{n}^{I}\right)^{2}}\right\\}$
(41)
$\displaystyle\zeta^{m,j,n}=\frac{(-1)^{m+1}sin\left(k_{j}L\right)}{k_{j}^{2}+\left(k_{n}^{I}\right)^{2}}\left\\{\frac{k_{m}}{k_{m}^{2}-k_{j}^{2}}-\frac{k_{m}}{k_{m}^{2}+\left(k_{n}^{I}\right)^{2}}\right\\}$
(42) $\displaystyle
A_{j}=cos\left(k_{j}d_{w}\right)-1+cos\left(k_{j}\left(x_{1}+d_{w}\right)\right)-cos\left(k_{j}\left(x_{1}+3d_{w}\right)\right)+cos\left(k_{j}\left(x_{1}+x_{2}+4d_{w}\right)\right)-+cos\left(k_{j}\left(x_{1}+x_{2}+3d_{w}\right)\right)$
(43) $\displaystyle
B_{j}=sin\left(k_{j}d_{w}\right)+sin\left(k_{j}\left(x_{1}+d_{w}\right)\right)-sin\left(k_{j}\left(x_{1}+3d_{w}\right)\right)+sin\left(k_{j}\left(x_{1}+x_{2}+4d_{w}\right)\right)-sin\left(k_{j}\left(x_{1}+x_{2}+3d_{w}\right)\right)$
(44)
$\displaystyle\phi_{\rightleftharpoons}^{fe}(x,y)=\frac{2}{L}\sum_{m=1}^{\infty}\frac{sin\left(k_{m}x\right)sinh(k_{m}y)D_{sf}^{m,n}}{\epsilon_{fe}k_{m}cosh(k_{m}t_{fe})}+\frac{2}{L}\sum_{m=1}^{\infty}sin\left(k_{m}x\right)\left(\frac{\left(1-(-1)^{m}\right)\left(V_{gs}-\Delta\phi\right)cosh\left(k_{m}\left(t_{fe}-y\right)\right)}{k_{m}cosh(k_{m}t_{fe})}\right)$
(45)
$\displaystyle\phi_{\pm}^{fe}(x,y)=\frac{2}{t_{fe}}\sum_{n}\frac{sin\left(k_{n}^{I}y\right)\left(A_{1}^{n}sinh\left(k_{n}^{I}\left(L-x\right)\right)+A_{2}^{n}sinh\left(k_{n}^{I}x\right)\right)}{sinh(k_{n}^{I}L)}$
(46)
$\displaystyle\phi^{ox}(x,y)=\phi_{\rightleftharpoons}^{ox}(x,y)+\phi_{\pm}^{ox}(x,y)+\phi^{ox}_{n=0}(x)$
(47)
$\displaystyle\phi_{\rightleftharpoons}^{ox}(x,y)=\frac{2}{t_{fe}}\sum_{n=1}^{\infty}\frac{cos\left(k_{n}^{II}\left(t_{fe}-y\right)\right)}{sinh(k_{n}^{II}L)}\left(B_{1}^{n}sinh\left(k_{n}^{II}\left(L-x\right)\right)-B_{2}^{n}sinh\left(k_{n}^{II}x\right)\right)$
(48)
$\displaystyle\phi_{\pm}^{ox}(x,y)=\frac{2}{L}\sum_{n=1}^{\infty}\frac{sin\left(k_{m}x\right)\
}{\epsilon_{ox}k_{m}sinh(k_{m}t_{ox})}\left\\{D_{sf}^{m,n}cosh\left(k_{m}\left(t_{1}-y\right)\right)-D_{sb}^{m,n}cosh\left(k_{m}\left(t_{fe}-y\right)\right)\right\\}$
(49)
$\displaystyle\phi^{ox}_{n=0}(x)=\frac{1}{2t_{ox}L}\left(t_{1}^{2}-t_{fe}^{2}\right)\left(b_{1s}\left(L-x\right)+xb_{1d}\right)+\frac{1}{2L}\left(\left(L-x\right)b_{2s}+xb_{2d}\right)$
(50) where, $A_{1}^{n},A_{2}^{n},B_{1}^{n}$, and $B_{2}^{n}$ are the Fourier
coefficients evaluated at the oxide gaps. $\displaystyle
A_{1}^{n}=\frac{\left(V_{gs}-\Delta\phi\right)}{k_{n}^{I}}+a_{1s}\frac{sin\left(k_{n}^{I}t_{fe}\right)}{\left(k_{n}^{I}\right)^{2}\left({1+\epsilon_{y}/\epsilon_{x}}\right)};\;B_{1}^{n}=b_{1s}\frac{\left((-1)^{n}-1\right)}{\left(k_{n}^{II}\right)^{2}}$
(51) $\displaystyle
a_{1s}=\frac{E_{g}/2q-V_{gs}+\Delta\phi+P(0,t_{fe}/2)t_{ox}/\epsilon_{ox}}{t_{fe}+t_{ox}\left(\epsilon_{y}/\epsilon_{x}\right)};\;b_{1s}=a_{1s}\frac{\epsilon_{y}}{\epsilon_{x}}-\frac{P_{s}}{\epsilon_{ox}};\;b_{2s}=\frac{E_{g}}{2q}-b_{1s}t_{1};\;a_{0}=\frac{2P_{s}}{T_{D}}\left(x_{a}-x_{b}\right)$
(52) $b_{1d},$ $b_{2d}$, $A_{2}^{n},\text{and}B_{2}^{n}$ are calculated by the
same manner by replacing $E_{g}/2q$ with $E_{g}/2q+V_{ds}$. $\displaystyle
k_{m}=m\pi/L;\;k_{n}^{I}=\frac{\left(2n-1\right)\pi}{t_{fe}};\;k_{n}^{II}=\frac{n\pi}{t_{ox}};\;k_{j}=\frac{2\pi
j}{T_{D}};\;t_{1}=t_{fe}+t_{ox};\;t_{2}=t_{1}+t_{si}.$
$\displaystyle\Delta\phi=\text{metal-work function, default value = 4.4 eV.}$
Green’s function for the various region are derived in our earlier work[51].