# Predictions of Electromotive Force of Magnetic Shape Memory Alloy (MSMA)
Using Constitutive Model & Generalized Regression Neural Network

Md Esharuzzaman Emu1 1Md Emu is a Graduate Student, Dept. of Mechanical
Engineering, University of Illinois at Chicago, 842 West Taylor Street
Chicago, IL 60607-7022

###### Abstract

Ferromagnetic shape memory alloys (MSMAs), such as Ni-Mn-Ga single crystals,
can exhibit the shape memory effect due to an applied magnetic field at room
temperature. Under a variable magnetic field and a constant bias stress
loading, MSMAs have been used for actuation applications. Under variable
stress and a constant bias field, MSMAs can be used in power harvesting or
sensing devices, e.g., in structural health monitoring applications. This
behavior is primarily a result of the approximately tetragonal unit cell whose
magnetic easy axis is approximately aligned with the short axis of the unit
cell within the Ni-Mn-Ga single crystals. Under an applied field, the magnetic
easy axis tends to align with the external field. Similarly, under an applied
compressive force, the short side of the unit cell tends to align with the
direction of the force.

This work introduced a new feature to the existing macro-scale magneto-
mechanical model for Ni-Mn-Ga single crystal. This model includes the fact
that the magnetic easy axis in the two variants is not exactly perpendicular
as observed by D’silva et al. [1]. This offset helps explain some of the power
harvesting capabilities of MSMAs.

Model predictions are compared to experimental data collected on a Ni-Mn-Ga
single crystal. The experiments include both stress-controlled loading with
constant bias magnetic field load (which mimics power harvesting or sensing)
and field-controlled loading with constant bias compressive stress (which
mimics actuation). Each type of test was performed at several different load
levels, and the applied field was measured without the MSMA specimen present
so that demagnetization does not affect the experimentally measured field as
suggested by Eberle et al. [2]. Results show decent agreement between model
predictions and experimental data.

Although the model predicts experimental results decently, it does not capture
all the features of the experimental data. In order to capture all the
experimental features, finally, a generalized regression neural network (GRNN)
was used to train the experimental data (stress, strain, magnetic field, &
emf) so that it can make a reasonably better prediction.

##  I INTRODUCTION

Magnetic shape memory alloys (MSMAs) are a type of smart material that
exhibits magnetic field and stress-induced shape memory effect at room
temperature. The most common alloy in this class is NiMnGa, which is the
material used in this work. The twin variants in NiMnGa alloys reorient at
low-stress levels, which changes the material’s magnetization. This phenomenon
makes MSMAs suitable for power harvesting and sensing applications [3]. Also,
the twin variants reorient under the applied field. This, along with the
relatively fast response ($\approx$ 1 kHz) of MSMAs, makes these materials
suitable for actuation [4] and micropump applications [5].
The microstructure of the material shown in Fig. 1 can explain the unique
behaviour. The magnetic easy axis is approximately aligned with the short side
of the unit cell (Fig. 1 left). When a compressive stress is applied to the
material, the short side of the unit cell tends to align with the direction of
the stress (Fig. 1, top, right). Once fully aligned, this configuration is
called the stress preferred variant. The magnetization vector M begins to
align with the direction of the external field when magnetic field H is
applied (Fig. 1 right, bottom). Once fully aligned, this configuration is
known as the field preferred variant. As can be seen in Fig. 1, transitions
between stress preferred and field preferred variants or vice versa causes a
change in magnetization vector inside the material and the bulk magnetic field
around it which can be harvested into useful voltage in conjunction with the
proper setup [6].

Figure 1: Left: Tetragonal variants in Ni-Mn-Ga samples. The solid arrows
represent the internal magnetization vectors, which are approximately aligned
with the short side of the unit cell. The direction of the short side (c) of
the unit cell defines each variant. Right: Typical behavior of MSMA under
applied stress (top) or field (bottom). Notice the short side of the unit cell
tends to align with the direction of an external compressive stress or the
direction of an external applied magnetic field.

For proper implementation of MSMA as a sensor, energy harvester, and actuator,
a generalized constitutive model is required that can predict the behavior of
MSMA subject to general magneto-mechanical loading. Several constitutive
models for the macro scale behavior of MSMAs have been proposed, including but
not limited to, those by Kiefer and Lagoudas [7], Wang and Steinmann [8],
LaMaster et al. [9], and others. But no model has been shown to satisfy all
the laws of thermodynamics and accurately capture the magnetization and
applied field (M vs. H) curves as well as the stress-strain and field-strain
response under a wide array of loading conditions. Nelson et al. [6] attempted
to predict the electromotive force (emf) generated by MSMA under certain
loading by using LaMaster’s model [9]. However, there was a noticeable
difference between predicted peak to peak emf and experimental peak to peak
emf. Further, D Silva et al. [1] observed that the field preferred and stress
preferred variants are not exactly perpendicular to each other. This offset
from perpendicular has not been included in any models.

In order to predict mechanical behavior, besides using constitutive modeling,
machine learning-based modeling is getting popular [10, 11, 12, 13]. Machine
learning-based modeling, specifically neural networks can predict nonlinear
behavior very well. In the present work, a general constitutive model is
proposed that can predict the M vs. H curves correctly and give a decent
prediction of the stress-strain and field-strain curve under various load
conditions by using the manufacturer given anisotropy constant. The offset
angle from perpendicular to the variants is also included in the model. After
that, the model is used to predict emf. Finally, a generalized regression
neural network (GRNN) [14] based model is used to predict the emf. The
predictions of emf from both models are compared.

##  II Constitutive Modeling

Among the many potential applications of MSMA, power harvesting is one of the
most promising ones. MSMAs are becoming popular in energy harvesting
applications. Their high fatigue life around $10^{9}$ cycles [15], ease of
design, particularly compared to piezo and magneto strictive material [16] and
the range of operational frequencies, 2-100 Hz [16]. In power harvesting
applications, an MSMA specimen is exposed to variable mechanical stress in
combination with a constant external magnetic field. During mechanical
loading, the internal magnetization vector of the MSMA changes. This
facilitates power generation in a pickup coil surrounded on MSMA specimen (see
Fig. 2). In the configuration shown in Fig. 2, the MSMA specimen is kept under
a constant lateral magnetic field $H_{2}$ (0.3-0.8 Tesla) with a varying
compressive load of $\sigma_{1}$ between 0 to 5 MPa to generate emf.
Previously Guiel et al. [3] found that a small constant axial magnetic field
$H_{1}$ (0.03-0.1 T) dramatically increased the emf output of the power
harvesting system shown in Fig. 2. In order to design and optimize such power
harvesting devices, a model is required to predict emf under the loading
conditions shown in Fig. 2.

Figure 2: Schematic of power harvesting setup with MSMA specimen warped by
pickup coil under varying compressive stress $\sigma_{1}$, constant axial
field $H_{1}$, and lateral field $H_{2}$.

In order to give physical grounding to this work, the predictions of emf will
be based on the model by LaMaster et al. [9]. This model is based on
sufficient conditions to satisfy the laws of thermodynamics and has been shown
to be able to predict stress-strain experiments and some field-strain
experiments reasonably well. However an offset angle observed experimentally
by D Silva et al. [1] is introduced to the model. Light microscopy on MSMA
surface reveals that the field preferred and stress preferred variants are not
exactly perpendicular to each other (see Fig. 3). The offset from
perpendicular between the two variants will be called $\beta$. The value of
$\beta$ can be arbitrarily positive or negative depending on how the specimen
is oriented as shown in Fig. 4. By performing the stress-controlled test,
Nelson et al. [6] experimentally showed that the orientations shown in Fig. 4
(a) $\&$ Fig. 4(b) produced similar emf as well as Fig. 4(c) $\&$ Fig. 4(d)
produced similar emf. Fig. 4(a, b) and Fig. 4(c, d) presume that $\beta$ is
positive and negative, respectively. This offset angle may explain the
different responses of the material with different orientations of the
specimen as represented by Nelson et al. [6].

Figure 3: Magnetic domain pattern in-plane showing magnetic domains in the
field and stress preferred variant, as viewed through a Magneto-Optical
Indicator Film (MOIF). Reproduced from [1]. Figure 4: All the possible
orientations of the MSMA specimen. TF stands for the top and front side of the
specimen, TB stands for the top and back side of the specimen. $V_{1}$ and
$V_{2}$ represent stress preferred, and field preferred variants,
respectively. The arrows indicate possible orientations of the internal
magnetization vector. The possible orientations include the offset angle
$\beta$ in $V_{1}$. Reproduced from [6].

###  II-A Model Derivation

The model is homogenized [17], which means that the model does not capture
specific features of the microstructure but rather captures the average
behaviour of the material at a point. In this case, it is also assumed that
the load on the MSMA is uniform, and thus, all points in the material in the
same sense. The model captures the average behaviour through three internal
variables which are each associated with its microstructure: magnetization
vector rotation ($\theta$), volume fraction of the magnetic domain ($\alpha$),
and volume fraction of variants reorientation ($\xi$). The model assumes that
the MSMA is consisted of only two possible variants (Fig. 5) $\xi_{1}$ and
$\xi_{2}$, which as volume fraction must sum to 1:

$\xi_{1}+\xi_{2}=1$ (1)

and are subjected to the constraint:

$0<\xi_{i}<1$ (2)

where i = 1, 2. The other two internal variables shown in Fig. 5, $\alpha_{i}$
and $\theta_{i}$ are independent of each other and bound by the following
range:

$\displaystyle 0<\theta_{i}<\frac{\pi}{2}$ (3) $\displaystyle 0<\alpha_{i}<1$
(4)

The bounds on $\alpha_{i}$ are because it is also a volume fraction. The
bounds on $\theta_{i}$ are based on physical constrain shown in Fig. 5. When
$\theta_{i}=0$ there is no rotation, and when $\theta_{i}=\frac{\pi}{2}$ the
magnetization vector has fully rotate towards the hard axis in that variant.

Figure 5: Idealized microstructure with offset angle ($\beta$) added, assuming
only two variants, distinct domain walls in each variant and magnetization
vector rotation.

With the offset angle $\beta$ added, the rotation of the internal
magnetization vector of one of the variants will be changed, as the no field
applied position of the internal magnetization vector changes. Fig. 5, shows a
schematic of Cartesian component of the magnetization vector, which can be
written as

$\displaystyle\textbf{{M}}=M_{s}\,\begin{Bmatrix}\xi_{2}\,\sin\left(\theta_{2}\right)+\xi_{1}\,{\left(\alpha_{1}-1\right)}\,\cos\left(\theta_{1}-\beta\right)\,+\xi_{1}\,\alpha_{1}\,\cos\left(\theta_{1}+\beta\right)\\\
\xi_{1}\,\alpha_{1}\,\sin\left(\theta_{1}+\beta\right)+\xi_{1}\,(1-\alpha_{1})\,\sin\left(\theta_{1}-\beta\right)+\xi_{2}\,\left(2\,\alpha_{2}-1\right)\,\cos\left(\theta_{2}\right)\end{Bmatrix}$
(5)

where $M_{s}$ is the magnetization saturation value given by the manufacturer.
Obeying the rules of thermodynamics and considering mechanical and magnetic
energy sources, the Gibbs free energy takes the following form:

$\displaystyle\begin{split}g\left(\boldsymbol{\sigma},\boldsymbol{H^{app}},\xi_{i},\alpha_{i},\theta_{i}\right)=\frac{1}{2\rho}\boldsymbol{\sigma}:\boldsymbol{S\sigma}-\frac{\mu_{0}}{\rho}\boldsymbol{M.H}\\\
+g^{an}+\frac{1}{\rho}f^{\xi}+\frac{1}{\rho}f^{\alpha}+\frac{1}{\rho}f^{\theta}\end{split}$
(6)

where $\boldsymbol{\sigma}$ is stress, S is the compliance tensor, M is the
magnetization vector, H is the magnetic field, $g^{an}$ is the anisotropy
energy, and $f^{\xi}$, $f^{\alpha}$, $f^{\theta}$ are assumed hardening
functions. These hardening functions capture how the evolution of the internal
variables effect energy storage in the material. These terms aim to capture
micro-structural interactions in a phenomenological sense. In this model, the
values of $f^{\alpha}$ and $f^{\theta}$ are assumed to be zero. The $f^{\xi}$
term is used to help predict changes seen in the stress-strain and field-
strain curves during reorientation.
The strain $\boldsymbol{\varepsilon}$ in MSMA can be expressed as

$\displaystyle\boldsymbol{\varepsilon}=\boldsymbol{\varepsilon^{e}}+\boldsymbol{\varepsilon^{r}}$
(7)

where $\boldsymbol{\varepsilon^{e}}$ is the elastic strain tensor and
$\boldsymbol{\varepsilon^{r}}$ is the reorientation strain tensor.
The reorientation strain tensor $\boldsymbol{\varepsilon^{r}}$ can be written
as

$\displaystyle\boldsymbol{\varepsilon^{r}}=\begin{bmatrix}\xi_{1}&0\\\
0&\xi_{2}\end{bmatrix}\varepsilon^{r,max}$ (8)

where $\varepsilon^{r,max}$ is the maximum achievable reorientation strain
obtained from the experimental data.
The compliance tensor S can be expressed as

$\displaystyle\boldsymbol{S}=\boldsymbol{S}^{\xi_{1}}\xi_{1}+\boldsymbol{S}^{\xi_{2}}\xi_{2}$
(9)

where $S^{\xi_{1}}$ and $S^{\xi_{2}}$ are the compliance tensors for variants
1 and 2, respectively.
Following [18, 19], and neglecting the higher order term, the anisotropy
energy $g^{an}$ [20] can be written as

$\displaystyle
g^{an}=\xi_{1}\,K_{1}\,\sin^{2}\,\theta_{1}+\xi_{2}\,K_{1}\,\sin^{2}\,\theta_{2}$
(10)

where $K_{1}$ is the anisotropy coefficient which is the measurement of of the
magnetization vectors’ resistance to rotation.
The internal magnetic field is assumed as follows:

$\textbf{{H}}=\textbf{{H}}^{app}-\textbf{{D}}\,\textbf{{M}}$ (11)

where H is the internal magnetic field experienced by the MSMA specimen,
$H^{app}$ is the external magnetic field applied on the specimen, D is the
demagnetization tensor [2].
Using Eq. (7) $\&$ (8) in Eq. (6), the Clausius-Duhem inequality can be
written as

$\begin{split}\left(\sigma_{11}\,\varepsilon^{r,max}-\rho\frac{\partial
g}{\partial\xi_{1}}\right)\dot{\xi_{1}}+\left(\sigma_{22}\,\varepsilon^{r,max}-\rho\frac{\partial
g}{\partial\xi_{2}}\right)\dot{\xi_{2}}\\\ -\rho\frac{\partial
g}{\partial\theta_{1}}\dot{\theta_{1}}-\rho\frac{\partial
g}{\partial\theta_{2}}\dot{\theta_{2}}-\rho\frac{\partial
g}{\partial\alpha_{1}}\dot{\alpha_{1}}-\rho\frac{\partial
g}{\partial\alpha_{2}}\dot{\alpha_{2}}\,\geq\,0\end{split}$ (12)

Introducing $\pi$ as the driving force for change of each internal variable,
Eq. (12) can be written as

$\pi^{\xi_{1}}\dot{\xi_{1}}+\pi^{\xi_{2}}\dot{\xi_{2}}+\pi^{\theta_{1}}\dot{\theta_{1}}+\pi^{\theta_{2}}\dot{\theta_{2}}+\pi^{\alpha_{1}}\dot{\alpha_{1}}+\pi^{\alpha_{2}}\dot{\alpha_{2}}\,\geq\,0$
(13)

where the $\pi$ terms are defined as

$\displaystyle\pi^{\xi_{1}}=\left(\sigma_{11}\,\varepsilon^{r,max}-\rho\frac{\partial
g}{\partial\xi_{1}}\right)\,\,$ (14)
$\displaystyle\pi^{\xi_{2}}=\left(\sigma_{22}\,\varepsilon^{r,max}-\rho\frac{\partial
g}{\partial\xi_{2}}\right)$ (15)
$\displaystyle\pi^{\theta_{1}}=-\rho\frac{\partial
g}{\partial\theta_{1}},\,\,\pi^{\theta_{2}}=-\rho\frac{\partial
g}{\partial\theta_{2}}$ (16)
$\displaystyle\pi^{\alpha_{1}}=-\rho\frac{\partial
g}{\partial\alpha_{1}},\,\,\pi^{\alpha_{2}}=-\rho\frac{\partial
g}{\partial\alpha_{2}}$ (17)

In this model, $\alpha_{i}$ is piece-wise defined. When the applied field
reaches a certain threshold, the value of alpha is considered 1, and less than
that value of alpha is 0.5. This threshold value is called $H_{cri}^{\alpha}$,
which is calculated from the M-H curve given by the manufacturer. The equation
of alpha can be expressed as follows:

$\alpha_{i}=\begin{cases}0.5,&\text{if }H_{i}<H_{cri}^{\alpha}\\\
1,&\text{otherwise}\end{cases}$ (18)

As $\alpha_{i}$ is piece-wise constant in the Eq. (18), $\dot{\alpha_{i}}$ is
zero almost everywhere. Thus, $\pi^{\alpha_{i}}\dot{\alpha_{i}}$ is also zero,
which satisfies the thermodynamics. However, there is a possibility of the
thermodynamics being violated at the single instant where $\alpha_{i}$
changes, but as it is only a single instant, this potential violation of
thermodynamic requirements will be ignored.

Following [17], the hardening function $f^{\xi}$ used in this model is

$f^{\xi}=\frac{1}{6}C_{1}\,\xi_{1}^{6}+\frac{1}{2}C_{2}\,\xi_{1}^{2}$ (19)

where $C_{1}$ & $C_{2}$ are constants that need to be calibrated.

Finally, the model is completed by using following Kuhn-Tucker conditions to
ensure positive dissipation during variant reorientation:

$\begin{cases}\pi^{\xi_{1}}-\pi^{\xi_{2}}=Y,&\dot{\xi_{1}}>0\\\
\pi^{\xi_{1}}-\pi^{\xi_{2}}=-Y,&\dot{\xi_{1}}<0\end{cases}$ (20)

where Y is a positive constant and calibrated from experimental data.

###  II-B Model Calibration

The model needs to be calibrated with the experimental data to find the
hardening constants $C_{1}$, $C_{2}$, and the positive constant of the Kuhn-
Tucker condition Y. Following [9], the simplified loading condition is used
for calibration. For calibration, the specimen is brought into a fully
elongated state by applying a lateral field. Then the lateral field is
removed, and the specimen is gradually compressed to 4.5 MPa no field applied.
As there is no magnetic field, the magnetic domains $\alpha_{i}$ are 0.5, and
there is no variant rotation, i.e., $\theta_{i}=0$. Using these values in Eq.
(20) and considering $\dot{\xi_{1}}>0$, as the specimen is transitioning from
variant 2 to variant 1, the following equation is obtained:

$Y+2\,C_{1}\,{\xi_{1}^{5}}_{i}+2\,C_{2}\,{\xi_{1}}_{i}=e^{r,max}{\sigma_{11}}_{i}+{\sigma_{11}^{2}}_{i}(S^{\xi_{1}}-S^{\xi_{1}})$
(21)

where i=1, 2, 3; ${\sigma_{11}}_{i}$ and ${\xi_{1}}_{i}$ are the compressive
stress and variant volume fraction of the $i-{th}$ point on the experimental
stress-strain curve with zero field shown in Fig. 6(a).

(a) Experimental curve

(b) Predicted curve

Figure 6: Calibration test with uni-axial load and zero field, experimental
stress-strain curve is on the left side (a), black dots are the points for
calculating the compliance tensor and the red three dots on the reorientation
region are the points for calculating the constants $C_{1},C_{2},Y$, predicted
stress-strain curve (solid red line) is shown on the right side (b).

The variant volume fraction (${\xi_{1}}_{i}$) shown in Eqn. (21) can be found
from the experimental data, using the Eqns. (7), (8) $\&$ (9), as

${\xi_{1}}_{i}=\frac{{\varepsilon_{11}}_{i}-S^{\xi_{2}}\,{\sigma_{11}}_{i}}{(S^{\xi_{1}}-S^{\xi_{2}}){\sigma_{11}}_{i}\,\varepsilon^{r,max}}$
(22)

The compliance tensor $S^{\xi_{i}}$ is calculated by taking the inverse of the
modulus of elasticity, which is measured directly from the experimental
stress-strain curve. Two different moduli of elasticity are found; one is in
variant 2 (fully elongated state), and the other is in variant 1 (fully
compressed state) from the two distinct elastic regions in Fig. 6(a).
Specifically, the two points in each elastic region (see Fig. 6(a), black
dots) are used to calculate the slope, which gives the modulus of elasticity.
The maximum reorientation $\varepsilon^{r,max}$ is calculated from the
experimental data by using the following equations:

$\varepsilon^{r,max}=\varepsilon_{min}-\frac{\sigma_{min}}{E_{1}}$ (23)

where $\varepsilon_{min}$ is the maximum negative strain of the data set,
$\sigma_{min}$ is the maximum compressive stress of the experimental data, and
$E_{1}$ is the modulus of elasticity of the variant 1. The origins of this
equation are shown in Fig. 7.

Figure 7: Schematic of measuring maximum reorientation strain
$\varepsilon^{r,max}$.

Finally, all the parameters of Eq. (21) are known except $C_{1}$, $C_{2}$ $\&$
Y. These constants can be calculated by solving a system of three linear
equations, which is obtained by taking three points on the calibration curve
(see Fig. 6(a), red dots). These calibrated constants are listed in Table II.

TABLE I: Material Parameters. $\rho\,K_{1}$ | $0.89\text{\times}{10}^{5}\text{\,}\mathrm{J}\mathrm{/}\mathrm{m}^{3}$
---|---
$M_{s}$ | $5.44\text{\times}{10}^{5}\text{\,}\mathrm{A}\text{\,}{\mathrm{m}}^{-1}$
${H^{\alpha}}_{cri}$ | $23\,000\text{\,}\mathrm{A}\text{\,}{\mathrm{m}}^{-1}$
$\mu_{0}$ | 4$\pi$$\times 10^{7}$ $Hm^{-1}$
$D_{11}$ | 0.0674
$D_{22}$ | 0.466
TABLE II: Calibrated Constants. $C_{1}$ | 39408.25
---|---
$C_{2}$ | 9962.91
$Y$ | 44181.38
$S^{\xi_{1}}$ | $1.957\text{\times}{10}^{-10}\text{\,}\mathrm{m}^{2}\mathrm{/}\mathrm{N}$
$S^{\xi_{2}}$ | $2.247\text{\times}{10}^{-9}\text{\,}\mathrm{m}^{2}\mathrm{/}\mathrm{N}$
$e^{r,max}$ | -0.057867

Using the material parameters and constants listed in Table I $\&$ II and the
derived model, the simulation is run, and the predicted calibration curve is
shown in Fig. 6(b). As expected, the calibration curve is well matched by the
simulation and the simulation exactly matches the experiment at the red points
used for calibration.

###  II-C Introducing EMF to the Model

From Faraday’s law of induction, the following equations can be used to
calculate emf:

$E=-NA\,\frac{d\,B}{d\,t}\\\ $ (24) $B=\mu_{0}\left(H+M\right)$ (25)

where E is emf, N is the number of turns in the coil, A is cross-section area,
$\mu_{0}$ is magnetic permeability constant, and t is time.

Using Eq. (11), (24) & (25), emf can be simplified to

$E=-NA\,\mu_{0}\left(1-D_{11}\right)\frac{d\,M_{1}}{d\,t}$ (26)

where $D_{11}$ is the 11 component of demagnetization tensor, and $M_{1}$ is
the magnetization vector component in 1 direction shown in Fig. 2 obtained
from the model. $\frac{d\,M_{1}}{d\,t}$ can be calculated numerically by
assuming

$\frac{d\,M_{1}}{d\,t}\simeq\frac{\Delta\,M_{1}}{\Delta\,t}$ (27)

where $\Delta\,M_{1}$ is the difference in magnetization between two
consecutive time steps, and $\Delta\,t$ is the time difference between those
steps.

##  III Experiments

The MSMA specimen under the experimental configuration shown in Fig. 2 was
used to gather emf data from the coil. Several experiments were conducted by
Nelson et al. [6] and Guiel et al. [3] to get voltage output from the MSMA
specimen. Nelson et al. [6] performed the stress-controlled test to be
consistent with the LaMaster’s model [9]. But instead of using the axial
field, Nelson et al. [6] tilted the lateral field at a certain angle. The
tilted lateral field was not ideal because, in their setup, it was not
possible to get a uniform magnetic field over the MSMA specimen. Moreover,
more recent experimental result with a similar setup suggested that these
experiments might have had significant variability of the stress and this
input might not have been controlled as well as originally thought.

On the other hand, Guiel et al. [3] performed the strain-controlled test,
i.e., the Instron Machine was kept in position-controlled mode. In this setup,
the maximum strain was measured experimentally before each cycle. But in this
case, the specimen was experienced higher compressive stress ($\approx$ 8
MPa), which caused cracks in the specimen.

For several reasons these experiments could not be directly compared with
model predictions. First, the Guiel et al. experiments were strain-controlled
while the model is stress controlled. Next, the field was measured at the
surface in these tests, which accordingly to Eberle et al. [2] should not be
used as the input applied field. Finally, the Guiel et al. experiments, which
were performed at 15 Hz frequency showed evidence that inertial effects were
evident (i.e. the strain obtained from the tests at 15 Hz seemed different
than the strain obtained from those tests done very slowly at 2 Hz).

Therefore, to be consistent with the model, additional stress-controlled
experiments under the conditions shown in Fig. 2 were performed at a slower
rate of 10 Hz in order to compare the results with the model. At this
frequency level, the strain seems consistent with the strain obtained from the
test of 2 Hz frequency. The power harvesting experimental data used to
validate test this model was acquired through stress-controlled tests that
were performed with improved PID control parameters, correcting the errors
encountered in the data reported by Nelson et al. [6]. During these tests, the
axial compressive stress was varied between 0 to 5 MPa, while the specimen is
exposed to constant axial and lateral fields.

As shown by Guiel et al. [3], the maximum emf can be obtained with a certain
angle ($\phi$) of the resultant magnetic field with the horizontal axis (see
Fig 8). The angle $\phi$ is calculated by using the following equation:

$\phi=\tan^{-1}\left(\frac{H_{1}}{H_{2}}\right)$ (28)

The resultant field $H_{res}$ is calculated as:

$H_{res}=\sqrt{H_{1}^{2}+H_{2}^{2}}$ (29)

Figure 8: Schematic of MSMA with magnetic field and angle. $H_{2}$ is the
lateral field, $H_{1}$ is the axial magnetic field, $H_{res}$ is the resultant
magnetic field, and $\phi$ is the angle between the lateral field and the
resultant field.

The emf has been measured experimentally at several angles ranging from 0 to
7.2$\mathrm{\SIUnitSymbolDegree}$ while keeping the resultant field unchanged
at 0.578 T. Using these angles and associated field values, simulations have
been performed to calculate emf from the revised model.

The experimental setup for power harvesting is shown in Fig. 9. The axial
field ($H_{1}$) was applied using the permanent magnets as shown in Fig. 9,
lateral field ($H_{2}$) was applied by using an electromagnet shown in Fig. 9.
The magnitude of these magnetic fields were measured by using a Hall probe and
the angle was calculated by using Eq. (28). The compressive stress was applied
to the sample through the grip shown in Fig. 9 by using an 8874 Instron
testing machine. The machine was kept in load controlled mode and the applied
forces on the specimen were recorded via Instron Machine associated data
acquisition software Wavematrix 2. The force applied on the specimen and the
displacement of the top part of the sample were recorded automatically by the
software controlling the machine. The strain could not be measured during the
experiment with the available equipment because the coil surrounded the sample
(see Fig. 9). The strain could be calculated, although not very accurately,
from the position data collected. Before running the experiment, a white dot
was placed on one of the sides of the specimen through which the lateral filed
was applied. When the dot was kept in the top left side of the specimen, this
configuration was assumed as $+\beta$, and when the dot was held at the top
right side of the specimen, the configuration was assumed as $-\beta$.

Figure 9: Experimental setup for the power harvesting test. Axial field
$H_{1}$ is applied on the specimen by using permanent magnet and lateral field
$H_{2}$ is applied by using electromagnet. Different values of $\phi$ are
obtained by varying $H_{1}$ $\&$ $H_{2}$. Note that $H_{res}$ is always kept
same.

##  IV Generalized Regression Neural Network

Generalized regression neural network (GRNN) which was proposed by D.F. Specht
[14], was adapted from radial basis neural networks [21] so that it can be
used for regression, prediction, and classifications.

Figure 10: Simplified architecture of Radial Basis Function Network
(RBFN).Here $x_{n}$ is the input vector and $\mu$ is the prototype. The
prototype vector may be defined as the center of the neuron since the
prototype is the value of the center of the bell curve.

GRNN represents an enhanced technique in the neural networks on the basis of
non-parametric regression [22] which bounds every training sample to represent
a mean to a radial basis neuron. Mathematically it can be expressed as

$Y(x)=\dfrac{\sum_{k=1}^{N}y_{k}K(x,x_{k})}{\sum_{k=1}^{N}K(x,x_{k})}$ (30)

where $Y(x)$ is the predicted value of input $x$, $y_{k}$ is the activation
weight for the pattern layer neuron at $k$, $K(x,x_{k})$ is the Radial basis
function kernel (RBFK) [23]. Here Gaussian kernel is used as RBFK which can be
expressed as

$K(x,x_{k})=e^{-d_{k}/{2\sigma^{2}}}$ (31)

where $d_{k}$ is the squared euclidean distance between the training samples
($x_{k}$) and the input ($x$) which is expressed as

$d_{k}=(x-x_{k})^{T}(x-x_{k})$ (32)

Eq. 30 can be modified as [14]:

$E(F,H_{1},H_{2})=\dfrac{\sum_{k=1}^{N}A^{k}e^{-d_{k}/{2\sigma^{2}}}}{\sum_{k=1}^{N}B^{k}e^{-d_{k}/{2\sigma^{2}}}}$
(33)

where E indicates the estimator of emf, $A^{k}$ $\&$ $B^{k}$ represent the
coefficients for the cluster.

Figure 11: Simplified architecture of Generalized Regression Neural Network
(GRNN) for predicting emf.

In this work, the GRNN is chosen for the prediction because of its advantages
over other neural networks. In our study, the data set is small and GRNN is
efficient working with less volume of data and it utilizes Gaussian functions
which increase the prediction accuracy. Moreover, GRNN does not require any
backpropagation to update the weights which makes it possible to avoid the
limitation of backpropagation [24]. In Fig. 11 a feed-forward neural network
is shown which can be used to estimate the emf E from the measurement vector,
force (F), and magnetic field (H).

Figure 12: Simplified Flow Chart of Machine Learning Based Modeling.

In Fig. 12, a simplified workflow chart of the machine learning-based
modelling is shown. We had a total of four variables, axial force ($F$), axial
magnetic field ($H_{1}$), lateral magnetic field ($H_{2}$), and emf, among
these three, are independent variables and one is the dependent variable
(emf). 75$\%$ of the total data was used for training the model and the rest
25$\%$ was used to test the model using the builtin test-train split function
of scikit-learn. All the values of the data set were scaled in-between 0 to 1
to increase the model efficiency. k-fold cross-validation is used in this
model where the original sample is randomly separated into k equal-sized
subgroups. Among k subgroups, a single group is picked as the validation data
for testing the model, and the rest $(k-1)$ subgroups are used as training
data.

##  V Results

The experimental results shown here were obtained using the loading frequency
of 10 Hz.

Figure 13: Peak to peak emf vs. angle plots. The solid blue curve represents
the experimental emf with positive $\beta$. The dashed red, black, & orange
curve represent the model predicted emf with modified model, LaMaster et. al.
model [9], & GRNN model, respectively.

The predicted curve shown in Fig. 13 follows the same trend as the experiment,
suggesting that this model reasonably predicts the angle for maximum emf.
However, discrepancies between the model and experiment suggest that further
improvement is needed in constitutive modeling. However, GRNN predicts the
experimental curve very well. In Fig. 14, experimental and machine learning
model predicted emf is compared with a function of time. The GRNN model
captures the experimental curve with a accuracy of 99$\%$.

Figure 14: Emf vs. time plots. The solid blue curve represents the
experimental emf with positive $\beta$ while the dashed orange curve
represents the model predicted emf with GRNN model. The lateral and axial
magnetic fields were 0.575 & 0.057 T, respectively.

##  VI Discussion

The primary purpose of this study was to accurately predict the open circuit
emf output of MSMA in power harvesting applications. This prediction can be
used to design and optimize power harvesting devices. Towards that end, emf
was predicted from the LaMaster et al. [9], a modified model where an
experimentally observed offset angle was added to the LaMaster et al. model,
and generalized regression neural network (GRNN) model. All models were used
to predict experimental emf data and the results are shown in Fig. 13. The
experimental emf shown in Fig. 13 &14 were obtained by performing a load
control test with a magnetic field applied at various angles. Note that the
models all assume that the applied field was uniform, but that was hard to
achieve experimentally.

One of the major contributions of this study is to add an offset angle to the
model which helps to explain some of the power-harvesting capabilities of
MSMAs. Using Fig. 15, the reason for getting different peak to peak emf from
different orientations [6] of the MSMA can be hypothesized. There seem to be
competing factors explaining how $\beta$ affects emf output. On the one hand,
Fig. 15(a) with $\beta$ positive seems to make reorientation easier. As the
$\xi_{i}$ vs. time plot in Fig. 16 shows, there is more reorientation, and the
slope is stiffer with $+\beta$ than with $-\beta$ because of the requirements
of less energy to reorient. Thus, change in magnetization within a certain
time period is more in Fig. 16(a) than Fig. 16(b). As the difference in
magnetization is higher, so the emf is also higher (Eq. (26)). On the other
hand, in Fig. 15(b) there is more overall change in magnetization as the
specimen changes from the stress preferred state to the field preferred state.
Thus, when the full reorientation is obtained (see Fig. 17), $-\beta$ gives
more emf than $+\beta$. In this case, the stress and field level is very high,
more than 8 MPa and 0.7 T, respectively. As the variants are fully reoriented,
there is no additional reorientation of the case of $+\beta$ and instead, the
more change in magnetization in $-\beta$ case leads to more emf with $-\beta$.

Figure 15: Schematic of MSMA with two different orientations, i.e., positive
$\beta$ (a) and negative $\beta$ (b). Constant lateral ($H_{2}$) and axial
($H_{1}$) magnetic field are applied on the specimen. Varying axial
compressive stress ($\sigma_{11}$) is also applied. $V_{1}$ and $V_{2}$
represent stress preferred, and field preferred variants, respectively. Figure
16: Emf and variant volume fraction prediction with associate load condition
from the model with offset angle added for both (a) positive and (b) negative
$\beta$. Peak to peak voltage in (a) $\&$ (b) are 0.371 $\&$ 0.315 volts,
respectively. Note that full reorientation is not obtained here due to low
axial stress (0-4.8 MPa) and resultant field (0.578 T) level. The axial field
$H_{1}$ is 0.029 T. Figure 17: Emf and variant volume fraction prediction with
associate load condition from the model with offset angle added for both (a)
positive and (b) negative $\beta$. Peak to peak voltage in (a) $\&$ (b) are
0.5131 $\&$ 0.6021 volts, respectively. Note that full reorientation is
obtained here due to high axial stress (0-8.3 MPa) and resultant field (0.7 T)
level. The axial field $H_{1}$ is 0.018 T.

In this work, the experiments performed were stress-controlled, and the stress
varies between 0 to 4.8 MPa. Therefore, full reorientation was not obtained.
Thus, the ease of reorientation dominates, and as seen earlier, the $+\beta$
configuration (see Fig. 15(a)) produces more emf than $-\beta$ configuration
(see Fig. 15(b)). The configuration which produces more emf experimentally
than the other is assumed as $+\beta$ configuration in Fig. 13. The idea that
$-\beta$ produces more emf with very large fields and stresses was not
verified experimentally, however, it should be as part of future work. The
hypothesizes about the sign of $\beta$ and the maximum emf suggests that two
different types of energy harvester are possible to design, one is for low
stress level, and another one is for high stress level or strain-controlled.
In either case, Fig. 13 suggests that the peak to peak emf is maximum at an
angle between 3-7$\mathrm{\SIUnitSymbolDegree}$ and this power harvesting
devices of the type in Fig. 9 should have a small axial field to maximize emf
output. However, as the angle depends on the magnitude of the axial and
lateral field, it is necessary to pick suitable magnets for maximizing the emf
output.

In Fig. 13 the predicted emf obtained from the modified model (red dashed
line) falls closer to the experimental emf than the LaMaster et al. model
(black dashed line). This suggests that adding a new feature, i.e, offset
angle to the model increases the model accuracy. However, there is still a
significant difference between experimental and the predicted results which
suggests that some key features are still missing from the constitutive model.
Model calibration can be the possible reason behind the inconsistency of poor
prediction. All the constitutive models were calibrated at zero axial fields
while the experiments were performed with the presence of different axial
field values. This may be the possible reason for the deviation of predicted
trend extremely with increasing of the angle (Fig. 13).

In order to overcome all the constrain of constitutive modeling, a generalized
regression neural network (GRNN) model is implemented in this study to make
accurate predictions. GRNN model is extensively dependent on experimental data
since the model learns from the data and finds the hidden co-relation. In our
study, the GRNN model was able to predict the experimental curve with a high
level of accuracy. Besides capturing peak to peak emf (Fig. 13), GRNN model
can predict emf as a function of time (Fig. 14) exceptionally well.

##  VII CONCLUSIONS

Since the constitutive models in this work were not particularly accurate
enough at predicting emf, their use in design is limited, and significant
future work is needed to find the missing features. Some future work is needed
to test some of the hypothesizes and understanding of the models in this work
regarding the emf at the high-stress level. Further study is also recommended
to find any other form of energy that may be included in Gibbs free energy to
improve the prediction of variant reorientation. One such kind of energy is
magnetic exchange energy which is the fundamental property of ferromagnetic
materials that favors alignment of the magnetization vector spins along a
common direction [25]. Although constitutive modeling does not produce good
predictions, the GRNN model was successful to predict the behavior and these
machine learning approaches can be useful in designing MSMA-based sensors and
power harvesters.

## ACKNOWLEDGMENT

This work has been supported by the National Science Foundation under Grant
No. 1101108 and Grant No. 1561866. Any opinions, findings, and conclusions or
recommendations expressed in this article are those of the author(s) and do
not necessarily reflect the views of the National Science Foundation. The
author would like to thank Glen D’Silva, for his constructive feedback and
helping him doing the experiments. The author would also like to thank
Professor Sybil Derrible for helping him to implement the machine learning
algorithm.

## References

  * [1] G. J. D. Silva, C. Ciocanel, and H. P. Feigenbaum, “Visualization of magnetic domains and magnetization vectors in magnetic shape memory alloys under magneto-mechanical loading.,” Shape memory and superelasticity, vol. 6, March 2020.
  * [2] J. L. Eberle, H. P. Feigenbaum, and C. Ciocanel, “Demagnetizing field in single crystal ferromagnetic shape memory alloys,” Smart Materials and Structures, vol. 28, p. 025022, jan 2019.
  * [3] R. Guiel, H. Feigenbaum, and C. Ciocanel, “The effect of magnetic field orientation on the open-circuit voltage of Ni-Mn-Ga based power harvesters,” Smart Materials and Structures, vol. 27, jul 2018.
  * [4] B. Minorowicz, G. Leonetti, F. Stefanski, G. Binetti, and D. Naso, “Design, modelling and control of a micro-positioning actuator based on magnetic shape memory alloys,” Smart Materials and Structures, vol. 25, p. 075005, may 2016.
  * [5] K. Ullakko, L. Wendell, A. Smith, P. Müllner, and G. Hampikian, “A magnetic shape memory micropump: Contact-free, and compatible with PCR and human DNA profiling,” Smart Materials and Structures, vol. 21, no. 11, 2012.
  * [6] I. Nelson, J. Dikes, H. Feigenbaum, and C. Ciocanel, “Numerical predictions versus experimental findings on the power-harvesting output of a NiMnGa alloy,” Behavior and Mechanics of Multifunctional Materials and Composites 2014, vol. 9058, no. March 2014, p. 905815, 2014.
  * [7] B. Kiefer and D. C. Lagoudas, “Magnetic field-induced martensitic variant reorientation in magnetic shape memory alloys,” Philosophical Magazine, vol. 85, no. 33-35, pp. 4289–4329, 2005.
  * [8] J. Wang and P. Steinmann, “A variational approach towards the modeling of magnetic field-induced strains in magnetic shape memory alloys,” Journal of the Mechanics and Physics of Solids, vol. 60, no. 6, pp. 1179–1200, 2012.
  * [9] D. H. LaMaster, H. P. Feigenbaum, I. D. Nelson, and C. Ciocanel, “A full two-dimensional thermodynamic-based model for magnetic shape memory alloys,” Journal of Applied Mechanics, Transactions ASME, vol. 81, no. 6, pp. 1–12, 2014.
  * [10] M. Mozaffar, R. Bostanabad, W. Chen, K. Ehmann, J. Cao, and M. A. Bessa, “Deep learning predicts path-dependent plasticity,” Proceedings of the National Academy of Sciences of the United States of America, vol. 116, no. 52, pp. 26414–26420, 2019.
  * [11] G. A. Holzapfel, K. Linka, S. Sherifova, and C. J. Cyron, “Predictive constitutive modelling of arteries by deep learning,” Journal of the Royal Society Interface, vol. 18, no. 182, 2021.
  * [12] K. Linka, M. Hillgärtner, K. P. Abdolazizi, R. C. Aydin, M. Itskov, and C. J. Cyron, “Constitutive artificial neural networks: A fast and general approach to predictive data-driven constitutive modeling by deep learning,” Journal of Computational Physics, vol. 429, p. 110010, 2021.
  * [13] Z. Zhang, Q. Liu, and D. Wu, “Predicting stress–strain curves using transfer learning: Knowledge transfer across polymer composites,” Materials and Design, vol. 218, p. 110700, 2022.
  * [14] D. Specht, “A general regression neural network,” IEEE Transactions on Neural Networks, vol. 2, no. 6, pp. 568–576, 1991.
  * [15] J. Tellinen, I. Suorsa, I. Aaltio, and K. Ullakko, “Basic Properties of Magnetic Shape Memory Actuators,” 8th international conference ACTUATOR 2002, pp. 10–12, June 2002.
  * [16] A. J. Niskanen and I. Laitinen, “Design and simulation of a magnetic shape memory (msm) alloy energy harvester,” in State-of-the-Art Research and Application of SMAs Technologies (4th CIMTEC), vol. 78 of Advances in Science and Technology, pp. 58–62, Trans Tech Publications Ltd, January 2013\.
  * [17] J. L. Eberle, “Predicting the magneto-mechanical behavior of single crystal magnetic shape memory alloys using homogenized models.,” Master’s thesis, Northern Arizona University, 2018.  ProQuest13423996.
  * [18] B. Kiefer and D. Lagoudas, “Modeling the coupled strain and magnetization response of magnetic shape memory alloys under magnetomechanical loading,” Journal of Intelligent Material Systems and Structures, vol. 20, pp. 143–170, 01 2009.
  * [19] A. Waldauer, H. Feigenbaum, C. Ciocanel, and N. Bruno, “Improved thermodynamic model for magnetic shape memory alloys,” Smart Materials and Structures \- SMART MATER STRUCT, vol. 21, 09 2012.
  * [20] S. P. Farrell, R. A. Dunlap, L. M. Cheng, R. Ham-Su, M. A. Gharghouri, and C. V. Hyatt, “Magnetic properties of single crystals of Ni-Mn-Ga magnetic shape memory alloys,” Smart Structures and Materials 2004: Active Materials: Behavior and Mechanics, vol. 5387, no. July 2004, p. 186, 2004.
  * [21] D. BROOMHEAD, “Multivariable functional interpolation and adaptive networks,” Complex Systems, vol. 2, pp. 321–355, 1988.
  * [22] A. Bowman and A. Azzalini, Applied Smoothing Techniques for Data Analysis: The Kernel Approach with S-Plus Illustrations.  Oxford Statistical Science Series, OUP Oxford, 1997.
  * [23] J. Vert, K. Tsuda, and B. Schölkopf, A Primer on Kernel Methods, pp. 35–70.  Cambridge, MA, USA: MIT Press, 2004.
  * [24] Y. Lecun, Y. Bengio, and G. Hinton, “Deep learning,” Nature, vol. 521, no. 7553, pp. 436–444, 2015.
  * [25] C. J. García-Cervera, “Numerical micromagnetics: A review,” 2007.


# User Pairing and Outage Analysis in Multi-Carrier NOMA-THz Networks

Sadeq Bani Melhem, and Hina Tabassum S. B.Melhem and H. Tabassum are with the
York University, Canada (e-mail<EMAIL_ADDRESS>and hina@eecs.yorku.ca). This
work is supported by the Discovery Grant from the Natural Sciences and
Engineering Research Council of Canada.

###### Abstract

This paper provides a comprehensive framework to analyze the performance of
non-orthogonal multiple access (NOMA) in the downlink transmission of a
single-carrier and multi-carrier terahertz (THz) network. Specifically, we
first develop a novel user pairing scheme for the THz-NOMA network which
ensures the performance gains of NOMA over orthogonal multiple access (OMA)
for each individual user in the NOMA pair and adapts according to the
molecular absorption. Then, we characterize novel outage probability
expressions considering a single-carrier and multi-carrier THz-NOMA network in
the presence of various user pairing schemes, Nakagami-$m$ channel fading, and
molecular absorption noise. We propose a moment-generating-function (MGF)
based approach to analyze the outage probability of users in a multi-carrier
THz network. Furthermore, for negligible thermal noise, we provide simplified
single-integral expressions to compute the outage probability in a multi-
carrier network. Numerical results demonstrate the performance of the proposed
user-pairing scheme and validate the accuracy of the derived expressions.

###### Index Terms:

THz, NOMA, user-pairing, multi-carrier, outage probability, molecular
absorption.

##  I Introduction

Fueled by the emergence of machine-type communications in a variety of
wireless applications, the provisioning of massive connectivity becomes
instrumental. On the other hand, accommodating trillions of devices within the
extremely congested and limited sub-6GHz spectrum is becoming challenging. In
this context, shifting to higher frequency terahertz (THz) communication is
under consideration to obtain the data rates in the order of hundreds of Gbps
[1]. Also, non-orthogonal multiple access (NOMA) is becoming popular to
support multiple users in the same frequency and time while leveraging on
efficient interference cancellation mechanisms.

To date, most of the research focuses on analyzing the outage probabilities
(OPs) of users in a single-carrier sub-6GHz NOMA network [2, 3, 4]. An
exception is [5] where the authors considered the OP analysis for multi-
carrier NOMA in sub-6GHz or radio frequency (RF) networks. The derived
expressions rely on approximations and are in the form of Fox H’s function.
Nevertheless, the performance of NOMA in THz networks is not well-understood
neither in single-carrier nor in multi-carrier set-up. Different from RF, the
THz transmissions are susceptible to unique challenges such as molecular
absorption noise111Molecular absorption noise causes signal loss as the
electromagnetic (EM) energy gets partially transformed into internal energy of
the molecules., molecular absorption at different frequencies leading to
serious path-loss peaks, and a sophisticated Beer’s Lambert law-based channel
model.

Recently, a handful of research works considered analyzing the performance of
NOMA in THz networks. In [6], the authors optimized the transmission powers to
maximize energy efficiency in a single carrier THz network. The authors in [7]
present an energy-efficient cooperative NOMA strategy for multi-user indoor
multi-input-single-output (MISO) THz network that assures the minimum required
rate for cell-edge users. Lately, preliminary research demonstrated the gains
of NOMA over orthogonal multiple access (OMA) in THz networks using computer
simulations [8]. The authors considered a single THz channel, no fading,
random user pairing, and the molecular absorption model was limited to 250 -
450 GHz. In [9], the authors optimized power allocation, user clustering, and
hybrid precoding to maximize energy efficiency in a THz-NOMA system. In [10],
the authors considered maximizing the network sum-rate while optimizing the
beamforming weights, sub-array selection, power allocation, and sub-band
assignment in a THz-NOMA network subject to the user’s quality of service
(QoS) requirements.

To the best of our knowledge, there is no research work that provides a
comprehensive framework for the OP analysis of users in a single-carrier and
multi-carrier THz-NOMA network and/or develops a low-complexity user pairing
scheme for the THz-NOMA network with a guaranteed gain over OMA for each
individual user. It is also noteworthy that analyzing the performance of
multi-carrier NOMA in RF networks is straight-forward as all sub-channels
experience identical channel statistics. However, in THz networks, each sub-
channel experiences a different molecular absorption indicated by its
molecular absorption coefficient. Thus the channel statistics on each sub-
channel are non-identical. Furthermore, the mathematical structure of the
channel propagation model based on Beer’s-Lambert law adds to the challenge.
Subsequently, characterizing OP expressions in a multi-carrier THz-NOMA
network is challenging. In the sequel, the main contributions of this paper
include:

$\bullet$ We develop a novel low-complexity user pairing scheme in a THz-NOMA
network. The proposed scheme ensures NOMA outperform OMA for each individual
user in the NOMA pair and adapts according to the molecular absorption.

$\bullet$ We characterize the exact outage expressions in a single-carrier and
multi-carrier THz-NOMA network considering Nakagami-$m$ fading to capture the
line-of-sight feature of THz transmissions, and molecular absorption noise.
Our expressions can be customized for various user-pairing schemes and are
applicable to the entire THz frequency range.

$\bullet$ For multi-carrier THz-NOMA network, we propose a moment-generating-
function (MGF)-based exact approach and derive simplified single integral
expressions to compute the OP, as opposed to Fox-H’s based expressions in [5].

$\bullet$ Our numerical results validate the accuracy of our derived
expressions and demonstrate the performance of the proposed user-pairing
scheme compared to the conventional random pairing and nearest-farthest
schemes.

##  II System model and Assumptions

We consider the downlink NOMA transmission of a single-antenna access point
(AP) operating at THz frequency. We consider two users222 Our framework can be
generalized for multiple users by making multiple two-user NOMA pairs in
orthogonal time/frequency resource blocks. However, the consideration of other
NOMA pairs would not effect the performance of users in a specific NOMA pair.
Two-user NOMA has been standardized as Multi-user Superposition Transmission
(MUST) in 3GPP [11]. in each NOMA cluster that are located at distances
$d_{1}$ and $d_{2}$ from the AP such that $d_{1}<d_{2}$. The users located at
$d_{1}$ and $d_{2}$ are referred to as user 1 ($U_{1}$) and user 2 ($U_{2}$),
respectively.

####  II-1 Channel Model

The line-of-sight (LoS) channel power between the AP and user $i$ is
formulated as follows [12]:

$|{h_{L}(d_{i})}|^{2}=\left(\frac{c}{4{\pi}fd_{i}}\right)^{2}\
e^{-k(f)d_{i}}=\zeta d_{i}^{-2}e^{-k(f)d_{i}},$ (1)

where the molecular absorption coefficient $k(f)$ is defined as:

${k(f)}\mathrm{=}\sum_{(i,g)}{\frac{p^{2}T_{\mathrm{sp}}q^{(i,g)}{N_{A}}S^{(i,g)}f\tanh{\left(\frac{hcf}{2k_{b}T}\right)}}{p_{0}VT^{2}f^{(i,g)}_{c}\tanh{\left(\frac{hcf^{(i,g)}_{c}}{2k_{b}T}\right)}}}F^{(i,g)}\left(f\right),$
(2)

where $p$ and $p_{0}$ indicate the ambient pressure of the transmission medium
and the reference pressure, respectively, $T$ is the temperature of the
transmission medium, $T_{\mathrm{sp}}$ denotes the temperature at standard
pressure, $q^{(i,g)}$ indicates the mixing ratio of the isotopologue $i$ of
gas $g$, $N_{A}$ refers to the Avogadro number, and $V$ is the gas constant.
The line intensity $S^{(i,g)}$ defines the strength of the absorption by a
specific type of molecules and is directly obtained from the HITRAN database
[13]. In addition, $f$ and $f^{(i,g)}_{c}$ denote the THz frequency and the
resonant frequency of gas $g$, respectively, $c$ is the speed of light, $h$ is
the Planck’s constant, and $k_{b}$ is the Boltzmann constant. For the
frequency band $f$, we consider the Van Vleck-Weisskopf asymmetric line shape
to evaluate:

$F^{(i,g)}(f)=\frac{100\>c\>\alpha^{(i,g)}f}{\pi\>f_{c}^{(i,g)}}\left(\frac{1}{Y^{2}+(\alpha^{(i,g)})^{2}}+\frac{1}{Z^{2}+(\alpha^{(i,g)})^{2}}\right),$

where $Y=f+f_{c}^{(i,g)}$ and $Z=f-f_{c}^{(i,g)}$, and the Lorentz half-width
is given as follows:

$\alpha^{(i,g)}=\left(\left(1-q^{(i,g)}\right)\alpha_{\mathrm{air}}^{(i,g)}+q^{(i,g)}\alpha_{0}^{(i,g)}\right)\left(\frac{p}{p_{0}}\right)\left(\frac{T_{0}}{T}\right)^{\gamma},$

where $T_{0}$ indicates the reference temperature, the parameters air half-
widths, $\alpha_{\mathrm{air}}^{(i,g)}$, self-broadened half-widths,
$\alpha_{0}^{(i,g)}$, and temperature broadening coefficient, $\gamma$, are
obtained directly from the HITRAN database [13]. The resonant frequency of gas
$g$ at reference pressure $p_{0}$ is determined as
$f_{c}^{(i,g)}=f_{{c_{0}}}^{(i,g)}+{\delta}^{{(i,g)}}{(\frac{p}{p_{0}})}$,
where $\delta^{(i,g)}$ is the linear pressure shift [14].

####  II-2 SINR - NOMA Model

The signal-to-interference-plus-noise ratio (SINR) of $U_{1}$ and $U_{2}$ with
perfect successive-interference cancellation are modeled, respectively, as
follows:

$\mathrm{SINR}_{\mathrm{1}}^{(\mathrm{noma})}=\frac{a_{1}G_{t}G_{r}P|{h_{L}(d_{1})}|^{2}\chi_{1}}{{N_{1}^{\mathrm{(noma)}}}},$
(3)
$\mathrm{SINR}_{\mathrm{2}}^{(\mathrm{noma})}=\frac{a_{2}G_{t}G_{r}P|{h_{L}(d_{2})}|^{2}\chi_{2}}{a_{1}G_{t}G_{r}P|{h_{L}(d_{2})}|^{2}\chi_{2}+{N_{2}^{\mathrm{(noma)}}}},$
(4)

where $\chi$ is Nakagami-$m$ fading channel, and $G_{t}$ and $G_{r}$ are the
directional antenna gains of AP and users, respectively. Beam alignment
strategies are assumed that align the main lobes of the users and the THz AP.
The noise at the receivers of $U_{1}$ and $U_{2}$ comprises of thermal noise
$N_{0}$ and molecular absorption noise as defined, respectively, below:

$\displaystyle N_{1}^{\mathrm{(noma)}}=N_{0}+a_{1}G_{t}G_{r}\zeta
Pd_{1}^{-2}(1-e^{-k(f)d_{1}})\chi_{1},$ (5)
$\displaystyle{N_{2}^{\mathrm{(noma)}}}$
$\displaystyle=N_{0}+(a_{1}+a_{2})G_{t}G_{r}\zeta
Pd_{2}^{-2}(1-e^{-k(f)d_{2}})\chi_{2},$ $\displaystyle=N_{0}+G_{t}G_{r}\zeta
Pd_{2}^{-2}(1-e^{-k(f)d_{2}})\chi_{2},$ (6)

where $a_{1}$ and $a_{2}$ depict the fraction of the AP transmit power
allocated for $U_{1}$ and $U_{2}$, respectively, such that $a_{1}+a_{2}=1$.
Also, $P$ denotes the total transmit power budget of the AP. The spectral
efficiency of $U_{1}$ and $U_{2}$ (in bps/Hz) is computed for a duration of
time $\hat{T}$ as follows:

$\mathrm{C}_{\mathrm{i}}^{(\mathrm{noma})}=\hat{T}\log_{2}(1+\mathrm{SINR}^{(\mathrm{noma})}_{i}),\quad\forall
i=\\{1,2\\}.$

####  II-3 SINR - OMA Model

The SINR for $U_{1}$ and $U_{2}$ in OMA, where each user receives its
transmission for a predefined duration of $\hat{T}/2$, is modeled as follows:

$\mathrm{SINR}_{{i}}^{(\mathrm{oma})}=\frac{G_{t}G_{r}P|{h_{L}(d_{i})}|^{2}\chi_{i}}{N^{(\mathrm{oma})}_{i}},\quad\forall
i=\\{1,2\\},$ (7)

where $N_{i}^{(\mathrm{oma})}=N_{0}+G_{t}G_{r}\zeta
Pd_{i}^{-2}\left(1-e^{-k(f)d_{i}}\right)\chi_{i}.$ Thus, the spectral
efficiency of $U_{1}$ and $U_{2}$ can be computed as

$\mathrm{C}_{{i}}^{(\mathrm{oma})}={\frac{\hat{T}}{2}\log_{2}}{\left(1+\mathrm{SINR}_{i}^{(\mathrm{oma})}\right)},\quad\forall
i=\\{1,2\\}.$

Without loss of generality, the duration $\hat{T}$ is taken as unity.

##  III Outage Analysis: single-Carrier THz-NOMA

In this section, we first describe the proposed user-pairing scheme along with
the two benchmark user-pairing schemes, describe the distance distributions of
the users, and present a framework to calculate the outage of users for the
proposed and benchmark user-pairing schemes.

For bench-marking purposes, we consider a random and nearest-farthest user
pairing schemes. In the random pairing scheme, we pick only two users randomly
with independent and identically distributed distances $r_{1}$ and $r_{2}$
from AP. The near and far user’s distance can thus be defined as
$d_{1}=\mathrm{min}(r_{1},r_{2})$ and $d_{2}=\mathrm{max}(r_{1},r_{2})$,
respectively. Therefore, the PDF and CDF of the distances of $U_{1}$ and
$U_{2}$ are given, respectively, as:

$f_{d_{1}}(d_{1})=2f_{r}\left(d_{1}\right)(1-F_{r}\left(d_{1}\right)),\>\>F_{d_{1}}(d_{1})=1-[1-F_{r}\left(d_{1}\right)]^{2},$
(8) $f_{d_{2}}(d_{2})=2f_{r}\left(d_{2}\right)F_{r}\left(d_{2}\right),\quad
F_{d_{2}}\left(d_{2}\right)=[F_{r}\left(d_{2}\right)]^{2},$ (9)

On the other hand, in the nearest-farthest scheme, we select two users out of
$N$ users with minimum and maximum distances from the AP. The near and far
user’s distance can thus be defined as
$d_{1}=\mathrm{min}(r_{1},r_{2},\cdots,r_{N})$ and
$d_{2}=\mathrm{max}(r_{1},r_{2},\cdots,r_{N})$, respectively. The PDF and CDF
of $d_{1}$ and $d_{2}$ can thus be given, respectively, as follows:

$f_{d_{1}}(d_{1})=N[1-F_{r}\left(d_{1}\right)]^{N-1}f_{r}\left(d_{1}\right),F_{d_{1}}(d_{1})=1-\left[1-F_{r}\left(d_{1}\right)\right]^{N},$
(10)
$f_{d_{2}}\left(d_{2}\right)=N\left[F_{r}\left(d_{2}\right)\right]^{N-1}f_{r}\left(d_{2}\right),F_{d_{2}}\left(d_{2}\right)=\left[F_{r}\left(d_{2}\right)\right]^{N},$
(11)

where the PDF and CDF of $r$ are given, respectively, as
$f_{r}\left(r\right)=\frac{2r}{R^{2}},\mathrm{and}\>F_{r}\left(r\right)=\frac{r^{2}}{R^{2}},$
since all users are uniformly distributed in a circular region of radius $R$.

Figure 1: An illustration of the proposed user pairing scheme.

###  III-A Proposed User Grouping Scheme

Different from the conventional approach where user-pairing is typically
performed to maximize the sum-rate performance of users in a NOMA pair, the
proposed user scheme selects only those users who can benefit from NOMA
transmissions compared to OMA transmission. The scheme enables fairness among
users by allowing each of them to join a NOMA pair if and only if they have a
benefit over OMA.  To ensure the gains of NOMA over OMA for each user, we
derive the necessary condition to ensure
${C}_{{i}}^{\mathrm{(noma)}}>{C}_{{i}}^{\mathrm{(oma)}},\>\forall i=\\{1,2\\}$
and $d_{1}<d_{2}$ as shown in the following.

######  Lemma 1 (Proposed User Pairing Scheme).

The proposed pairing scheme is a sufficient condition for NOMA to outperform
OMA, where the near user should be within the distance $R_{\mathrm{th1}}$ and
the far user should be located beyond the distance $R_{\mathrm{th2}}$, such
that $d_{1}<d_{2}$, where:

$R_{\mathrm{th1}}=\frac{1}{\mathrm{}{k(f)}}\
{\mathrm{ln}\left(\frac{1-a_{1}}{1-2a_{1}}\right)},$ (12)
$R_{\mathrm{th2}}=\frac{1}{\mathrm{}{k(f)}}\
{\mathrm{ln}\left(\frac{{a_{1}}^{2}}{1-2a_{1}}+1\right)}.$ (13)

Comparing (12) to (13), it is straight-foward to verify that
$R_{\mathrm{th1}}>R_{\mathrm{th2}}$ always, since $a_{1}+a_{2}=1$.

###### Proof.

Starting from the condition
${C}_{{i}}^{\mathrm{(noma)}}>{C}_{{i}}^{\mathrm{(oma)}},\>\forall i=\\{1,2\\}$
and taking $\hat{T}=1$ without loss of generality, we have:

$\log_{2}{\left(1+\mathrm{SINR}_{{i}}^{(\mathrm{noma})}\right)}>{0.5\log}_{2}{\left(1+\mathrm{SINR}_{{i}}^{(\mathrm{oma})}\right)},\>\forall
i=\\{1,2\\}.$ (14)

Now, we substitute (3) and (7) in (14) for near user, and (4) and (7) in (14)
for far user. After basic algebraic manipulations, we note that
$d_{1}<R_{\mathrm{th}1}$ and $d_{2}>R_{\mathrm{th}2}$ are necessary conditions
to guarantee the gains of NOMA over OMA for each individual user, where the
value of $R_{\mathrm{th}1}$ and $R_{\mathrm{th}2}$ can be written as in (12)
and (13), respectively. ∎

Evidently, as can also be seen in Fig. 1, four user-pairing cases are possible
to guarantee the gains of NOMA over OMA for each individual user. The three
regions can vary as a function of $k(f)$ and the allocated powers $a_{1}$ and
$a_{2}$, e.g., increasing $a_{1}$ will increase both $R_{\mathrm{th1}}$ and
$R_{\mathrm{th2}}$. However, increasing $k(f)$ will decrease both
$R_{\mathrm{th1}}$ and $R_{\mathrm{th2}}$. This is different from RF NOMA,
wherein there is only one case possible due to only one threshold distance (or
two regions).

######  Lemma 2.

The PDF and CDF of $d_{1}$ and $d_{2}$ in the proposed scheme are given,
respectively, as follows:

$f_{d_{1}}\left(d_{1}\right)=\frac{2d_{1}}{{R^{2}_{\mathrm{th1}}}},F_{d_{1}}\left(d_{1}\right)=\frac{{d^{2}_{1}}}{{R^{2}_{\mathrm{th1}}}},$
(15)
$f_{d_{2}}\left(d_{2}\right)=\frac{2d_{2}}{R^{2}-R_{\mathrm{th2}}^{2}},F_{d_{2}}\left(d_{2}\right)=\frac{d_{2}^{2}-R_{\mathrm{th2}}^{2}}{R^{2}-R_{\mathrm{th2}}^{2}}.$
(16)

###### Proof.

Since all users are uniformly distributed in a circular region of radius $R$,
thus the PDF and CDF of a user from the BS $r$ are given, respectively, as
$f_{r}\left(r\right)=\frac{2r}{R^{2}}\mathrm{and}\>F_{r}\left(r\right)=\frac{r^{2}}{R^{2}}$.
In the proposed scheme, the range of the near user is $d_{1}\in\
[0,R_{\mathrm{th}1}]$. Subsequently, $d_{1}$ follows the truncated
distribution of $r$ and its CDF can be calculated as follows:

$F_{d_{1}}\left(d_{1}\right)=\frac{F_{r}\left(d_{1}\right)-F_{r}\left(0\right)\
}{F_{r}\left(R_{\mathrm{th1}}\right)-F_{r}\left(0\right)}=\frac{\frac{{d_{1}}^{2}}{{R}^{2}}-\frac{({0)}^{2}}{{R}^{2}}}{\frac{R^{2}_{\mathrm{th1}}}{{R}^{2}}-\
\frac{({0)}^{2}}{{R}^{2}}\ }=\frac{{d^{2}_{1}}}{{R^{2}_{\mathrm{th1}}}},$

Now, we can calculate the PDF by taking the derivative of the CDF as shown in
Lemma 2. Similarly, the far user is located in the range $d_{2}\in\
[R_{\mathrm{th}2},R]$ and the CDF of $d_{2}$ is given as:

$F_{d_{2}}\left(d_{2}\right)=\frac{F_{r}\left(d_{2}\right)-F_{r}\left(R_{\mathrm{th2}}\right)\
}{F_{r}\left(R\right)-F_{r}\left(R_{\mathrm{th2}}\right)}=\frac{d_{2}^{2}-R_{\mathrm{th2}}^{2}}{R^{2}-R_{\mathrm{th2}}^{2}},$

Now, we can calculate the PDF of $d_{2}$ by taking the derivative of the CDF
as shown in Lemma 2. ∎

######  Corollary 1 (Enhanced Proposed Scheme).

The proposed scheme can be enhanced further in terms of spectral efficiency by
selecting the nearest user as $U_{1}$. In this case, $f_{d_{1}}(d_{1})$ and
$F_{d_{1}}(d_{1})$ can be given as in (10).

###### Proof.

To enhance the proposed scheme, we consider that the near user is the user
with minimum distance, thus, $f_{d_{1}}(d_{1})$ and $F_{d_{1}}(d_{1})$ can be
given as in (10). On the other hand, the far user is located beyond
$R_{\mathrm{th2}}$. So the range of the far user will be as $d_{2}\in\
[R_{\mathrm{th2}},R]$ , then in this case, $f_{d_{2}}(d_{2})$ and
$F_{d_{2}}(d_{2})$ can be given as in (16). ∎

###  III-B OP Analysis

The OP is defined as the probability that user $i$ does not achieve its
target-spectral efficiency $\tau_{i}$, i.e.,

${\mathcal{O}^{({\cdot})}_{i}}=\Pr{\left(\mathrm{C}^{({\cdot})}_{\mathrm{i}}\leq\
\mathrm{\tau}_{i}\right)},$ (17)

where ${({\cdot})}$ denotes NOMA or OMA transmission. Now we formally derive
the OP of users in NOMA and OMA separately.

####  III-B1 NOMA (Near User)

The OP of the near user in the downlink NOMA mode is given as:

$\displaystyle\mathcal{O}^{\mathrm{(noma)}}_{1}$
$\displaystyle{=}\mathrm{Pr}{\left(\frac{B_{1}d_{1}^{-2}{e^{{-k}\left(f\right)d_{1}}}{\chi}_{1}}{N_{0}+B_{1}d_{1}^{-2}\left(1-e^{-k\left(f\right)d_{1}}\right){\chi}_{1}}\leq\
{y}_{1}\right)}$
$\displaystyle{=}\int_{0}^{R}\frac{\gamma\left[m,\frac{{y_{1}}N_{0}d_{1}^{2}}{{\Theta
B}_{1}\left(e^{-k\left(f\right)d_{1}}\left(1+y_{1}\right)-{y_{1}}\right)}\right]}{\mathrm{\Gamma}(m)}f_{d_{1}}(d_{1})\;dd_{1},$
(18)

where $\gamma(\cdot)$ is the lower incomplete Gamma function, $\Gamma(\cdot)$
is the complete Gamma function, $m$ is the fading severity, $\Theta$ is the
fading power, $y_{1}=2^{\tau_{1}}-1$ and ${B_{1}{=a}_{1}G}_{t}G_{r}P\zeta$.
Since the thermal noise $N_{0}$ is negligible compared to the molecular
absorption noise in THz networks, (III-B1) is simplified as:

$\hat{\mathcal{O}}^{\mathrm{(noma)}}_{1}=1-\mathrm{Pr}{\left(d_{1}\leq\frac{\ln{(\frac{1+y_{1}}{{y}_{1}})}}{k\left(f\right)}\right)}=1-F_{d_{1}}\left(\frac{\ln{(\frac{1+y_{1}}{{y}_{1}})}}{k\left(f\right)}\right).$
(19)

Now substituting (8) in (19) or (III-B1) for random scheme, (10) in (19) or
(III-B1) for nearest-farthest scheme, and (15) in (19) or (III-B1) for
proposed scheme, gives us the outage expressions. The outage of $U_{1}$ in the
proposed scheme is given as:

$\hat{\mathcal{O}}^{\mathrm{(noma)}}_{1}=1-{{\left(\frac{\ln{(\frac{1+y_{1}}{{y}_{1}})}}{R_{\mathrm{th1}}\>k\left(f\right)}\right)^{2}}}.$
(20)

####  III-B2 OMA (Near User)

The OPs with and without noise can be given by replacing $B_{1}$ with
${A=G}_{t}G_{r}P\zeta$ and $y_{1}$ with ${{x_{1}}=2}^{2{\tau_{1}}}-1$ in
(III-B1), (19), and (20) for all schemes.

####  III-B3 NOMA (Far User)

The OP of the far user in the downlink NOMA mode is formulated as follows:

$\displaystyle\mathcal{O}^{\mathrm{(noma)}}_{2}=\mathrm{Pr}{\left(\frac{B_{2}|{h_{L}(d_{2})}|^{2}\chi_{2}}{B_{1}|{h_{L}(d_{2})}|^{2}\chi_{2}+{N_{2}^{(\mathrm{noma})}}}\leq
y_{2}\right)}$ $\displaystyle{=}\int_{0}^{R}\frac{\gamma\left[m,\
\frac{y_{2}N_{0}d_{2}^{2}}{{}_{\Theta}\left(e^{{-k}\left(f\right)d_{2}}({B_{2}}\
{+\ y_{2}B}_{2})-y_{2}\ ({B_{1}+B}_{2})\ \
\right)}\right]}{\Gamma(m)}f_{d_{2}}(d_{2})dd_{2},$ (21)

where, $B_{2}=(a_{1}+a_{2})G_{t}G_{r}P\zeta$. since $N_{0}$ is negligible
compared to molecular absorption noise, (III-B3) is simplified as:

$\begin{split}\hat{\mathcal{O}}^{\mathrm{(noma)}}_{2}=1-F_{d_{2}}\left(\frac{\ln{\left[\frac{a_{2}+a_{2}y_{2}}{y_{2}}\right]}}{k\left(f\right)}\right).\end{split}$
(22)

Now substituting (9) in (22) for random scheme, (11) in (22) for nearest-
farthest, and (16) in (22) for proposed scheme, gives us the outage
expressions. The outage of the far user in the proposed scheme is given as:

$\hat{\mathcal{O}}^{\mathrm{(noma)}}_{2}=1-\frac{\left(\frac{\ln{\left[\frac{a_{2}+a_{2}y_{2}}{y_{2}}\right]}}{k\left(f\right)}\right)^{2}-R_{\mathrm{th2}}^{2}}{R^{2}-R_{\mathrm{th2}}^{2}}.$
(23)

####  III-B4 OMA (Far User)

The OPs can be given by replacing $d_{1}$ with $d_{2}$, $B_{1}$ with
${A=G}_{t}G_{r}P\zeta$ and $y_{2}$ with ${{x_{2}}=2}^{2{\tau_{2}}}-1$ in
(III-B1), (19), and (20) for all schemes.

##  IV Outage Analysis: Multi-Carrier THz-NOMA

In this section, we present moment-generating function (MGF)-based approach to
derive tractable outage expressions of the users in a multi-carrier THz-NOMA
network. Both the near and far users will get $N$ subcarriers allocated. The
spectral efficiency of $U_{1}$ and $U_{2}$ is given, respectively, as:

$C_{1}^{\mathrm{(noma)}}=\sum_{n=1}^{N}\log_{2}{\left(1+\frac{{B_{1}{\chi}_{1,n}d_{1}^{-2}}{{e^{{-k}\left(f_{n}\right)d_{1}}}}}{N_{1,n}^{\mathrm{(noma)}}}\right)},$
(24)
$C_{2}^{\mathrm{(noma)}}=\sum_{n=1}^{N}\log_{2}{\left(1+\frac{B_{2}\chi_{2,n}d_{2}^{-2}e^{-k(f_{n})d_{2}}}{B_{1}\chi_{2,n}d_{2}^{-2}e^{-k(f_{n})d_{2}}+N_{2,n}^{(\mathrm{noma})}}\right)}.$
(25)

The OP can thus be formulated as follows:

$\displaystyle\mathcal{O}_{i}^{(\mathrm{noma})}$
$\displaystyle=\mathrm{Pr}\left(C_{i}^{\mathrm{(noma)}}<\tau_{i}\right)=\mathrm{Pr}{\left(\sum_{n=1}^{N}\mathrm{log_{2}}{\left({W}_{i,n}\right)}<\tau_{i}\right)},$
$\displaystyle=\mathrm{Pr}{\left(X_{i}=\sum_{n=1}^{N}\mathrm{\ln}{\left({W}_{i,n}\right)}<\tau_{i}\ln(2)\right)}.$
(26)

Note that the MGF is a useful tool to deal with the sum of random variables.
For instance, characterizing the PDF becomes analytically intractable for a
sum of random variables due to multiple convolutions required. On the other
hand, the MGF of a sum of random variables can be derived by deriving the
product of the MGF of all random variables. Therefore, to analyze the OP of
users in a multi-carrier THz-NOMA network, we resort to an MGF-based approach.
Our methodology is as follows: (i) we first derive the PDF of $W_{i,n}$
conditional on $d_{1}$ for near user and $d_{2}$ for far user, (ii) derive the
conditional MGF of $X_{i,n}=\mathrm{ln}\,W_{i,n}$, (iii) compute the
conditional cumulative MGF of $X_{i}=\sum_{n=1}^{N}\mathrm{ln}\,W_{i,n}$,
i.e., $M_{X_{i}|d_{i}}=\prod_{n=1}^{N}M_{X_{i,n}|d_{i}}$, and (iv) substitute
in the Gil-Pelaez inversion lemma to compute the OP as follows:

$\displaystyle\mathcal{O}_{i}^{(\mathrm{noma})}=\frac{1}{2}+\frac{1}{\pi}\int_{0}^{\infty}\frac{\mathrm{Im}[M_{X_{i}}(s)e^{j\omega\tau_{i}\ln(2)}]}{\omega}d\omega,$
(27)

where $s=j\omega$,
$M_{X_{i}}(s)=\mathbb{E}_{d_{i}}[\prod_{n=1}^{N},M_{{X_{i,n}|d_{i}}}]$, and

$\displaystyle M_{X_{i,n}|d_{i}}(s)$
$\displaystyle={\color[rgb]{0,0,0}E[e^{{-s}\ln{\left({W}_{i,n}\right)}}|d_{i}]=E[{{W_{i,n}}^{-s}}|d_{i}]},$
$\displaystyle=\int_{-\infty}^{\infty}{{W}_{i,n}}^{-s}{f}_{{W}_{i,n}}\left({W}_{i,n}\right){d}{{W}_{i,n}}.$
(28)

From (24) and (25), we have the following:

$\chi_{1,n}=g_{1}(W_{1,n})=\frac{N_{0}(W_{1,n}-1)}{W_{1,n}\left({e}^{-k\left(f_{n}\right)d_{1}}-1\right)B_{1}d_{1}^{-2}+B_{1}d_{1}^{-2}},$
${\chi}_{2,n}=g_{2}(W_{2,n})=\frac{N_{0}-N_{0}{W}_{2,n}}{{{B}_{2}d}_{2}^{-2}({W}_{2,n}-{a_{2}e}^{{-k}\left(f_{n}\right)d_{2}}{{\
W}_{2,n}}-1)}.$

Since $\chi$ is Gamma distributed, we apply random variable transformation to
get the PDF of $W_{i,n}$:

$\displaystyle
f_{W_{i,n}}\left(W_{i,n}\right)=\frac{B_{i}d_{i}^{-2}A_{i}}{N_{0}{e}^{-k\left(f_{n}\right)d_{i}}}\left.\left(\frac{\chi_{i,n}^{m-1}e^{\frac{-{\chi}_{i,n}}{\Theta}}}{\Theta^{m}\Gamma\left(m\right)}\right)\right|_{{{\chi}_{i,n}=g_{1}(W_{i,n})}},$
(29)

where
$A_{1}=\left({W}_{1,n}\left({e}^{-k\left(f_{n}\right)d_{1}}-1\right)+1\right)^{2}$
and
$A_{2}=\left({W}_{2,n}\left(a_{2}{e}^{-k\left(f_{n}\right)d_{2}}-1\right)+1\right)^{2}/a_{2}$.
Finally, we can obtain $M_{X_{1,n}}(s)$ for $U_{1}$ and $M_{X_{2,n}}(s)$ for
$U_{2}$ by substituting (29) for $i=1$ in (IV) and $i=2$ in (IV),
respectively.

######  Lemma 3 (OP in Multi-carrier Network).

The outage can be computed of user $i$ in a multi-carrier THz-NOMA network
with the proposed scheme can be derived as follows:

$\displaystyle\mathcal{O}_{i}^{(\mathrm{noma})}$
$\displaystyle=\mathrm{Pr}{\left(X_{i}=\sum_{n=1}^{N}\mathrm{\ln}{\left({W}_{i,n}\right)}<\tau_{i}\ln(2)\right)}$
$\displaystyle=\int_{0}^{R}\mathbb{U}(X_{i}-\tau_{i}\mathrm{ln}(2))f_{d_{i}}(d_{i})dd_{i}.$
(30)

Now substituting (8) in (3) for random scheme and (10) in (3) for nearest-
farthest scheme gives the respective OPs.

###### Proof.

From (24) and (25), we can write

$\displaystyle{W}_{1,n}=1+\frac{{B_{1}{\chi}_{1,n}d_{1}^{-2}}{{e^{{-k}\left(f_{n}\right)d_{1}}}}}{N_{1,n}^{\mathrm{(noma)}}}$
(31)
$\displaystyle{W}_{2,n}=1+\frac{B_{2}\chi_{2,n}d_{2}^{-2}e^{-k(f_{n})d_{2}}}{B_{1}\chi_{2,n}d_{2}^{-2}e^{-k(f_{n})d_{2}}+N_{2,n}^{(\mathrm{noma})}}$
(32)

For enhanced tractability, we ignore the thermal noise and obtain the
following simplified results after algebraic manipulations for near and far
users, respectively.

${W}_{1,n}=\frac{1}{1-e^{-k\left(f_{n}\right)d_{1}}}\quad\mathrm{or}\quad
d_{1}=\frac{1}{k\left(f_{n}\right)}{\ln{\left[\frac{{W}_{1,n}}{{W}_{1,n}-1}\right]}},$
${W}_{2,n}=\frac{1}{1-{a_{2}e}^{-k\left(f_{n}\right)d_{2}}}\quad\mathrm{or}\quad
d_{2}=\frac{1}{k\left(f_{n}\right)}{\ln{\left[\frac{a_{2}{W}_{2,n}}{{W}_{2,n}-1}\right]}}.$

From (IV) and using the fact that $X_{i}$ is a constant conditional on $d_{i}$
and the CDF of a constant is a unit-step function, we got to (3) as in Lemma
3. ∎

From (12), we note that the $R_{\mathrm{th}1}$ and $R_{\mathrm{th}2}$ will
vary for each carrier, since they depend on the frequency-dependent absorption
coefficient $k(f_{n})$ which is computed using (2). Thus, each carrier should
be allocated to a user selected from a different region based on the proposed
scheme. Therefore, we extended our scheme for multi-carrier networks by
choosing a user from a region that is valid for all subcarriers.

######  Corollary 2.

We consider
$R_{\mathrm{th}1}^{\mathrm{min}}=\mathrm{min}(R_{\mathrm{th_{n}}}),\forall
n\in\\{1,2,\cdots,N\\}$ and
$R_{\mathrm{th}2}^{\mathrm{max}}=\mathrm{max}(R_{\mathrm{th_{n}}}),\forall
n\in\\{1,2,\cdots,N\\}$ for near user and far user, respectively.
Subsequently, after applying $R_{\mathrm{th1}}^{\mathrm{min}}$ and
$R_{\mathrm{th2}}^{\mathrm{max}}$, the PDF of the distance of near user and
far user can be given, respectively, as follows:

$\displaystyle
f_{d_{1}}\left(d_{1}\right)=\frac{2d_{1}}{{(R_{\mathrm{th1}}^{\mathrm{min}})^{2}}},f_{d_{2}}\left(d_{2}\right)=\frac{2d_{2}}{R^{2}-{(R_{\mathrm{th2}}^{\mathrm{max}})^{2}}}.$
(33)

###### Proof.

The near user will need to be located inside $R_{\mathrm{th}1}$. However, each
subcarrier will observe a different molecular absorption coefficient resulting
in a different threshold distance $R_{\mathrm{th}1}$ at each subcarrier.
Choosing the smallest threshold distance will not violate the threshold
requirement of all other subcarriers, thus we consider
$R_{\mathrm{th}1}^{\mathrm{min}}=\mathrm{min}(R_{\mathrm{th_{n}}}),\forall
n\in\\{1,2,\cdots,N\\}$ for near user. The range of $d_{1}$ becomes $d_{1}\in\
[0,R_{\mathrm{th}1}^{\mathrm{min}}]$, thus replacing $R_{\mathrm{th}1}$ with
$R_{\mathrm{th}1}^{\mathrm{min}}$ with in Lemma 2, we can obtain the result in
Corollary 2. Similarly, the far user will need to be located outside
$R_{\mathrm{th}2}$. Therefore, choosing the maximum threshold distance will
not violate the threshold requirement of all other subcarriers, i.e.,
$R_{\mathrm{th}2}^{\mathrm{max}}=\mathrm{max}(R_{\mathrm{th_{n}}}),\forall
n\in\\{1,2,\cdots,N\\}$ for far user. The range of $d_{2}$ becomes $d_{2}\in\
[0,R_{\mathrm{th}1}^{\mathrm{min}}]$, thus replacing $R_{\mathrm{th}2}$ with
$R_{\mathrm{th}2}^{\mathrm{max}}$ in Lemma 2, we can obtain the result in
Corollary 2. ∎

Now substituting (33) in (3) for the proposed scheme gives the outage
expressions of users in our proposed scheme.

Table I: Simulation Parameters for Calculating $k(f)$ in (2)  Symbol |  Value |  Symbol |  Value
---|---|---|---
$p_{0},p$ |  1 atm, 1 atm |  $q^{(i,g)}$ |  0.05 [%]
$T_{0},T$ |  296 K, 396 K |  $k_{b}$ |  1.3806$\times 10^{-23}$ J/K
$f_{{c_{0}}}^{(i,g)}$ |  276 Hz |  $T_{\mathrm{sp}}$ |  273.15 K
$\gamma$ |  0.83 |  $N_{A}$ |  6.0221 $\times 10^{23}$
$S^{(i,g)}$ |  2.66-25Hz-m${}^{2}/$mol |  h |  6.6262$\times 10^{-34}$ J s
$\alpha_{0}^{(i,g)},\alpha_{\mathrm{air}}^{(i,g)}$ |  0.916Hz, 0.1117Hz |  c |  2.9979 $\times 10^{8}$ m/s
${\delta}^{{(i,g)}}$ |  0.0251 Hz |  V |  8.2051$\times 10^{-5}$m3atm/K/mol

##  V Numerical Results and Discussions

In this section, we compare the performance of $U_{1}$ and $U_{2}$ in a
single-carrier and multi-carrier THz-NOMA network, considering a variety of
user-pairing schemes (i) random scheme, (ii) proposed scheme, (iii) enhanced
scheme, and (iv) near-far user pairing schemes.

Unless stated otherwise, the parameters are listed herein. We consider 300
users are uniformly distributed in a circular disc of radius 60 m. The antenna
gains $G_{t}$ and $G_{r}$ are set as 20 dB. The AP transmit power is 1W and
the power allocation coefficients $a_{1}=0.33$ and $a_{2}=1-a_{1}$. Note that
our framework is general for any arbitrary value of $a_{1}$ in the range
$0\leq a_{1}<0.5$. Nakagami-$m$ fading parameter is set as 2 and $\Omega=1$.
In multi-carrier NOMA, we consider six subcarriers, where each subcarrier has
the same transmission bandwidth, with the frequencies [0.85, 0.9, 0.95, 1.0,
1.05, 1.1] THz and their respective $k(f)=$[0.0357, 0.04, 0.0446, 0.0494,
0.0545, 0.0598] m-1 computed using (2) considering water vapour molecules. We
list the numerical values of parameters in Table I which are taken from [15].

Figure 2: Outage performance of near and far users as a function of the
molecular absorption coefficient in THz spectrum considering a single-carrier
network, $R$ = 60 m, $\tau_{1}$ = 3 bps/Hz, and $\tau_{2}$ = 0.5 bps/Hz.

Fig. 2 depicts the OP of $U_{1}$ and $U_{2}$ for different values of $k(f)$
and validates (III-B1) and (III-B3) through Monte-Carlo simulations. It is
shown that the values obtained through derived expressions (shown in circles),
exactly match those obtained through simulations (shown by lines). The
enhanced scheme (Corollary 1) significantly outperforms the random and
nearest-farthest schemes. It is interesting to note that our scheme adapts
user selection based on the molecular absorption coefficient $k(f)$. That is,
$R_{\mathrm{th1}}$ and $R_{\mathrm{th2}}$ reduce with the increase in $k(f)$
as can be seen from (12) and (13), respectively. Thus, closeby users are
selected to combat the effect of increased molecular absorption. Finally, it
can be seen that the OP increases with the increase in frequency and $k(f)$
due to increased molecular absorption; therefore, a lower THz frequency with
lower $k(f)$ is preferred.

Figure 3: Outage performance of near and far users as a function of the power
allocation coefficient of near user considering a single-carrier network,
$k(f)=0.03$ $R$ = 60 m, $\tau_{1}$ = 3 bps/Hz and $\tau_{2}$ = 0.5 bps/Hz.

Fig. 3 demonstrates the OP of $U_{1}$ and $U_{2}$ as a function of the power
allocation coefficient of near user $(a_{1})$ and highlights the gain of NOMA
over OMA. With the increase in $a_{1}$, the OP at $U_{2}$ increases due to the
increased interference from $U_{1}$ and increasing values of
$R_{\mathrm{th2}}$ as can be seen in (13). On the other hand, with the
increase in $a_{1}$, the OP at $U_{1}$ increases due to the increasing values
of $R_{\mathrm{th1}}$ as can be seen in (12). Our scheme adapts the user
pairing according to the change in $a_{1}$. The gain of NOMA over OMA improves
for low values of $a_{1}$ and is significant for the enhanced and proposed
schemes. Besides, we can observe that NOMA’s performance gets close to OMA’s
when $a_{1}$ is close to 0.5.

Figure 4: Outage performance of near and far users as a function of the number
of subcarriers in THz spectrum, $\tau_{1}$=8 bps/Hz and $\tau_{2}$=0.5 bps/Hz.

Fig. 4 shows the OP of $U_{1}$ and $U_{2}$ as a function of the number of
channels for different schemes in a multi-carrier network. The OP decreases
significantly when the number of channels allocated to the user increases due
to aggregate spectral efficiency. We note that the increase in subchannels
benefits the near user more as the outage decreases exponentially; whereas,
the outage decreases linearly for the far user due to its channel conditions.
Our analytical results match well with the simulations, and it is evident that
the enhanced scheme outperforms benchmark schemes for both users.

##  VI Conclusion

We analyzed the performance of users considering single-carrier and multi-
carrier set-up in both THz-NOMA and THz-OMA network. The derived OP
expressions are general to capture the entire range of THz spectrum,
Nakagami-$m$ fading and molecular absorption noise. We have developed an
adaptive pairing scheme for THz-NOMA network where user selection adapts
according to molecular absorption and the gains of NOMA are guaranteed for
each individual user. The framework can be extended to analyze the diversity
order or optimize network parameters, e.g., transmit power allocations.

## References

  * [1] M. Rasti, S. K. Taskou, H. Tabassum, and E. Hossain, “Evolution toward 6G wireless networks: A resource management perspective,” _arXiv preprint arXiv:2108.06527_ , 2021.
  * [2] B. Xia et al, “Outage performance analysis for the advanced SIC receiver in wireless NOMA systems,” _IEEE Trans. on Veh. Technol._ , vol. 67, no. 7, pp. 6711–6715, 2018.
  * [3] H. Tabassum et al, “Uplink vs. downlink NOMA in cellular networks: Challenges and research directions,” in _2017 IEEE 85th Veh Technol. Conf. (VTC Spring)_ , 2017, pp. 1–7.
  * [4] Z. Ding, Z. Yang, P. Fan, and H. V. Poor, “On the performance of non-orthogonal multiple access in 5G systems with randomly deployed users,” _IEEE Signal Process. Lett._ , vol. 21, no. 12, pp. 1501–1505, 2014\.
  * [5] S. Li et al., “Outage probability analysis for the multi-carrier NOMA downlink relying on statistical CSI,” _IEEE Trans. on Commun._ , vol. 68, no. 6, pp. 3572–3587, 2020.
  * [6] S. R. Sabuj, A. M. S. Khan, and M. Hamamura, “Application of non-orthogonal multiple access for machine type communication in sub-Terahertz band,” _Comput. Netw._ , vol. 182, p. 107508, 2020.
  * [7] O. Maraqa et al, “Energy-efficient coverage enhancement of indoor THz-MISO systems: An FD-NOMA approach,” in _2021 IEEE 32nd Annu. Int. Symp. Pers., Indoor, Mobile Radio Commun (PIMRC)_ , 2021, pp. 483–489.
  * [8] O. Ülgen, et al, “Non-orthogonal multiple access for Terahertz communication networks,” in _11th IEEE Annu. Ubiquitous Comput. Electron. Mobile Commun. Conf. (UEMCON)_ , 2020, pp. 0737–0742.
  * [9] H. Zhang, H. Zhang, W. Liu, K. Long, J. Dong, and V. C. M. Leung, “Energy efficient user clustering, hybrid precoding and power optimization in Terahertz mimo-noma systems,” _IEEE Journal on Sel. Areas in Commun._ , vol. 38, no. 9, pp. 2074–2085, 2020.
  * [10] X. Zhang, C. Han, and X. Wang, “Joint beamforming-power-bandwidth allocation in Terahertz NOMA networks,” in _16th Annu. IEEE Intl. Conf. on Sens., Commun, and Netw. (SECON)_ , 2019, pp. 1–9.
  * [11] 3GPP, “Study on downlink multiuser superposition transmission (must) for lte (release 13),” in _TR36.859, Dec. 2015_.
  * [12] J. Sayehvand and H. Tabassum, “Interference and coverage analysis in coexisting RF and dense Terahertz wireless networks,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 10, pp. 1738–1742, 2020.
  * [13] “The HITRAN2016 molecular spectroscopic database,” _J. of Quantitative Spectroscopy and Radiative Transfer_ , vol. 203, pp. 3 – 69, 2017, HITRAN2016 Special Issue. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0022407317301073
  * [14] J. M. Jornet and I. F. Akyildiz, “Channel modeling and capacity analysis for electromagnetic wireless nanonetworks in the Terahertz band,” _IEEE Trans. on Wireless Commun._ , vol. 10, no. 10, pp. 3211–3221, 2011.
  * [15] A. Faisal et al., “Ultramassive MIMO systems at Terahertz bands: Prospects and challenges,” _IEEE Veh Technol. Mag._ , vol. 15, no. 4, pp. 33–42, 2020.