# Finite temperature effects on the structural stability of Si-doped HfO2
using first-principles calculations

Y. Harashima<EMAIL_ADDRESS>Division of Materials Science, Nara
Institute of Science and Technology, Ikoma 630-0192, Japan  Center for
Computational Sciences, University of Tsukuba, Tsukuba 305-8577, Japan  H.
Koga  Department of Computational Science and Technology, Research
Organization for Information Science and Technology, Tokyo 105-0013, Japan  Z.
Ni  S-Technology Development Center, Tokyo Electron Technology Solutions Ltd.
T. Yonehara  Department of Computational Science and Technology, Research
Organization for Information Science and Technology, Tokyo 105-0013, Japan  M.
Katouda  Department of Computational Science and Technology, Research
Organization for Information Science and Technology, Tokyo 105-0013, Japan  A.
Notake  Advanced Data Planning Department, Tokyo Electron Ltd.  H. Matsui
S-Technology Development Center, Tokyo Electron Technology Solutions Ltd.  T.
Moriya  Advanced Data Planning Department, Tokyo Electron Ltd.  M. K. Si
Center for Computational Sciences, University of Tsukuba, Tsukuba 305-8577,
Japan  R. Hasunuma  Faculty of Pure and Applied Science, University of
Tsukuba, Tsukuba 305-8573, Japan  A. Uedono  Faculty of Pure and Applied
Science, University of Tsukuba, Tsukuba 305-8573, Japan  Y. Shigeta  Center
for Computational Sciences, University of Tsukuba, Tsukuba 305-8577, Japan

###### Abstract

The structural stabilities of the monoclinic and tetragonal phases of Si-doped
HfO2 at finite temperatures were analyzed using a computational scheme to
assess the effects of impurity doping. The finite temperature effects
considered in this work represented lattice vibration and impurity
configuration effects. The results show that 6% Si doping stabilizes the
tetragonal phase at room temperature, although a higher concentration of Si is
required to stabilize the tetragonal phase at zero temperature. These data
indicate that lattice vibration and impurity configuration effects are
important factors determining structural stability at finite temperatures.

##  I Introduction

Hafnium-based materials, typically HfO2, have found numerous applications in
the semiconductor industry as high-k dielectric compounds Wilk  _et al._
(2001) and are also promising ferroelectric and antiferroelectric materials
Park  _et al._ (2015). These compounds were first employed in the 200s, such
as during the 45 nm technology node in 2007 Mistry  _et al._ (2007) and in
certain earlier applications in dynamic random-access memory Kim and Popovici
(2018). The application of high-k materials allowed the continued scaling of
various devices by suppressing the gate leakage associated with greater
physical thicknesses exhibited by conventional SiO2 dielectric materials.
Since then, HfO2 and analogues obtained by doping have become crucial
components of many advanced device structures, ranging from gate-all-around
transistors Ernst  _et al._ (2006) to post-silicon 2D channels such as
transition metal dichalcogenides Li  _et al._ (2019); Xia  _et al._ (2017).

The exceptional electronic properties of HfO2 appear primarily in the
tetragonal and orthorhombic phases. Notably, the dielectric constant,
$\kappa$, of tetragonal HfO2 is $\sim 30$ and so is approximately twice that
of stable monoclinic HfO2 Tomida  _et al._ (2006). However, the tetragonal and
orthorhombic phases are not stable at room temperature Wang  _et al._ (1992);
Tomida  _et al._ (2006). Consequently, to allow the application of these
materials in various devices, it is important to develop techniques to
stabilize the tetragonal phase.

Techniques for controlling the phase have been developed through doping of
impurities at different concentration. As an example, a prior study produced a
Si-doped HfO2 thin film having a ferroelectric orthogonal structure Böscke
_et al._ (2011a, b). Additionally, recent experimental work fabricated
heterogeneous Hf1-xZrxO2 almost entirely composed of a non-monoclinic phase
using physical vapor deposition with a composition depth profile suggested by
a machine learning model (that is, using Bayesian optimization) Ni and Matsui
(2022). Fischer et al. studied doping with the group IV elements Si, C, Ge,
Sn, Ti and Ce into bulk HfO2 Fischer and Kersch (2008a, b) and found that Si
and Ge stabilize the tetragonal structure. Doping with Si, P, Ge, Al, Y, Ti,
Zr, Gd and Sc was also investigated in HfO2 by Lee et al. Lee  _et al._
(2008), who identified the microscopic mechanisms by which these dopants
stabilize the tetragonal structure on the basis of ionic radii. In other work
Barabash (2017), possible metastable phases of HfO2 were analyzed using first-
principles calculations. As well, a systematic screening of dopants capable of
inducing ferroelectric characteristics in HfO2 has been performed Batra  _et
al._ (2017).

When assessing the viability of synthesizing a particular phase, a theoretical
analysis of the structural stability at finite temperatures is required. In
addition, the total energies obtained from first-principles calculations based
on density functional theory indicate structural stability at zero
temperature. To investigate the stabilization imparted by doping with
impurities at finite temperatures, it is necessary to consider both lattice
vibration effects and impurity configuration effects. In a previous study
Fischer and Kersch (2008b), calculations were performed that incorporated the
lattice vibrational free energy based on considering phonon dispersion.
Typically, the impurity configuration is assumed to be homogeneous, such that
the impurity atoms are periodically aligned within the host lattice. Recently
the authors performed a systematic search for elements capable of stabilizing
HfO2 and ZrO2 using first-principles calculations Harashima  _et al._ (2022).
In this prior work, several impurity atoms were included in a supercell and
the total energies for all possible impurity configurations in the supercell
were estimated. In this manner, the most stable configuration at zero
temperature was ascertained. The results indicated that Si substituted at 6 %
was aligned along linearly and thus not homogeneously distributed throughout
the host matrix, suggesting that the dopant configuration cannot be neglected
when calculating energy values.

The present study investigated the finite temperature phase stability of Si-
doped HfO2, based on a computational scheme that simultaneously considered
lattice vibration and impurity configuration effects. On the basis of the
results, the transformation between the tetragonal and monoclinic phases of
HfO2 in association with Si doping is discussed. The data indicate that 6% Si
doping stabilizes the tetragonal structure at room temperature while the
monoclinic structure is the most stable at zero temperature.

##  II Methods

Figure 1: Schematic of the model used to approximate impurity configurations,
showing a configuration $\Lambda$ is shown. The left-hand side of the figure
presents two types of supercells, $q_{1}$ (red) and $q_{2}$ (green), having
different impurity configurations. The impurities are indicated by blue
spheres. The overall system is modelded as a combination of supercells and the
contribution from each cell is approximated by that of a piece of a system in
which the pieces are aligned periodically. Two periodically aligned systems
($q_{1}$-type and $q_{2}$-type) are considered, as shown in the right-hand
side of the figure, and the electronic states and vibrational modes for these
two systems are calculated. Here, $N_{q_{1}}$ and $N_{q_{2}}$ are the numbers
of the two types of cells in the entire system.

This work determined the free energies of impurity substituted systems
considering impurity configurations and lattice vibrations. The calculations
are based on a system having impurity configuration $\Lambda$ and lattice
vibrational degrees of freedom $\Gamma$. The state of $\Gamma$ is dependent on
$\Lambda$ but this relationship is not explicitly written. In this process,
the partition function $Z$ is expressed as the sum over all possible of
$\Lambda$, written as

$\displaystyle Z$ $\displaystyle=\sum_{\Lambda}\sum_{\Gamma}\exp(-\beta
E_{\Lambda}(\Gamma))$ (1) $\displaystyle=\sum_{\Lambda}Z_{\Lambda},$ (2)
$\displaystyle Z_{\Lambda}$ $\displaystyle\equiv\sum_{\Gamma}\exp(-\beta
E_{\Lambda}(\Gamma)).$ (3)

Here, $E_{\Lambda}(\Gamma)$ denotes the electronic and vibrational energy of
the configuration $\Lambda$ and $Z_{\Lambda}$ is the partition function for
this configuration. It is not feasible to calculate the exact value of $Z$ or
even of $Z_{\Lambda}$ because the system size is very large (or even infinite
at the thermodynamic limit). Therefore, the system is separated into
supercells and the energy of the overall system is approximated by summing
over these supercells. This process can be expressed in terms of $Z_{\Lambda}$
as

$\displaystyle Z_{\Lambda}$ $\displaystyle=\prod_{I}z_{I}.$ (4)

In this case, the impurities are distributed in each supercell $I$. The
energies for a single supercell $I$ can be approximated as the energies for a
system of periodically aligned impurities (see Fig. 1). Considering the free
energy $f_{I}$ of a periodically aligned supercell, we can write

$\displaystyle Z_{\Lambda}$
$\displaystyle\approx\prod_{I}\sum_{\gamma_{I}}\exp(-\beta e_{I}(\gamma_{I}))$
(5) $\displaystyle=\prod_{I}\exp(-\beta f_{I})$ (6)
$\displaystyle=\exp\left(-\beta\sum_{I}f_{I}\right).$ (7)

Here, $e_{I}$ denotes the energy of supercell $I$. These energies include
contributions from vibrational modes $\gamma_{I}$ in the approximated periodic
system, such that the free energy $f_{I}$ in cell $I$ is

$f_{I}\equiv-\dfrac{1}{\beta}\ln\sum_{\gamma_{I}}\exp(-\beta
e_{I}(\gamma_{I})).$ (8)

Note that this $f_{I}$ is the value per supercell. Because the supercell size
is finite, there is only a finite number of impurity configuration patterns
that need to be considered. The partition function for the overall system is
therefore written in terms of the free energy $f_{q}$ of each pattern $q$ as

$\displaystyle Z$
$\displaystyle\approx\sum_{\Lambda}\exp\left(-\beta\sum_{I}f_{I}\right)$ (9)
$\displaystyle=\sum_{\Lambda}\exp\left(-\beta N\sum_{q}n_{q}f_{q}\right)$ (10)
$\displaystyle=\sum_{\\{n_{q}\\}}W_{\\{n_{q}\\}}\exp\left(-\beta
N\sum_{q}n_{q}f_{q}\right),$ (11)

where $n_{q}$ is the proportion of cells having pattern $q$ in the entire
system and $N$ is the number of supercells which, as noted, will be huge or
even infinite at the thermodynamic limit. The energy value is approximated by
summing over the patterns and the summation over $\Lambda$ can be obtained by
summation over $\\{n_{q}\\}$, as when going from Eq. (10) to (11). Using
Stirling’s approximation, the weight $W_{\\{n_{q}\\}}$ can be estimated as

$W_{\\{n_{q}\\}}=\dfrac{N!}{\prod_{q}N_{q}!}\approx\exp\left(-N\sum_{q}n_{q}\ln
n_{q}\right).$ (12)

$N_{q}$ denotes $Nn_{q}$. Substituting Eq. (12) to Eq. (11), we obtain

$\displaystyle Z\approx\sum_{\\{n_{q}\\}}\exp\left(-\beta
N\sum_{q}n_{q}\left(f_{q}+\dfrac{1}{\beta}\ln n_{q}\right)\right).$ (13)

At the thermodynamic limit ($N\rightarrow\infty$), the most probable set
$\\{n_{q}\\}$ is realized and this set satisfies the conditions

$\displaystyle\dfrac{\partial}{\partial
n_{q}}\sum_{q}n_{q}\left(f_{q}+\dfrac{1}{\beta}\ln n_{q}\right)=0,$ (14)
$\displaystyle\sum_{q}n_{q}=1.$ (15)

By solving this equation, we obtain

$n_{q}=\dfrac{\exp(-\beta f_{q})}{\sum_{q^{\prime}}\exp(-\beta
f_{q^{\prime}})}.$ (16)

The dominant contribution in Eq. (13) is given by the set of Eq. (16). Taking
the logarithm, the free energy per supercell can be written as

$f\equiv\dfrac{F}{N}=-\dfrac{1}{\beta N}\ln
Z\approx-\dfrac{1}{\beta}\ln\left[\sum_{q}\exp(-\beta f_{q})\right].$ (17)

This formula includes the contributions of configurational entropy, by summing
over $q$, and of vibrational free energy via $f_{q}$. The derivation of this
formula is similar to that of the canonical ensemble but uses $f_{q}$, instead
of energy.

##  III Results and Discussion

The relationships derived in the preceding section were employed to analyze
the phase stability of Si-doped HfO2. HfO2 can exist as a monoclinic phase
having a low dielectric constant or a tetragonal phase having a high
dielectric constant. The former is stable at room temperature and the
transformation from the monoclinic to the tetragonal phases occurs at 2052 K
Wang  _et al._ (2006). Here, the free energy values for both phases at finite
temperatures were compared.

The $f_{q}$ term was obtained using first-principles calculations, employing
the Vienna Ab initio Simulation Package (VASP) based on the projector
augmented-wave method Blöchl (1994); Kresse and Furthmüller (1996); Kresse and
Joubert (1999). The generalized gradient approximation was also employed for
the exchange-correlation energy functional Perdew  _et al._ (1996). The
crystal structures were numerically optimized based on preserving the lattice
symmetry of the pristine compound and applying $a=b$ and
$\alpha=\beta=\gamma=\pi/2$ for the tetragonal phase and $a=b$ and
$\alpha=\gamma=\pi/2$ for the monoclinic phase. Each dimension of the
supercell was twice that of the conventional cell ($2\times 2\times 2$ of the
conventional cell), which contains 96 atoms in total. The impurity
configuration patterns could be classified into several groups based on
spatial symmetry. As an example, 2 Hf substitutions in the tetragonal $2\times
2\times 2$ supercell can produce ${}_{32}C_{2}=496$ impurity configurations
that can be classified into nine symmetrically equivalent groups. The
supercell code Okhotnikov  _et al._ (2016) was used to identify these groups.
The $496$ configurations are separated into 128 patterns $\times$ 2 groups, 64
patterns $\times$ 2 groups, 32 patterns $\times$ 2 groups, and 16 patterns
$\times$ 3 groups. 2 Hf substitution into the monoclinic supercell produces 23
groups which contains $32\times 8$ and $16\times 15$ groups. In the case of
practical calculations (that is, for 2 Hf substitutions in the tetragonal
structure), Eq. (17) is estimated as

$\displaystyle f$
$\displaystyle=-\dfrac{1}{\beta}\ln\left[\sum_{q^{\prime}=1}^{9}M_{q^{\prime}}\cdot\exp(-\beta
f_{q^{\prime}})\right],$ (18) $\displaystyle M_{q^{\prime}}$
$\displaystyle=\\{128,128,64,\cdots\\},$ (19)

such that nine $f_{q^{\prime}}$ values can be calculated using first-
principles calculations as explained below. By comparing $f$ for the
tetragonal and monoclinic phases, the effect of temperature on phase stability
and the structural transformation temperature could be obtained.

In this study, the vibrational free energy of each impurity configuration,
$f_{q^{\prime}}$, was estimated using the Debye-Grüneisen model Moruzzi  _et
al._ (1988) based on the equations

$\displaystyle F_{q^{\prime}}(V,T)=E_{\mathrm{el}}(V)$ $\displaystyle\quad-
k_{\mathrm{B}}T\left[D(\dfrac{\Theta_{\mathrm{D}}}{T})-3\ln\left(1-\exp\left(-\dfrac{\Theta_{\mathrm{D}}}{T}\right)\right)\right]+\dfrac{9}{8}k_{\mathrm{B}}\Theta_{\mathrm{D}},$
(20)
$\displaystyle\Theta_{\mathrm{D}}=(6\pi^{2})^{\frac{1}{3}}\dfrac{\hbar}{k_{\mathrm{B}}}v^{\frac{1}{6}}\left(\dfrac{B}{m}\right)^{\frac{1}{2}},$
(21)

and

$\displaystyle D(x)=\dfrac{3}{x^{3}}\int_{0}^{x}dx^{\prime}\dfrac{x^{\prime
3}}{e^{x^{\prime}}-1},$ (22)

where $\Theta_{\mathrm{D}}$ and $D(x)$ are the Debye temperature and Debye
function, respectively. $v$ and $m$ are averaged volume and mass per atom in
supercell, and $B$ is a bulk modulus. The free energy term $F_{q^{\prime}}$ in
Eq. (20) is calculated by first ascertaining the bulk modulus. This can be
accomplished by considering the Murnaghan equation of states

$\displaystyle
B(p)=B_{0}+B_{0}^{\prime}p=B_{0}\left(\dfrac{V}{V_{0}}\right)^{-B_{0}^{\prime}},$
(23) $\displaystyle
E_{\mathrm{el}}(V)=E_{0}-\dfrac{B_{0}}{B_{0}^{\prime}}V_{0}\left[\dfrac{1}{1-B_{0}^{\prime}}\left(\dfrac{V}{V_{0}}\right)^{1-B_{0}^{\prime}}-\dfrac{V}{V_{0}}-\dfrac{B_{0}^{\prime}}{1-B_{0}^{\prime}}\right].$
(24)

The total energies were calculated for volumes expanded and contracted from
the stable structure by several percent and $B_{0}$ and $B_{0}^{\prime}$ were
obtained by regression of these total energies. The stable lattice structure
for a given temperature was determined as the minimum value of Eq. (20) and
the free energy values were estimated for an arbitrary temperature. It should
be noted that the present approach involved a systematic error. Specifically,
the structural transformation temperatures for pristine HfO2 and ZrO2 were
estimated using this process as benchmarks and values of 1000 K and 750 K were
obtained. These values are scaled by a factor of 2 with respect to the
experimental data, and so the calculated transformation temperatures for the
doped systems should possible be scaled by a factor of 2 as well. The accuracy
of these calculations could possibly be improved by replacing the process used
to assess vibrational modes with first-principles phonon calculations Togo and
Tanaka (2015); Tadano  _et al._ (2014); Kuwabara  _et al._ (2005); Luo  _et
al._ (2009) More sophisticated calculations such as these will be examined in
future work.

Figure 2: The transformation temperature as a function of the extent of Si
doping in Hf1-xSixO2.

The transformation temperature was estimated based on the crossing point of
the free energies of the tetragonal and monoclinic phases. Figure 2 plots the
structural transformation temperature $T_{\mathrm{s}}$ of HfO2 versus the
extent of Si doping and demonstrates that the correlation between the two is
not linear. The $T_{\mathrm{s}}$ value for 3% Si doping is almost the same as
that of the pristine compound but rapidly decreases as the extent of
substitution increases from 3% to 6%. The decreasing trend against Si
concentration is consistent with the results of the previous studies. Fischer
and Kersch (2008b) The $T_{\mathrm{s}}$ value determined for 6% substitution
was 130 K. Although, as noted, this temperature may be underestimated and
could possibly be scaled by a factor of 2, the value would still be lower than
room temperature. This finding indicates that a 6% Si-doped HfO2 dielectric
layer in a typical device will exist in the tetragonal phase at a normal
operating temperature.

Figure 3: Effect of temperature on the proportion $n_{q}$ of the monoclinic
phase (including the multiplicity $M_{q^{\prime}}$). The most probable
configurations are shown in the upper panels at the various temperatures. The
gold, red and blue spheres denote Hf, O and Si, respectively.

The temperature dependence of the impurity configurations was also examined.
The proportion of each impurity configuration in a system defined as in Eq.
(16) can be estimated from the $f_{q}$ calculated for each temperature. In the
tetragonal phase, a single configuration (and symmetrically equivalent
configurations) was dominant up to $\simeq 1500$ K. In contrast, in the case
of the monoclinic phase, three impurity configurations appeared depending on
temperature. Figure 3 summarizes the effect of temperature on $n_{q}$ for the
monoclinic phase. It is evident that the most stable configuration varied with
temperature, indicating that the local structure of the impurities affected
phase stability. Even though the calculated transformation temperature was
lower than the temperatures at which several configurations were mixed,
$\simeq 250$ K or $\simeq 800$ K, the temperatures at which such compounds are
synthesized will be similar to or much higher than these values. Hence, it is
important to take into account the appearance of several configurations during
the synthetic process.

##  IV Conclusion

This work analyzed the structural stability of Si-doped HfO2 using the scheme
that estimated the free energy and considered impurity configurations and
lattice vibrations. The present process is evidently a useful means of
assessing the phase stability of doped polymorphic compounds. The structural
transformation from the monoclinic phase to the tetragonal phase was also
investigated based on analyzing the free energy difference. The transformation
temperature values indicate that a tetragonal structure with 6% Si doping will
be stable at room temperature while a monoclinic structure is more stable at
zero temperature. In the case of a device operating at room temperature,
doping with 6% Si was found to be sufficient to stabilize the tetragonal HfO2
phase.

In general, the physical properties of a compound will be highly correlated
with its crystal structure, as demonstrated by the HfO2 used in this work as a
model case. The proposed method in this paper is expected to be applicable to
a variety of scenarios. As an example, in permanent magnets, certain compounds
exhibit excellent magnetic properties but are unstable and so it is important
to stabilize these substances Miyake  _et al._ (2014); Harashima  _et al._
(2016, 2018, 2020). In particular, stabilization at the high temperatures at
which these compounds are synthesized has recently attracted much attention
Xing  _et al._ (2021). We expect that the method demonstrated herein will be
effective in such cases and could be helpful in the search for structural
stabilizers.

###### Acknowledgements.

This work was supported by JSPS KAKENHI grant number JP22K03449. The
computation was partly conducted using the facilities of the Supercomputer
Center at the Institute for Solid State Physics in the University of Tokyo,
the Multidisciplinary Cooperative Research Program at the Center for
Computational Sciences in University of Tsukuba, Subsystem B of the ITO system
in Kyushu University, and Research Center for Computational Science, Okazaki,
Japan (Project: 22-IMS-C104).

## References

  * Wilk  _et al._ (2001) G. D. Wilk, R. M. Wallace, and J. M. Anthony, Journal of Applied Physics 89, 5243 (2001), https://doi.org/10.1063/1.1361065 .
  * Park  _et al._ (2015) M. H. Park, Y. H. Lee, H. J. Kim, Y. J. Kim, T. Moon, K. D. Kim, J. Müller, A. Kersch, U. Schroeder, T. Mikolajick, and C. S. Hwang, Advanced Materials 27, 1811 (2015), https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.201404531 .
  * Mistry  _et al._ (2007) K. Mistry, C. Allen, C. Auth, B. Beattie, D. Bergstrom, M. Bost, M. Brazier, M. Buehler, A. Cappellani, R. Chau, C.-H. Choi, G. Ding, K. Fischer, T. Ghani, R. Grover, W. Han, D. Hanken, M. Hattendorf, J. He, J. Hicks, R. Huessner, D. Ingerly, P. Jain, R. James, L. Jong, S. Joshi, C. Kenyon, K. Kuhn, K. Lee, H. Liu, J. Maiz, B. McIntyre, P. Moon, J. Neirynck, S. Pae, C. Parker, D. Parsons, C. Prasad, L. Pipes, M. Prince, P. Ranade, T. Reynolds, J. Sandford, L. Shifren, J. Sebastian, J. Seiple, D. Simon, S. Sivakumar, P. Smith, C. Thomas, T. Troeger, P. Vandervoorn, S. Williams, and K. Zawadzki, in  _2007 IEEE International Electron Devices Meeting_ (2007) pp. 247–250.
  * Kim and Popovici (2018) S. K. Kim and M. Popovici, MRS Bulletin 43, 334 (2018).
  * Ernst  _et al._ (2006) T. Ernst, C. Dupre, C. Isheden, E. Bernard, R. Ritzenthaler, V. Maffini-Alvaro, J.-C. Barbe, F. De Crecy, A. Toffoli, C. Vizioz, S. Borel, F. Andrieu, V. Delaye, D. Lafond, G. Rabille, J.-M. Hartmann, M. Rivoire, B. Guillaumot, A. Suhm, P. Rivallin, O. Faynot, G. Ghibaudo, and S. Deleonibus, in  _2006 International Electron Devices Meeting_ (2006) pp. 1–4.
  * Li  _et al._ (2019) W. Li, J. Zhou, S. Cai, Z. Yu, J. Zhang, N. Fang, T. Li, Y. Wu, T. Chen, X. Xie, K. Ma, Haibo Yan, N. Dai, X. Wu, H. Zhao, Z. Wang, D. He, L. Pan, Y. Shi, P. Wang, W. Chen, K. Nagashio, X. Duan, and X. Wang, Nature Electronics 2, 563 (2019).
  * Xia  _et al._ (2017) P. Xia, X. Feng, R. J. Ng, S. Wang, D. Chi, C. Li, Z. He, X. Liu, and K.-W. Ang, Scientific Reports 7, 40669 (2017).
  * Tomida  _et al._ (2006) K. Tomida, K. Kita, and A. Toriumi, Applied Physics Letters 89, 142902 (2006), https://doi.org/10.1063/1.2355471 .
  * Wang  _et al._ (1992) J. Wang, H. P. Li, and R. Stevens, Journal of Materials Science 27, 5397 (1992).
  * Böscke  _et al._ (2011a) T. S. Böscke, J. Müller, D. Bräuhaus, U. Schröder, and U. Böttger, Applied Physics Letters 99, 102903 (2011a), https://doi.org/10.1063/1.3634052 .
  * Böscke  _et al._ (2011b) T. S. Böscke, S. Teichert, D. Bräuhaus, J. Müller, U. Schröder, U. Böttger, and T. Mikolajick, Applied Physics Letters 99, 112904 (2011b), https://doi.org/10.1063/1.3636434 .
  * Ni and Matsui (2022) Z. Ni and H. Matsui, Japanese Journal of Applied Physics 61, SH1009 (2022).
  * Fischer and Kersch (2008a) D. Fischer and A. Kersch, Applied Physics Letters 92, 012908 (2008a), https://aip.scitation.org/doi/pdf/10.1063/1.2828696 .
  * Fischer and Kersch (2008b) D. Fischer and A. Kersch, Journal of Applied Physics 104, 084104 (2008b), https://doi.org/10.1063/1.2999352 .
  * Lee  _et al._ (2008) C.-K. Lee, E. Cho, H.-S. Lee, C. S. Hwang, and S. Han, Phys. Rev. B 78, 012102 (2008).
  * Barabash (2017) S. V. Barabash, J. Comput. Electron. 16, 1227 (2017).
  * Batra  _et al._ (2017) R. Batra, T. D. Huan, G. A. J. Rossetti, and R. Ramprasad, Chemistry of Materials 29, 9102 (2017), https://doi.org/10.1021/acs.chemmater.7b02835 .
  * Harashima  _et al._ (2022) Y. Harashima, H. Koga, Z. Ni, T. Yonehara, M. Katouda, A. Notake, H. Matsui, T. Moriya, M. K. Si, R. Hasunuma, A. Uedono, and Y. Shigeta, in  _2022 International Symposium on Semiconductor Manufacturing (ISSM)_ (2022) pp. 1–3.
  * Wang  _et al._ (2006) C. Wang, M. Zinkevich, and F. Aldinger, Journal of the American Ceramic Society 89, 3751 (2006), https://ceramics.onlinelibrary.wiley.com/doi/pdf/10.1111/j.1551-2916.2006.01286.x .
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Okhotnikov  _et al._ (2016) K. Okhotnikov, T. Charpentier, and S. Cadars, J. Cheminform. 8, 17 (2016).
  * Moruzzi  _et al._ (1988) V. L. Moruzzi, J. F. Janak, and K. Schwarz, Phys. Rev. B 37, 790 (1988).
  * Togo and Tanaka (2015) A. Togo and I. Tanaka, Scripta Materialia 108, 1 (2015).
  * Tadano  _et al._ (2014) T. Tadano, Y. Gohda, and S. Tsuneyuki, Journal of Physics: Condensed Matter 26, 225402 (2014).
  * Kuwabara  _et al._ (2005) A. Kuwabara, T. Tohei, T. Yamamoto, and I. Tanaka, Phys. Rev. B 71, 064301 (2005).
  * Luo  _et al._ (2009) X. Luo, W. Zhou, S. V. Ushakov, A. Navrotsky, and A. A. Demkov, Phys. Rev. B 80, 134119 (2009).
  * Miyake  _et al._ (2014) T. Miyake, K. Terakura, Y. Harashima, H. Kino, and S. Ishibashi, J. Phys. Soc. Jpn. 83, 043702 (2014).
  * Harashima  _et al._ (2016) Y. Harashima, K. Terakura, H. Kino, S. Ishibashi, and T. Miyake, J. Appl. Phys. 120, 203904 (2016).
  * Harashima  _et al._ (2018) Y. Harashima, T. Fukazawa, H. Kino, and T. Miyake, J. Appl. Phys. 124, 163902 (2018).
  * Harashima  _et al._ (2020) Y. Harashima, T. Fukazawa, and T. Miyake, Scripta Materialia 179, 12 (2020).
  * Xing  _et al._ (2021) G. Xing, T. Ishikawa, Y. Miura, T. Miyake, and T. Tadano, Journal of Alloys and Compounds 874, 159754 (2021).


# A General Model for Pointing Error of High Frequency Directional Antennas

Mohammad. T. Dabiri, Mazen Hasna, Senior Member, IEEE, Nizar Zorba, Senior
Member, IEEE, Tamer Khattab, Senior Member, IEEE, and Khalid Qaraqe, Senior
Member, IEEE Mohammad Taghi Dabiri, Mazen Hasna, Nizar Zorba, and Tamer
Khattab are with the Department of Electrical Engineering, Qatar University,
Doha, Qatar. (E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>tkhattab@ieee.org).Khalid A. Qaraqe is with the Department of Electrical and
Computer Engineering, Texas A$\&$M University at Qatar, Doha 23874, Qatar
(E-mail: khalid.qaraqe@qatar.tamu.edu).

###### Abstract

This paper focuses on providing an analytical framework for the quantification
and evaluation of the pointing error for a general case at high-frequency
millimeter wave (mmWave) and terahertz (THz) communication links. For this
aim, we first derive the the probability density function (PDF) and cumulative
distribution functions (CDF) of the pointing error between an unstable
transmitter (Tx) and receiver (Rx), that have different antenna patterns and
for which the vibrations are not similar in the Yaw and Pitch directions. The
special case where the Tx and Rx are both equipped with uniform linear array
antenna is also investigated. In addition, using $\alpha-\mu$ distribution,
which is a valid model for small-scale fading of mmWave/THz links, the end-to-
end PDF and CDF of the considered channel is derived for all the considered
cases. Finally, by employing Monte-Carlo simulations, the accuracy of the
analytical expressions is verified and the performance of the system is
studied.

###### Index Terms:

Antenna pattern, antenna misalignment, backhaul links, pointing errors,
mmWave, THz systems.

##  I Introduction

In the quest to increase network capacity, high-frequency millimeter wave
(mmWave) and terahertz (THz) communication links are widely considered to be
one of the next frontiers for future wireless systems, because of the wide
swaths of unused and unexplored spectrum. While the mmWave frequencies are
already considered in 5G, the vendors and operators have started to look at
the usage of the enormous bandwidth of higher THz frequencies [1] and[2]. To
mitigate the negative effects of the high path-loss at the mmWave/THz bands,
the small wavelength enables the realization of a compact form of highly
directive antenna arrays, which also allow a better spatial reuse [3] and [4].
This is an important feature for aerial and even space nodes such as low- and
high-altitude unmanned aerial vehicles (UAVs), and low earth orbit (LEO)
satellites, which are in great needs of large data rates. Therefore, one of
the attractive scenarios is to use the large available bandwidth at mmWave/THz
frequencies in order to provide extra data rate for point-to-point aerospace
communications.

Figure 1: A nominal illustration of the next generation of aerial networks
includes a large number of high-frequency and ultra-high data rates between
aerial nodes that require the use of directional antennas.

One of the main challenges of the tackled scenario is mobility, making the
high-directional mmWave/THz antennas suffer from transceivers antenna
misalignment, known as pointing error for aerial/mobile nodes. Therefore,
along with other research topics in mmWave/THz bands, studying the effect of
pointing errors on system performance is a necessity for establishing a high
reliable communication link, which has recently been the subject of several
studies [5, 6, 7, 8, 9, 10, 11, 12, 13]. In summary, in [5, 6, 7, 8, 9, 10,
11, 12], the well-known pointing error model provided in [14] is used, which,
although suitable for optical communication systems and a special case of
mmWave/THz systems, as we will show in this work, cannot be directly used for
the typical mmWave/THz communication systems.

More recently, in [13] a new pointing error model is provided for mmWave/THz
communication links in the presence of an unstable transmitter (Tx) and
receiver (Rx) nodes as a function of real antennas’ pattern. However, the
results of [13] are obtained for a very special case where the intensity of
the vibrations of both Tx and Rx is the same in the direction of Yaw and
Pitch, and also the Tx and Rx nodes have the same antenna pattern with the
same stability. However, as shown in Fig. 1, the Tx and Rx nodes are not
necessarily of the same type, not necessarily at the same altitude and
atmospheric conditions, and as a result, their instability is not the same. As
an example, for a communication link between a small multi-rotor UAV and a
balloon, it is clear that their instability is different, and more
importantly, because the balloon is capable of carrying more payload, it is
possible to use a larger antenna for the balloon. In addition, the intensity
of the angular vibrations of the nodes is not necessarily the same in the
directions of Yaw and Pitch. For example, although for some types of hovering
UAVs the intensity of vibrations in both directions of Yaw and Pitch is almost
the same, for fixed-wing UAVs or airplanes, the intensity of vibrations in the
direction of Pitch is greater than the intensity of vibrations in the
direction of Yaw. Therefore, the model provided in [13] is used only in the
special scenario where the Tx and Rx are perfectly symmetric, standing as a
major drawback if wanted to use in practical systems. In order to design and
analyze the mmWave/THz connections between the moving nodes, it seems
necessary to provide a general and tractable model for pointing errors based
on the severity of Tx and Rx vibrations, which is the main contribution of
this work.

In this work, we consider a general case between two unstable Tx and Rx nodes
that are completely asymmetric, meaning that they have different antenna
patterns and the vibrations of the nodes are not necessarily the same in the
directions of Yaw and Pitch. For the considered general scenario, we derive
the probability density function (PDF) and the cumulative distribution
functions (CDF) of pointing error by tacking into account the 3D real antenna
pattern. We show that the obtained closed-form models are completely different
from the pointing error models of an optical link, and therefore, the pointing
error models of the optical channels cannot be directly used to model the
pointing error of mmWave/THz links. Then, for a special case that the standard
deviation (SD) of vibrations are approximately equal in the directions of Yaw
and Pitch, we derive a more tractable model for the PDF and CDF of pointing
error. Notice that, for the special case where the Tx and Rx nodes are
completely symmetric (meaning that the SD of vibrations of Tx and Rx nodes
along with the antenna pattern are the same), the pointing error models
provided in this work reduce to the simple pointing error model provided in
[13]. In some situations, such as communications between Low-altitude UAVs,
the use of a uniform linear array (ULA) antenna in the direction perpendicular
to the ground can have several advantages compared to the planar array
antenna. Therefore, for the case where the Tx and Rx are equipped with ULA
antennas, we also derive the PDF and CDF of pointing error. In addition, using
$\alpha-\mu$ distribution, which is a valid model for small-scale fading of
higher frequences (i.e., mmWave and THz links), the end-to-end PDF and CDF of
the considered channel are derived for all the mentioned cases. Finally, by
employing Monte-Carlo simulations, the accuracy of the analytical expressions
is verified and the performance of the system is studied.

The remainder of the paper is structured as follows. In Section II, the
channel of a point-to-point communication link between two unstable Tx and Rx
is characterized. The analytical derivations for pointing error are provided
in Section III along with several numerical and simulation results to verify
the accuracy of the derivations. In Section IV, the analytical derivations for
the end-to-end channel are provided for different scenarios. At the end of
Section IV, by providing several simulation and analytical results, the
performance of the considered system is studied. Finally, conclusions are
drawn in Section V.

Figure 2:  Graphical illustration of pointing errors of two directional
antennas which are mounted on an unstable Tx and Rx. The relationship between
the axes of the Cartesian coordinate system and the roll, Yaw, and Pitch
directions are depicted in the considered system model.

##  II The System Model

Consider an optical receiver where a lens is used to focus the incoming
received signal onto the small detector or feed antenna located at the focal
point. In this case, if the center of received Gaussian beam deviates by
$d_{v}$ from the center of the receiver lens, the fraction of the collected
power by the the circular lens with radius $a$ (called the pointing error) is
calculated as [14, 15]

$\displaystyle
h_{\text{p}}=\int_{-a}^{a}\int_{-\sqrt{a^{2}-x^{2}}}^{\sqrt{a^{2}-x^{2}}}\frac{2}{\pi
w_{z}^{2}}\exp\left(-2\frac{(x-d_{v})^{2}+y^{2}}{w_{z}^{2}}\right)\text{d}x\text{d}y,$
(1)

where $w_{z}$ is the beamwidth of the received Gaussian beam. Using Eq. (1)
and after a series of manipulations, the distribution of $h_{p}$ is derived in
[14, Eq. (11)]. As mentioned earlier, in the pointing error studies of THz
communications, it is typical to refer to [14, Eq. (11)] in order to model the
pointing error. Although this model is suitable for FSO systems and/or a
specific class of high frequency communication links that use a circular lens
in the receiver, in this work, we show that this model is not an accurate
model for Lens-less mmWave/THz communication systems.

Consider Fig. 2, where we tackle a point-to-point high frequency communication
system between two unstable nodes, and where Tx and Rx are equipped with high
directional antennas. The channel model can be written as [16, 17]

$\displaystyle
h=h_{L}h_{a}\sqrt{G_{t}(\theta_{t},\phi_{t})G_{r}(\theta_{r},\phi_{r})},$ (2)

where $G_{q}(\theta_{q},\phi_{q})$ is the antenna radiation pattern defined in
Eq. (II), $h_{a}$ is the small scale fading, $h_{L}=h_{Lf}h_{Lm}$ is the
channel path loss, $h_{Lf}=\left(\frac{\lambda}{4\pi Z}\right)^{2}$ is the
free-space path loss, and $h_{Lm}=e^{-\frac{\mathcal{K}(f)}{2}Z}$ denotes the
molecular absorption loss where $\mathcal{K}(f)$ is the frequency dependent
absorption coefficient. The experimental results show that water vapor
dominates the molecular absorption loss at high frequencies [18]. In
mmWave/THz channels, the $\alpha-\mu$ distribution is one of the most widely
used models employed for the distribution function of random variables (RV)
$h_{a}$ which is given as [19]

$\displaystyle
f_{h_{a}}(h_{a})=\frac{\alpha\mu^{\mu}}{\hat{h_{a}}^{\alpha\mu}\Gamma(\mu)}h_{a}^{\alpha\mu-1}\exp\left(-\mu\frac{h_{a}^{\alpha}}{\hat{h_{a}}^{\alpha}}\right),$
(3)

where $\Gamma(\cdot)$ is the Gamma function, $\alpha>0$ is a fading parameter,
$\mu$ is the normalized variance of the fading channel envelope, and
$\hat{h_{a}}$ is the $\alpha$-root mean value of the fading channel envelope.

In Eq. (2), the parameter $G_{q}(\theta_{q},\phi_{q})$ is the antenna
radiation pattern in the directions of $\theta_{q}$ and $\phi_{q}$ where the
subscript $q\in\\{t,r\\}$ determines the Tx and Rx nodes. Here, we consider a
standard uniform $N\times N$ array antenna. By taking into account the effect
of all elements [17], the array radiation gain in the direction of
$\theta_{q}$ and $\phi_{q}$ will be
$G_{q}(\theta_{q},\phi_{q})=G_{0}(N_{q})G^{\prime}_{q}(\theta_{q},\phi_{q})$,
where

$\displaystyle G^{\prime}_{q}(N_{q},\theta_{q},\phi_{q})$
$\displaystyle=\left(\frac{\sin\left(\frac{N_{q}(kd_{x}\sin(\theta_{q})\cos(\phi_{q}))}{2}\right)}{N_{q}\sin\left(\frac{kd_{x}\sin(\theta_{q})\cos(\phi_{q})}{2}\right)}\right.$
$\displaystyle\times\left.\frac{\sin\left(\frac{N_{q}(kd_{y}\sin(\theta_{q})\sin(\phi_{q}))}{2}\right)}{N_{q}\sin\left(\frac{kd_{y}\sin(\theta_{q})\sin(\phi_{q})}{2}\right)}\right)^{2},$
(4)

and $d_{x}=d_{y}=\frac{\lambda}{2}$ are the spacing between the elements along
the $x$ and $y$ axes, respectively, $k=\frac{2\pi}{\lambda}$ denotes the wave
number, $\lambda=\frac{c}{f_{c}}$ represents the wavelength, $f_{c}$ denotes
the carrier frequency and $c$ is the speed of light. Also, in order to
guarantee that the total radiated power of antennas with different $N$ are the
same, the coefficient $G_{0}$ is defined as

$\displaystyle
G_{0}(N_{q})=\frac{4\pi}{\int_{0}^{\pi}\int_{0}^{2\pi}G^{\prime}_{q}(\theta_{q},\phi_{q})\sin(\theta_{q})d\theta_{q}d\phi_{q}}.$
(5)

Based on Eq. (II), the maximum value of the antenna gain is equal to
$G_{0}(N)$, which is obtained when $\theta_{q}=0$.

As shown in Fig. 2, we assume that the Tx and Rx are located on axis $z$ and
at a distance $Z$ from each other, and both the Tx and Rx try to place the
main lobe of the antenna pattern on the $z$ axis. The use of mmWave/THz high-
gain antennas makes them more sensitive to antenna misalignment or pointing
errors, especially for mobile or aerospace communications. The angular
fluctuations of each mobile node are modeled in three directions; namely: Yaw,
roll, and Pitch. Due to the symmetry in the main-lobe of the antenna pattern
in the direction of roll (angle $\phi_{q}$), we can neglect the effect of
orientation deviations in the roll direction on the pointing error. Hence, the
orientation fluctuations in the Yaw and Pitch directions cause a pointing
error. As shown in Fig. 2, orientation deviations in the directions of Yaw and
Pitch are equivalent to the orientation deviations in $x-z$ and $y-z$ planes,
respectively. Let $\theta_{tx}\sim\mathcal{N}(0,\sigma_{\theta_{tx}})$ and
$\theta_{ty}\sim\mathcal{N}(0,\sigma_{\theta_{ty}})$ denote the orientation
fluctuations of Tx in $x-z$ and $y-z$ planes, respectively, while
$\theta_{rx}\sim\mathcal{N}(0,\sigma_{\theta_{rx}})$ and
$\theta_{ry}\sim\mathcal{N}(0,\sigma_{\theta_{ry}})$ denote the orientation
fluctuations of Rx in $x-z$ and $y-z$ planes, respectively. Therefore, in our
model, the parameters $\theta_{q}$ can be defined as functions of RV
$\theta_{qx}$ and $\theta_{qy}$ as follows:

$\displaystyle\theta_{q}$
$\displaystyle=\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{qx})+\tan^{2}(\theta_{qy})}\right).$
(6)

##  III Pointing Error Modeling

In this section, we first derive the distribution function of pointing error
for the general case and then provide more tractable and precise analytical
expressions for some specific cases. Moreover, the accuracy of the derived
analytical expressions is verified by Monte-Carlo simulations.

Let’s rewrite Eq. (2) as:

$\displaystyle h=h_{L}h_{a}h_{p},$ (7)

where $h_{p}=G_{0}(N_{t},N_{r})h^{\prime}_{p}$ is the pointing error
coefficient, $G_{0}(N_{t},N_{r})=\sqrt{G_{0}(N_{t})G_{0}(N_{r})}$ is the
maximum antenna gain, and $0<h^{\prime}_{p}<1$ is the normalized pointing
error

$\displaystyle
h^{\prime}_{p}=h^{\prime}_{pt}h^{\prime}_{pr}=\sqrt{G^{\prime}_{t}(N_{t},\theta_{tx},\theta_{ty})}\sqrt{G^{\prime}_{r}(N_{t},\theta_{rx},\theta_{ry})}.$
(8)

From Eq. (8), $h^{\prime}_{p}$ is a function of $N_{t}$ and $N_{r}$ as well as
four RVs $\theta_{tx},\theta_{ty}$, $\theta_{rx}$, and $\theta_{ry}$.

Let us define parameters
$\beta_{qw}=\frac{\sigma_{\theta_{qw}}^{2}}{w_{B_{q}}}$ where the subscript
$w\in\\{x,y\\}$ determines the vibrations of Yaw (or $x-z$ plane) and Pitch
(or $y-z$ plane), and the subscript $q\in\\{t,r\\}$ determines the Tx and Rx
nodes. Based on this definition, the general scenarios refer to the condition
that $\sigma_{\theta_{qw}}$s and $w_{Bq}$s (or $N_{q}$s) have different
values.

Theorem 1. Under the general condition when the parameters
$\sigma_{\theta_{qw}}$s and $w_{Bq}$s have different values, the PDF of
pointing error is obtained as:

$\displaystyle
f_{h_{p}}(h_{p})=C_{g}\frac{h_{p}^{1/\beta_{1}-1}}{G_{0}^{1/\beta_{1}}}\sum_{k=0}^{K}\frac{\Delta_{k}\left(-\ln\left(\frac{h_{p}}{G_{0}}\right)\right)^{k+1}}{\Gamma(k+2)\beta_{1}^{k+2}}$
(9)

while its corresponding CDF results in

$\displaystyle
F_{h_{p}}(h_{p})=1-\Theta_{0}+C_{g}\left(\frac{h_{p}}{G_{0}}\right)^{1/\beta_{1}}$
$\displaystyle~{}~{}~{}\times\sum_{k=0}^{K}\sum_{j=0}^{k+1}\frac{\Delta_{k}}{\Gamma(k+2-j)\beta_{1}^{k+1-j}}\left(-\ln\left(\frac{h_{p}}{G_{0}}\right)\right)^{k+1-j}$
(10)

where

$\displaystyle\left\\{\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{rl}&C_{g}=\prod_{i=1}^{4}\sqrt{\beta_{1}/\beta_{i}},\\\
&\Delta_{k}=\frac{1}{k}\sum_{i=1}^{k}i\gamma_{i}\Delta_{k+1-i}~{}~{}\text{for}~{}~{}k=1,...,K,\\\
&\Delta_{0}=1,~{}~{}~{}\gamma_{k}=\sum_{i=1}^{4}\frac{\left(1-{\beta_{1}}/{\beta_{i}}\right)^{k}}{2k}.\end{array}\right.$
(14)

Also, $\bar{\beta}=[\beta_{1},...,\beta_{4}]$ where $\beta_{1}$ is the least
value of $\beta_{q_{w}}$s, $\beta_{4}$ is the maximum value of
$\beta_{q_{w}}$s, and
$\beta_{qw}=\frac{\sigma^{2}_{\theta_{q_{w}}}}{w^{2}_{Bq}}$.

###### Proof:

Please refer to Appendix A. ∎

In Fig. 3, by employing Monte-Carlo simulations, we compare the accuracy of
the proposed analytical derivations in Theorem 1 for the three different
states of transceiver instability. As mentioned, the vibrations of Tx are
characterized by $\sigma_{\theta_{tx}}$ and $\sigma_{\theta_{ty}}$, while the
vibrations of Rx are characterized by $\sigma_{\theta_{rx}}$ and
$\sigma_{\theta_{ry}}$. In the simulations, we use the standard array antenna
pattern introduced in Eqs. (II) and (5). To do this, we first generate
$5\times 10^{7}$ independent RVs $\theta_{tx}$, $\theta_{ty}$, $\theta_{rx}$,
and $\theta_{ry}$, and then, for each $5\times 10^{7}$ independent run, using
Eqs. (II), (5), and (6), we generate $5\times 10^{7}$ independent RVs of
$h_{p}$. Next, using the generated instantaneous coefficients $h_{p}$s, the
PDF and CDF of pointing error are obtained. As can be seen, the Monte-Carlo
simulations verify the accuracy of the analytical expressions provided in
Theorem 1, especially for low $\sigma_{\theta_{qw}}$ values. The higher the
value of $\sigma_{\theta_{qw}}$, the probability of deviation from the main
lobe to the side lobes increases, resulting in a gap between the simulation
results and the analytical results. For example, for $N_{t}=25$ and
$N_{r}=30$, although the analytical results for
$\sigma_{\theta_{tx}}=0.6^{o},\sigma_{\theta_{ty}}=0.5^{o},\sigma_{\theta_{rx}}=0.7^{o}$,
and $\sigma_{\theta_{ry}}=0.4^{o}$ completely match the simulations, by
increasing $\sigma_{\theta_{qw}}$s to
$\sigma_{\theta_{tx}}=1.2^{o},\sigma_{\theta_{ty}}=1.3^{o},\sigma_{\theta_{rx}}=1.0^{o}$,
and $\sigma_{\theta_{ry}}=1.1^{o}$, it is observed that the analytical results
deviate slightly from the simulations.

(a)

(b)

Figure 3: Comparison of the accuracy of the provided analytical results in
Theorem 1 with Monte-Carlo simulations for three different states of the
intensity of transceiver vibrations (a) for the PDF, and (b) for the CDF.

Theorem 2. When $\beta_{tx}=\beta_{ty}=\beta_{t}$ and
$\beta_{rx}=\beta_{ry}=\beta_{r}$, the PDF and CDF of pointing error are
obtained as:

$\displaystyle
f_{h_{p}}(h_{p})=\frac{1}{G_{0}(\beta_{t}-\beta_{r})}\left[\left(\frac{h_{p}}{G_{0}}\right)^{1/\beta_{t}-1}-\left(\frac{h_{p}}{G_{0}}\right)^{1/\beta_{r}-1}\right],$
(15) $\displaystyle
F_{h_{p}}(h_{p})=\frac{1}{\beta_{t}-\beta_{r}}\left[\beta_{t}\left(\frac{h_{p}}{G_{0}}\right)^{1/\beta_{t}}-\beta_{r}\left(\frac{h_{p}}{G_{0}}\right)^{1/\beta_{r}}\right].$
(16)

###### Proof:

Please refer to Appendix B. ∎

The results of Theorem 2 are for a particular case where the severity of
vibrations along the Yaw and Pitch directions are approximately the same or
$\sigma_{\theta_{tx}}\simeq\sigma_{\theta_{ty}}=\sigma_{\theta}$ and
$\sigma_{\theta_{rx}}\simeq\sigma_{\theta_{ry}}=\sigma_{\theta}$. For example,
for communication between UAVs, when both UAVs are of the same model and at
the same altitude, they are expected to have similar fluctuation intensity. As
we observe from Theorem 2, the closed-form expressions provided in Eqs. (15)
and (16) are very tractable and simple. In addition, the simulation results of
Fig. 4 verify the accuracy of the analytical expressions under the different
values $\sigma_{\theta_{qw}}$s.

(a)

(b)

Figure 4: Comparison of the accuracy of the provided analytical results in
Theorem 2 with Monte-Carlo simulations for three different states of the
intensity of transceiver vibrations (a) for the PDF, and (b) for the CDF.

The results of Theorems 1 and 2 are for the planner array antenna. In some
situations, such as communications between Low-altitude UAVs, the use of a
linear array antenna in the direction perpendicular to the ground (as
illustrated in Fig. 5) can have several advantages compared to the planar
array antenna. The first important feature of using linear area antenna is
that the considered communication link will be resistant to vibrations in the
direction of the Yaw or x-z plane. This topology has an approximately flat
pattern in the x-z plane and therefore, the vibrations of the UAVs in this
direction do not have a significant effect on the system performance.
Moreover, this topology has a high-gain pattern in the direction of the Pitch
or y-z plane, and thus, interference caused by ground transmitters is reduced
significantly, and vice versa. In addition, from an aerodynamic point of view,
it is a better option for mounting on UAVs than the planar antennas,
especially for fixed-wing UAVs that are in motion.

(a)

(b)

Figure 5: Illustration of considered configuration based on ULA antenna in
such a way that the horizontal pattern, which is almost flat, is horizontal to
the ground, and the vertical pattern of the antenna is perpendicular to the
ground devices: (a) a better illustration of vertical and horizontal patterns,
and (b) a real 3D antenna pattern.

Theorem 3. When using ULA antennas, the PDF and CDF of pointing error are
obtained as:

$\displaystyle
f_{h_{p}}(h_{p})=\frac{1}{\sqrt{\beta_{ty}\beta_{ry}}}G_{0}^{\frac{2\beta_{ty}\beta_{ry}}{\beta_{ty}+\beta_{ry}}}h_{p}^{\frac{\beta_{ty}+\beta_{ry}}{2\beta_{ty}\beta_{ry}}-1}$
$\displaystyle~{}~{}~{}~{}\times
I_{0}\left(\frac{\beta_{ty}-\beta_{ry}}{2\beta_{ty}\beta_{ry}}\ln\left(\frac{h_{p}}{G_{0}}\right)\right)$
(17) $\displaystyle
F_{h_{p}}(h_{p})=1-Q\left(\Re_{1}\sqrt{\frac{-2\ln\left(\frac{h_{p}}{G_{0}}\right)}{\beta_{ty}+\beta_{ry}}},\Re_{2}\sqrt{\frac{-2\ln\left(\frac{h_{p}}{G_{0}}\right)}{\beta_{ty}+\beta_{ry}}}\right)$
$\displaystyle+Q\left(\Re_{2}\sqrt{\frac{-2\ln\left(\frac{h_{p}}{G_{0}}\right)}{\beta_{ty}+\beta_{ry}}},\Re_{1}\sqrt{\frac{-2\ln\left(\frac{h_{p}}{G_{0}}\right)}{\beta_{ty}+\beta_{ry}}}\right)$
(18)

where $Q(a,b)$ is the Marcum Q-function [20], and

$\displaystyle\left\\{\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{rl}&\Re_{1}=\frac{\sqrt{1-T_{q}^{4}}}{2T_{q}}\sqrt{\frac{1+T_{q}}{1-T_{q}}},\\\
&\Re_{2}=\Re_{1}\frac{1-T_{q}}{1+T_{q}},~{}~{}~{}T_{q}=\frac{\sigma_{\theta_{\text{min}}}}{\sigma_{\theta_{\text{max}}}},\\\
&\sigma_{\theta_{\text{max}}}=\text{max}\\{\sigma_{\theta_{q_{x}}},\sigma_{\theta_{q_{y}}}\\},~{}~{}~{}~{}\sigma_{\theta_{\text{min}}}=\text{min}\\{\sigma_{\theta_{q_{x}}},\sigma_{\theta_{q_{y}}}\\}\end{array}\right.$
(22)

###### Proof:

Please refer to Appendix C. ∎

As can be seen, the modified Bessel function of the first kind and the Gamma
function have been used in Eqs. (III) and (III). In order to reduce the
complexity of calculations, an approximate model is derived below, which uses
only a series of simple multiplication and addition operators to model the PDF
and CDF of pointing errors.

Lemma 1. For linear array antennas, the PDF and CDF of pointing error can be
well approximated as:

$\displaystyle
f_{h_{p}}(h_{p})\simeq\sum_{n=1}^{N}\frac{\Re_{4}\Re_{5}}{\Re_{3}(n)(\beta_{ty}+\beta_{ry})G_{0}}\left(\frac{h_{p}}{G_{0}}\right)^{\frac{\Re_{4}}{\beta_{ty}+\beta_{ry}}-1}$
(23) $\displaystyle
F_{h_{p}}(h_{p})\simeq\sum_{n=1}^{N}\frac{\Re_{5}}{\Re_{3}(n)}\left(\frac{h_{p}}{G_{0}}\right)^{\frac{\Re_{4}}{\beta_{ty}+\beta_{ry}}}$
(24)

where

$\displaystyle\left\\{\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{rl}&\Re_{3}(n)=1+\frac{1-T_{q}^{2}}{1+T_{q}^{2}}\cos\left(\pi\frac{2n-1}{N}\right),\\\
&\Re_{4}=\frac{(1+T_{q}^{2})^{2}}{2T_{q}^{2}}\Re_{3}(n),~{}~{}~{}\Re_{5}=\frac{2T_{q}}{N(1+T_{q}^{2})}.\end{array}\right.$
(27)

###### Proof:

Please refer to Appendix C. ∎

The closed-form expressions provided in Lemma 1 are very tractable and simple.
In addition, the simulation results of Fig. 6 verify the accuracy of the
analytical expressions of Theorem 3 and Lemma 1 under different values
$\sigma_{\theta_{qw}}$s.

(a)

(b)

Figure 6: Comparison of the accuracy of the provided analytical results in
Theorem 3 and Lemma 1 for ULA antenna with Monte-Carlo simulations for three
different states of the intensity of transceiver vibrations (a) the PDF, and
(b) the CDF.

In some scenarios, we may have a communication link between a stable node and
an unstable node. For example, for a ground-to-UAV or a UAV-to-ground station,
the intensity of vibrations of the ground node is negligible compared to the
UAV node. In this case, the pointing error model provided in Theorem 3 and
Lemma 1 can be used with a series of minor changes that are described in
Remark 1.

Remark 1. By following the method used in Appendices A and C, it can be easily
shown that the analytical expressions in Theorem 3 and Lemma 1 can be used for
a ground-to-UAV link by substituting $\beta_{rx}$ instead of $\beta_{ty}$.
Also, by substituting $\beta_{tx}$ instead of $\beta_{ry}$, we can use Theorem
3 and Lemma 1 for modeling the pointing error of a UAV-to-ground link.

##  IV End-to-End channel Modeling and Performance Analysis

The channel path loss, due to free-space propagation and the molecular
absorption loss, is a deterministic parameter and does not change the
distribution function of the considered channel. However, in addition to the
random parameter due to the pointing error, in some high-frequency typologies,
we may also have the effect of small-scale fading. Small-scale fading can be
caused by various factors such as atmospheric turbulence in rainy weather or
multi-path fading for ground links or close to the ground. In this work, we
use the $\alpha-\mu$ distribution to model the small-scale fading, which is
the most common model for modeling the small-scale fading of high frequency
channels. In the sequel, for the pointing error models provided in Theorems
1-3, we derive the end-to-end channel distribution function of the considered
communication link in Theorems 4-6 by considering $\alpha-\mu$ distribution
for small-scale fading.

(a)

(b)

(c)

Figure 7: Comparison of the accuracy of the provided analytical results with
Monte-Carlo simulations: (a) for the general scenario provided in Theorem 4,
(b) for a special case provided in Theorem 5 and Lemma 2, and (c) for ULA
antenna provided in Theorem 6.

Theorem 4. Under the general condition when the parameters
$\sigma_{\theta_{qw}}$s and $w_{Bq}$s have different values, the PDF of $h$ is
obtained as:

$\displaystyle
f_{h}(h)\simeq\frac{A_{1}C_{g}}{2G_{0}^{\alpha\mu}h_{L}}\left(\frac{h}{h_{L}}\right)^{\alpha\mu-1}\sum_{k=0}^{K}\frac{\Delta_{k}}{\Gamma(k+2)\beta_{1}^{k+2}}e^{-\frac{\mathbb{A}_{1}(h)}{\alpha^{2}}}$
$\displaystyle\times\frac{1}{\mathbb{A}_{1}(h)^{k/2+1}}\times\left[\Gamma\left(\frac{k+2}{2}\right)~\\!\\!{}_{1}F_{1}\left(\frac{k+2}{2},\frac{1}{2},\frac{\mathbb{A}_{2}(h)}{4}\right)+\right.$
$\displaystyle\left.\frac{\alpha
A_{3}-\mathbb{A}_{1}(h)}{\alpha\sqrt{\mathbb{A}_{1}(h)}}\Gamma\left(\frac{k+3}{2}\right)~\\!\\!{}_{1}F_{1}\left(\frac{k+3}{2};\frac{3}{2};\frac{\mathbb{A}_{2}(h)}{4}\right)\right]$
(28)

where

$\displaystyle\left\\{\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{rl}&\mathbb{A}_{0}(h)=\left(\frac{h}{G_{0}h_{L}}\right)^{\alpha},~{}~{}~{}\mathbb{A}_{1}(h)=\frac{A_{2}\alpha^{2}}{2}\mathbb{A}_{0}(h),\\\
&\mathbb{A}_{2}(h)=\frac{\left(A_{3}-\mathbb{A}_{1}(h)/\alpha\right)^{2}}{\mathbb{A}_{1}(h)},\\\
&A_{1}=\frac{\alpha\mu^{\mu}}{\hat{h_{a}}^{\alpha\mu}\Gamma(\mu)},~{}~{}~{}A_{2}=\frac{\mu}{\hat{h_{a}}^{\alpha}},~{}~{}~{}A_{3}=\alpha\mu-1/\beta_{1}\end{array}\right.$
(32)

Also, the CDF of $h$ is derived as

$\displaystyle
F_{h}(h)=1-C_{g}\sum_{m=0}^{\mu-1}\sum_{k=0}^{K}\frac{A_{2}^{m}}{\Gamma(m+1)}\frac{\Delta_{k}}{\Gamma(k+2)\beta_{1}^{k+2}}\mathbb{A}_{0}(h)^{m}$
$\displaystyle\times
e^{-A_{2}\mathbb{A}_{0}(h)}\sum_{j=0}^{k+1}2^{j-k-2}\binom{k+1}{j}\left(m\alpha-\frac{1}{\beta_{1}}-A_{2}\mathbb{A}_{0}(h)\alpha\right)^{k-j+1}$
$\displaystyle\times\left(\frac{\alpha^{2}A_{2}\mathbb{A}_{0}(h)}{2}\right)^{\frac{j-3}{2}-k}\exp\left(\frac{\left(m\alpha-\frac{1}{\beta_{1}}-A_{2}\mathbb{A}_{0}(h)\alpha\right)^{2}}{2\left(\alpha^{2}A_{2}\mathbb{A}_{0}(h)\right)}\right)$
$\displaystyle\times\Gamma\left(\frac{j+1}{2},\frac{\left(m\alpha-\frac{1}{\beta_{1}}-A_{2}\mathbb{A}_{0}(h)\alpha\right)^{2}}{2\left(\alpha^{2}A_{2}\mathbb{A}_{0}(h)\right)}\right)$
(33)

###### Proof:

Please refer to Appendix D. ∎

Theorem 5. When $\beta_{tx}=\beta_{ty}=\beta_{t}$ and
$\beta_{rx}=\beta_{ry}=\beta_{r}$, the PDF and CDF of $h$ are obtained as:

$\displaystyle
f_{h}(h)=\frac{1}{\alpha(\beta_{t}-\beta_{r})}\frac{1}{A_{2}h_{L}G_{0}}A_{1}\left(\mathbb{A}_{0}(h)\right)^{\mu-\frac{1}{\alpha}-1}$
$\displaystyle\times\bigg{[}(A_{2}\mathbb{A}_{0}(h))^{\frac{V_{1}}{2}}e^{-\frac{A_{2}\mathbb{A}_{0}(h)}{2}}\mathbb{W}_{\frac{-V_{1}}{2},\frac{1-V_{1}}{2}}(A_{2}\mathbb{A}_{0}(h))$
$\displaystyle-(A_{2}\mathbb{A}_{0}(h))^{\frac{V_{2}}{2}}e^{-\frac{A_{2}\mathbb{A}_{0}(h)}{2}}\mathbb{W}_{\frac{-V_{2}}{2},\frac{1-V_{2}}{2}}(A_{2}\mathbb{A}_{0}(h))\bigg{]}$
(34)

and

$\displaystyle F_{h}(h)$
$\displaystyle=1-\frac{1}{\alpha(\beta_{t}-\beta_{r})}\sum_{k=0}^{\mu-1}\frac{\left(A_{2}\mathbb{A}_{0}(h)\right)^{k-1}}{\Gamma(k+1)}e^{-\frac{A_{2}\mathbb{A}_{0}(h)}{2}}$
(35)
$\displaystyle~{}~{}~{}\times\bigg{[}\left(A_{2}\mathbb{A}_{0}(h)\right)^{\frac{V_{3}}{2}}\mathbb{W}_{\frac{-V_{3}}{2},\frac{1-V_{3}}{2}}\left(A_{2}\mathbb{A}_{0}(h)\right)$
$\displaystyle~{}~{}~{}-\left(A_{2}\mathbb{A}_{0}(h)\right)^{\frac{V_{4}}{2}}\mathbb{W}_{\frac{-V_{4}}{2},\frac{1-V_{4}}{2}}\left(A_{2}\mathbb{A}_{0}(h)\right)\bigg{]}$
(36)

where $\mathbb{W}_{-\frac{\nu}{2},\frac{1-\nu}{2}}(u)$ is the Whittaker
function [21], $V_{1}=(1/\alpha\beta_{t}-\mu+1)$,
$V_{2}=(1/\alpha\beta_{r}-\mu+1)$, $V_{3}=(1/\alpha\beta_{t}-k+1)$, and
$V_{4}=(1/\alpha\beta_{r}-k+1)$.

###### Proof:

Please refer to Appendix E. ∎

Lemma 2. When $\beta_{tx}=\beta_{ty}=\beta_{t}$, and
$\beta_{rx}=\beta_{ry}=\beta_{r}$, and under the condition that
$\big{(}\alpha\mu\min\\{\beta_{t},\beta_{r}\\}\big{)}>1$, the PDF of $h$ can
be approximated as:

$\displaystyle f_{h}(h)$
$\displaystyle=\frac{1}{\alpha(\beta_{t}-\beta_{r})}\frac{A_{1}\left(\mathbb{A}_{0}(h)\right)^{\mu-\frac{1}{\alpha}-1}}{A_{2}h_{L}G_{0}}$
$\displaystyle~{}~{}~{}\times\Big{[}(A_{2}\mathbb{A}_{0}(h))^{1-s_{1}}\Gamma\left(s_{1},A_{2}\mathbb{A}_{0}(h)\right)$
$\displaystyle~{}~{}~{}-(A_{2}\mathbb{A}_{0}(h))^{1-s_{2}}\Gamma\left(s_{2},A_{2}\mathbb{A}_{0}(h)\right)\Big{]}$
(37)

where $\Gamma(s,x)$, is the upper incomplete gamma function [22],
$s_{1}=(-1/\alpha\beta_{t}+\mu)$, and $s_{2}=(-1/\alpha\beta_{r}+\mu)$.

###### Proof:

Please refer to Appendix E. ∎

The results of Eq. (IV) are based on the upper incomplete gamma function
$\Gamma(s,x)$, which has faster computation than the Whittaker function.
However, it should be noted that for $\Gamma(s,x)$, the parameter $s$ must be
positive. In Eq. (IV), both $s_{1}$ and $s_{2}$ are function of the parameters
$\alpha$ and $\mu$ related to the small-scale distribution, as well as are
function of the parameters $\beta_{t}$ and $\beta_{r}$. The parameters
$\beta_{t}$ and $\beta_{r}$ characterize the strength of pointing error. For
smaller values of $N_{q}$ and/or smaller values of $\sigma_{\theta_{qw}}$, the
parameters $\beta_{t}$ and $\beta_{r}$ decrease, and thus, the values of
parameters $s_{1}$ and $s_{2}$ becomes negative. As a results, smaller values
of $N_{q}$ and $\sigma_{\theta_{qw}}$, the results of Lemma 2 are not valid.
However, in a wide range of channel parameter values, the results of Lemma 2
are valid.

Theorem 6. For linear array antenna, the PDF of $h$ is obtained as:

$\displaystyle
f_{h}(h)=\sum_{n=1}^{N}\frac{A_{1}\Re_{4}\Re_{5}}{\alpha\Re_{3}(n)(\beta_{ty}+\beta_{ry})}\frac{1}{A_{2}h_{L}G_{0}}(\mathbb{A}_{0}(h))^{\mu-1/\alpha-1}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}\times(A_{2}\mathbb{A}_{0}(h))^{\frac{V_{5}}{2}}e^{-\frac{A_{2}\mathbb{A}_{0}(h)}{2}}\mathbb{W}_{\frac{-V_{5}}{2},\frac{1-V_{5}}{2}}(A_{2}\mathbb{A}_{0}(h))$
(38)

and the CDF is calculated as

$\displaystyle
F_{h}(h)=1-\sum_{k=0}^{\mu-1}\sum_{n=1}^{N}\frac{1}{\alpha\Gamma(k+1)}\frac{\Re_{4}\Re_{5}}{\Re_{3}(n)(\beta_{ty}+\beta_{ry})}$
$\displaystyle\times\left(A_{2}\mathbb{A}_{0}(h)\right)^{k+\frac{V_{6}}{2}-1}e^{-\frac{A_{2}\mathbb{A}_{0}(h)}{2}}\mathbb{W}_{-\frac{V_{6}}{2},\frac{1-V_{6}}{2}}\left(A_{2}\mathbb{A}_{0}(h)\right)$
(39)

where $V_{5}=\frac{\Re_{4}}{\alpha(\beta_{ty}+\beta_{ry})}-\mu+1)$, and
$V_{6}=\frac{\Re_{4}}{\alpha(\beta_{ty}+\beta_{ry})}-k+1)$.

###### Proof:

Please refer to Appendix F. ∎

In Fig. 7, the accuracy of the analytical expressions provided in Theorems 4-6
is examined using Monte-Carlo simulations. Here, we set $\alpha=2$, $\mu=4$,
$f_{c}=280$GHz, $\mathcal{K}(f)=2$, and other parameters related to the
antenna pattern and pointing errors, as specified in Fig. 7. According to the
result obtained from Fig. 7a, although the analytical results obtained from
Theorem 4 are close to the simulation results with an acceptable accuracy,
however, a small gap is observed between the simulation results and the
results obtained from the approximated analytical expressions in Theorem 2.
The results of Figs. 7b and 7c also show that the analytical results obtained
from Theorems 5 and 6 completely match the simulation results.

Finally, in Figs. 8 and 9, we evaluate the performance of the considered
communication link in term of the outage probability. In the simulations, we
set transmit power $P_{t}=10$mW, and SNR threshold $\gamma_{th}=5$dB. Also,
for computing thermal noise power, we consider bandwidth $\Delta_{f}=100$MHz,
and receiver’s temperature $T=300$K. The results of Fig. 8 are plotted for
different values of molecular absorption loss characterized by
$\mathcal{K}(f)=0.5$, 1, 2, and 4. For example, if we consider that for a
quality of service (QoS) with target outage probability lower than $10^{-2}$,
it is observed that by increasing $\mathcal{K}(f)$ from 0.5 to 4, the maximum
achievable link length to guarantee the requested QoS decreases from 3.1km to
lower than 1 km. There are several ways to increase the maximum achievable
link-length of the considered system under the presence of molecular
absorption loss. However, due to the limited transmission power at higher
frequencies, one of the best methods is to increase the gain of the antenna by
increasing parameter $N$. For this aim, the results of Fig. 9 are plotted for
different values of $N$. The results of Fig. 9 show an important point that
system performance does not necessarily improve with increasing $N$. The
reason is that, with increasing antenna gain, the beamwidth of antenna pattern
decreases and the system becomes more sensitive to the vibrations of Tx/Rx. To
better illustrate this point, the results in Fig. 9 are plotted for two
different values for $\sigma_{\theta_{qw}}$s. The first category is for
$\sigma_{\theta_{t}}=1^{o}$ and $\sigma_{\theta_{r}}=0.8^{o}$. For this state
and for target outage probability equal to $10^{-3}$, the optimal value is
$N=30$. The second category is for $\sigma_{\theta_{t}}=0.4^{o}$ and
$\sigma_{\theta_{r}}=0.3^{o}$ that we observe the maximum link length is
obtained for $N=60$ in the target outage probability $10^{-3}$.

##  V Conclusions

A general pointing error model was investigated for directional mmWave/THz
links by taking into account actual antenna patterns. Our derivation is based
on the practical assumption that the intensity of instability in Tx and Rx and
also the intensity of instability in the Yaw and Pitch directions are
different. For lower frequencies as well as higher heights, the small-scale
fading effect can be ignored with good accuracy, and as a result, the
presented models for pointing errors can be used for end-to-end system
analysis. For higher frequencies and in some scenarios, the effect of small-
scale fading should be considered. To this end, using $\alpha-\mu$
distribution, which is a common model for small-scale fading of THz links, the
end-to-end PDF and CDF of the considered channel are derived for both planer
and linear array antennas. The results of this work can be also used to
calculate the end-to-end channel distribution in the presence of newer small-
scale fading models. Finally, by employing Monte-Carlo simulations, the
accuracy of the analytical expressions was verified and the performance of the
system was studied.

Figure 8: Outage probability of the considered communication link versus link
length for different values of molecular absorption loss characterized by
$\mathcal{K}(f)$. Figure 9: Outage probability of the considered communication
link versus link length for two different levels of system instability and
different values for $N$.

## Appendix A

The main lobe of the antenna pattern can be well approximated by the Gaussian
distribution function as [13]

$\displaystyle G_{q}(N_{q},\theta_{qx},\theta_{qy})=G_{0}(N_{q})$
$\displaystyle\times\exp\left(-\frac{\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{qx})+\tan^{2}(\theta_{qy})}\right)\right)^{2}}{w^{2}_{B}(N_{q})}\right),$
(40)

where $G_{0}(N_{q})=\pi N_{q}^{2}$, and $w_{B}(N_{q})=\frac{1.061}{N_{q}}$ is
the angular beamwidth (called the beam divergence). For lower values of
$\theta_{qx}$ and $\theta_{qy}$, we can approximate RV $\theta_{q}$ as

$\displaystyle\theta_{q}\simeq\sqrt{\theta_{qx}^{2}+\theta_{qy}^{2}}.$ (41)

Substituting Eqs. (41) and (A) in Eq. (8), we have

$\displaystyle h_{p}\simeq
G_{0}(N)\exp\left(-\frac{\theta_{t_{x}}^{2}+\theta_{t_{y}}^{2}}{2w^{2}_{Bt}}-\frac{\theta_{r_{x}}^{2}+\theta_{r_{y}}^{2}}{2w^{2}_{Br}}\right)$
(42)

Let us define parameter $\Theta$ as

$\displaystyle\Theta=\Theta_{t_{x}}+\Theta_{t_{y}}+\Theta_{r_{x}}+\Theta_{r_{y}}$
(43)

where $\Theta_{q_{w}}=\frac{\theta_{qw}^{2}}{2w^{2}_{Bq}}$. Since
$\theta_{qw}\sim\mathcal{N}(0,\sigma_{\theta_{qw}})$, then $\Theta_{q_{w}}$
follows a gamma distribution as [23]

$\displaystyle
f_{\Theta_{q_{w}}}(\Theta_{q_{w}})=\frac{1}{\sqrt{\beta_{qw}\Theta_{q_{w}}}\Gamma(1/2)}\exp\left(-\frac{\Theta_{q_{w}}}{\beta_{qw}}\right)$
(44)

where $\beta_{qw}=\frac{\sigma_{\theta_{q_{w}}^{2}}}{w^{2}_{Bq}}$. Let’s
arrange the parameters $\beta_{qw}$s in ascending order and then define the
vector $\bar{\beta}=[\beta_{1},...,\beta_{4}]$, where $\beta_{1}$ is the least
value of $\beta_{q_{w}}$s and $\beta_{4}$ is the maximum value of
$\beta_{q_{w}}$s. Based on the results of [24, Theorem 1], the distribution of
RV $\Theta$ is obtained as:

$\displaystyle
f_{\Theta}(\Theta)=C_{g}\sum_{k=0}^{K}\frac{\Delta_{k}\Theta^{k+1}\exp\left(-\frac{\Theta}{\beta_{1}}\right)}{\Gamma(k+2)\beta_{1}^{k+2}}$
(45)

where

$\displaystyle\left\\{\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{rl}&C_{g}=\prod_{i=1}^{4}\sqrt{\beta_{1}/\beta_{i}},\\\
&\Delta_{k}=\frac{1}{k}\sum_{i=1}^{k}i\gamma_{i}\Delta_{k+1-i}~{}~{}\text{for}~{}~{}k=1,...,K,\\\
&\Delta_{0}=1,~{}~{}~{}\gamma_{k}=\sum_{i=1}^{4}\frac{\left(1-{\beta_{1}}/{\beta_{i}}\right)^{k}}{2k}.\end{array}\right.$
(49)

From Eq. (42), the distribution of $h_{p}$ can be obtained as [23]

$\displaystyle
f_{h_{p}}(h_{p})=f_{\Theta}\left(-\ln\left(\frac{h_{p}}{G_{0}}\right)\right)\frac{\text{d}\Theta}{\text{d}h_{p}}=\frac{1}{h_{p}}f_{\Theta}\left(-\ln\left(\frac{h_{p}}{G_{0}}\right)\right)$
(50)

Using Eqs. (45) and (50), the closed-form expression for PDF of $h_{p}$ is
derived in Eq. (9).

From Eq. (45), the CDF of RV $\Theta$ is obtained as

$\displaystyle
F_{\Theta}(\Theta)=C_{g}\sum_{k=0}^{K}\frac{\Delta_{k}}{\Gamma(k+2)\beta_{1}^{k+2}}\int_{0}^{\Theta}x^{k+1}\exp\left(-\frac{x}{\beta_{1}}\right)$
(51)

In the following derivation, we use an integral identity [21, eq. (2.321.2)]:

$\displaystyle\int
x^{n}e^{ax}\text{d}x=e^{ax}\left(\sum_{j=0}^{n}\frac{(-1)^{j}j!\binom{n}{j}}{a^{j+1}}x^{n-j}\right)$
(52)

Using Eq. (52), the CDF of RV $\Theta$ is derived as

$\displaystyle
F_{\Theta}(\Theta)=\Theta_{0}-C_{g}e^{-\frac{\Theta}{\beta_{1}}}\sum_{k=0}^{K}\sum_{j=0}^{k+1}\frac{\Delta_{k}\Theta^{k+1-j}}{\Gamma(k+2-j)\beta_{1}^{k+1-j}}$
(53)

where $\Theta_{0}=C_{g}\sum_{k=0}^{K}\Delta_{k}$. Based on Eq. (42), the CDF
of $h_{p}$ is obtained as

$\displaystyle
F_{h_{p}}(h_{p})=\\!\text{Prob}\left\\{\Theta>-\\!\ln\left(\frac{h_{p}}{G_{0}}\right)\right\\}=1\\!-\\!F_{\Theta}\left(\\!-\\!\ln\left(\frac{h_{p}}{G_{0}}\right)\right)$
(54)

Finally, using Eqs. (53) and (54), the closed-form expression for CDF of
$h_{p}$ is derived in Eq. (III).

## Appendix B

In this appendix, we obtain the PDF and CDF of $h_{p}$ for the scenario where
the intensity of the vibrations along the Yaw and Pitch axes is the same,
approximately. For this case, using Eqs. (8), (A), and (41), the RV $h_{p}$
can be formulated as

$\displaystyle h_{p}\simeq
G_{0}(N)\exp\left(-\frac{\theta_{t}^{2}}{2w^{2}_{Bt}}-\frac{\theta_{r}^{2}}{2w^{2}_{Br}}\right).$
(55)

where $\theta_{q}=\sqrt{\theta_{qx}^{2}+\theta_{qy}^{2}}$ and since
$\theta_{qw}\sim\mathcal{N}(0,\sigma_{\theta_{qw}})$, the RV $\theta_{qw}$
follows Rayleigh distribution as [23]

$\displaystyle
f_{\theta_{q}}(\theta_{q})=\frac{\theta_{q}}{\sigma_{\theta_{q}}^{2}}\exp\left(-\frac{\theta_{q}^{2}}{2\sigma_{\theta_{q}}^{2}}\right)$
(56)

Let us define
$\theta_{tr}=\frac{\theta_{t}^{2}}{2w^{2}_{Bt}}+\frac{\theta_{r}^{2}}{2w^{2}_{Br}}$.
The CDF of RV $\theta_{tr}$ conditioned on $\theta_{r}$ is obtained as:

$\displaystyle F_{\theta_{tr}|\theta_{r}}(\theta_{tr})$
$\displaystyle=\text{Prob}\left\\{\frac{\theta_{t}^{2}}{2w^{2}_{Bt}}+\frac{\theta_{r}^{2}}{2w^{2}_{Br}}<\theta_{tr}|\theta_{r}\right\\}$
$\displaystyle=F_{t}\left(\sqrt{2w^{2}_{Bt}\theta_{tr}-\frac{w^{2}_{Bt}}{w^{2}_{Br}}\theta_{r}^{2}}\right)$
(57)

Taking the derivative of Eq. (B) with respect to $\theta_{tr}$, the PDF of RV
$\theta_{tr}$ conditioned on RV $\theta_{r}$ is obtained as

$\displaystyle f_{\theta_{tr}|\theta_{r}}(\theta_{tr})$
$\displaystyle=\frac{w^{2}_{Bt}}{\sqrt{2w^{2}_{Bt}\theta_{tr}-\frac{w^{2}_{Bt}}{w^{2}_{Br}}\theta_{r}^{2}}}$
$\displaystyle~{}~{}~{}\times
f_{\theta_{t}}\left(\sqrt{2w^{2}_{Bt}\theta_{tr}-\frac{w^{2}_{Bt}}{w^{2}_{Br}}\theta_{r}^{2}}\right)$
(58)

Using Eqs. (56) and (B), and after some manipulations, the PDF of
$\theta_{tr}$ is derived as

$\displaystyle f_{\theta_{tr}}(\theta_{tr})$
$\displaystyle=\frac{w^{2}_{Bt}w^{2}_{Br}}{w^{2}_{Br}\sigma_{\theta_{t}}^{2}-w^{2}_{Bt}\sigma_{\theta_{r}}^{2}}$
$\displaystyle\times\left[\exp\left(-\frac{w^{2}_{Bt}}{\sigma_{\theta_{t}}^{2}}\theta_{tr}\right)-\exp\left(-\frac{w^{2}_{Br}}{\sigma_{\theta_{r}}^{2}}\theta_{tr}\right)\right]$
(59)

Now, using Eqs. (50) and (B), and after some derivations, the PDF of pointing
error is derived in Eq. (15). By integrating Eq. (B), the CDF of $\theta_{tr}$
is obtained as

$\displaystyle F_{\theta_{tr}}(\theta_{tr})$
$\displaystyle=\frac{w^{2}_{Bt}w^{2}_{Br}}{w^{2}_{Br}\sigma_{\theta_{t}}^{2}-w^{2}_{Bt}\sigma_{\theta_{r}}^{2}}\left[\frac{\sigma_{\theta_{t}}^{2}}{w^{2}_{Bt}}\left(1-\exp\left(-\frac{w^{2}_{Bt}}{\sigma_{\theta_{t}}^{2}}\theta_{tr}\right)\right)\right.$
$\displaystyle~{}~{}~{}\left.-\frac{\sigma_{\theta_{r}}^{2}}{w^{2}_{Br}}\left(1-\exp\left(-\frac{w^{2}_{Br}}{\sigma_{\theta_{r}}^{2}}\theta_{tr}\right)\right)\right].$
(60)

Finally, using Eqs. (54) and (B), and after some derivations, the CDF of
$h_{p}$ is derived in Eq. (16).

## Appendix C

In this section, we derive the PDF and CDF of $h_{p}$ for the scenario in
which the Tx and Rx are equipped with ULA antennas. For ULA antennas, the
antenna pattern model provided in Eq. (II) is simplified as [17]

$\displaystyle G^{\prime}_{q}(N_{q},\theta_{q},\phi_{q})$
$\displaystyle=\left(\frac{\sin\left(\frac{N_{q}(kd_{y}\sin(\theta_{q})\sin(\phi_{q}))}{2}\right)}{N_{q}\sin\left(\frac{kd_{y}\sin(\theta_{q})\sin(\phi_{q})}{2}\right)}\right)^{2}.$
(61)

For this scenario, using Eqs. (8), (A), and (61), the RV $h_{p}$ is formulated
as

$\displaystyle h_{p}\simeq
G_{0}\exp\left(-\frac{\theta_{q_{x}}^{2}+\theta_{q_{y}}^{2}}{2w^{2}_{Bq}}\right)$
(62)

Let us define
$\Theta=\frac{\theta_{q_{x}}^{2}}{2w^{2}_{Bq}}+\frac{\theta_{q_{y}}^{2}}{2w^{2}_{Bq}}$.
Using [25, eqs. (3) and (5)], the distribution of $\Theta$ is obtained as

$\displaystyle
f_{\Theta}(\Theta)=\frac{1}{\sqrt{\beta_{qx}\beta_{qy}}}\exp\left(-\frac{\beta_{qx}+\beta_{qy}}{2\beta_{qx}\beta_{qy}}\Theta\right)I_{0}\left(\frac{\beta_{qx}-\beta_{qy}}{2\beta_{qx}\beta_{qy}}\Theta\right)$
(63)

where $I_{0}(\cdot)$ is the modified Bessel function of the first kind with
order zero. Using Eqs. (50) and (63), the PDF of $h_{p}$ is derived as

$\displaystyle f_{h_{p}}(h_{p})$
$\displaystyle=\frac{1}{\sqrt{\beta_{qx}\beta_{qy}}}G_{0}^{\frac{2\beta_{qx}\beta_{qy}}{\beta_{qx}+\beta_{qy}}}h_{p}^{\frac{\beta_{qx}+\beta_{qy}}{2\beta_{qx}\beta_{qy}}-1}$
$\displaystyle~{}~{}~{}\times
I_{0}\left(\frac{\beta_{qy}-\beta_{qx}}{2\beta_{qx}\beta_{qy}}\ln\left(\frac{h_{p}}{G_{0}}\right)\right)$
(64)

Based on Eq. (C) and using [26], the CDF of $\Theta$ is derived as

$\displaystyle F_{\Theta}(\Theta)$
$\displaystyle=Q\left(\Re_{1}\sqrt{\frac{2\Theta}{\beta_{qx}+\beta_{qy}}},\Re_{2}\sqrt{\frac{2\Theta}{\beta_{qx}+\beta_{qy}}}\right)$
$\displaystyle~{}~{}~{}-Q\left(\Re_{2}\sqrt{\frac{2\Theta}{\beta_{qx}+\beta_{qy}}},\Re_{1}\sqrt{\frac{2\Theta}{\beta_{qx}+\beta_{qy}}}\right)$
(65)

where where $Q(a,b)$ is the Marcum Q-function, and

$\displaystyle\left\\{\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{rl}&\Re_{1}=\frac{\sqrt{1-T_{q}^{4}}}{2T_{q}}\sqrt{\frac{1+T_{q}}{1-T_{q}}},\\\
&\Re_{2}=\Re_{1}\frac{1-T_{q}}{1+T_{q}},~{}~{}~{}T_{q}=\frac{\sigma_{\theta_{\text{min}}}}{\sigma_{\theta_{\text{max}}}},\\\
&\sigma_{\theta_{\text{max}}}=\text{max}\\{\sigma_{\theta_{q_{x}}},\sigma_{\theta_{q_{y}}}\\},~{}~{}~{}~{}\sigma_{\theta_{\text{min}}}=\text{min}\\{\sigma_{\theta_{q_{x}}},\sigma_{\theta_{q_{y}}}\\}\end{array}\right.$
(69)

Substituting Eq. (C) in Eq. (54), the CDF of $h_{p}$ is obtained as

$\displaystyle F_{h_{p}}(h_{p})$
$\displaystyle=1-Q\left(\Re_{1}\sqrt{\frac{-2\ln\left(\frac{h_{p}}{G_{0}}\right)}{\beta_{qx}+\beta_{qy}}},\Re_{2}\sqrt{\frac{-2\ln\left(\frac{h_{p}}{G_{0}}\right)}{\beta_{qx}+\beta_{qy}}}\right)$
$\displaystyle+Q\left(\Re_{2}\sqrt{\frac{-2\ln\left(\frac{h_{p}}{G_{0}}\right)}{\beta_{qx}+\beta_{qy}}},\Re_{1}\sqrt{\frac{-2\ln\left(\frac{h_{p}}{G_{0}}\right)}{\beta_{qx}+\beta_{qy}}}\right)$
(70)

Now in the continuation of this section, we provide simpler analytical
expressions for the PDF and CDF of $h_{p}$. Using [27, eq. (5)], the CDF of RV
$\Theta$ can be well approximated as

$\displaystyle F_{\Theta}(\Theta)\simeq
1-\sum_{n=1}^{N}\frac{\Re_{5}}{\Re_{3}(n)}\exp\left(-\frac{\Re_{4}}{\beta_{qx}+\beta_{qy}}\Theta\right)$
(71)

where

$\displaystyle\left\\{\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{rl}&\Re_{3}(n)=1+\frac{1-T_{q}^{2}}{1+T_{q}^{2}}\cos\left(\pi\frac{2n-1}{N}\right),\\\
&\Re_{4}=\frac{(1+T_{q}^{2})^{2}}{2T_{q}^{2}}\Re_{3}(n),~{}~{}~{}\Re_{5}=\frac{2T_{q}}{N(1+T_{q}^{2})}.\end{array}\right.$
(74)

Using Eqs. (54) and (C), the CDF of $h_{p}$ is derived in Eq. (24). Finally,
by taking the derivative of Eq. (24), the PDF of $h_{p}$ is derived in Eq.
(23).

## Appendix D

In this appendix, we obtain the PDF and CDF of $h$ for the considered general
scenario. From Eq. (7) and [23], we have

$\displaystyle
f_{h}(h)=\int_{0}^{G_{0}}\frac{1}{h_{L}h_{p}}f_{h_{a}}\left(\frac{h}{h_{L}h_{p}}\right)f_{h_{p}}(h_{p})\text{d}h_{p}$
(75)

Substituting Eqs. (3) and (9) in Eq. (75), we obtain:

$\displaystyle
f_{h}(h)=\frac{A_{1}C_{g}}{G_{0}}\frac{1}{h_{L}}\left(\frac{h}{h_{L}}\right)^{\alpha\mu-1}\sum_{k=0}^{K}\frac{\Delta_{k}}{\Gamma(k+2)\beta_{1}^{k+2}}\times$
$\displaystyle\int_{0}^{G_{0}}\\!\\!\\!\exp\left(\\!\\!-A_{2}\left(\\!\frac{h}{h_{L}h_{p}}\\!\right)^{\alpha}\\!\right)\frac{h_{p}^{\frac{1}{\beta_{1}}-\alpha\mu-1}}{G_{0}^{1/\beta_{1}-1}}\left(\\!-\\!\ln\left(\frac{h_{p}}{G_{0}}\right)\\!\right)^{k+1}\\!\\!\text{d}h_{p}$
(76)

where $A_{1}=\frac{\alpha\mu^{\mu}}{\hat{h_{a}}^{\alpha\mu}\Gamma(\mu)}$, and
$A_{2}=\frac{\mu}{\hat{h_{a}}^{\alpha}}$. Applying a change of variable rule
$x=-\ln\left(\frac{h_{p}}{G_{0}}\right)$, we rewrite Eq. (D) as

$\displaystyle
f_{h}(h)=\frac{A_{1}C_{g}}{G_{0}^{\alpha\mu}}\frac{1}{h_{L}}\left(\frac{h}{h_{L}}\right)^{\alpha\mu-1}\sum_{k=0}^{K}\frac{\Delta_{k}}{\Gamma(k+2)\beta_{1}^{k+2}}$
$\displaystyle\times\int_{0}^{\infty}\exp\left(-A_{2}\left(\frac{h}{G_{0}h_{L}}\right)^{\alpha}e^{\alpha
x}\right)e^{(\alpha\mu-1/\beta_{1})x}x^{k+1}\text{d}x$ (77)

It can be easily shown that for $x>\alpha$, the above integral expression
tends to zero. Therefore, we can approximate $e^{\alpha x}=1+\alpha
x+\frac{\alpha^{2}x^{2}}{2}$. Based on this, we approximate Eq. (D) as

$\displaystyle
f_{h}(h)=\frac{A_{1}C_{g}}{G_{0}^{\alpha\mu}}\frac{1}{h_{L}}\left(\frac{h}{h_{L}}\right)^{\alpha\mu-1}\sum_{k=0}^{K}\frac{\Delta_{k}}{\Gamma(k+2)\beta_{1}^{k+2}}$
$\displaystyle\times\exp\left(-A_{2}\left(\frac{h}{G_{0}h_{L}}\right)^{\alpha}\right)\int_{0}^{\infty}e^{-A_{3}x^{2}+A_{4}x}x^{k+1}\text{d}x$
(78)

where
$A_{3}=\left(\frac{A_{2}}{2}\left(\frac{h}{G_{0}h_{L}}\right)^{\alpha}\alpha^{2}\right)$,
and
$A_{4}=\left(\alpha\mu-1/\beta_{1}-A_{2}\left(\frac{h}{G_{0}h_{L}}\right)^{\alpha}\alpha\right)$.
In the following derivation, we use an integral identity [21]

$\displaystyle\int_{0}^{\infty}x^{n}e^{-x^{2}+bx}\text{d}x=\frac{1}{2}\left[b\Gamma\left(\frac{n}{2}+1\right)~\\!\\!{}_{1}F_{1}\left(\frac{n}{2}+1;\frac{3}{2};\frac{b^{2}}{4}\right)\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.+\Gamma\left(\frac{n+1}{2}\right)~\\!\\!{}_{1}F_{1}\left(\frac{n+1}{2},\frac{1}{2},\frac{b^{2}}{4}\right)\right]$
(79)

Using Eqs. (D) and (D), and after some manipulations, the closed form
expressions for the PDF of $h$ is derived in Eq. (IV).

From Eq. (7) and [23], we have

$\displaystyle
F_{h}(h)=\int_{0}^{G_{0}}F_{h_{a}}\left(\frac{h}{h_{L}h_{p}}\right)f_{h_{p}}(h_{p})\text{d}h_{p}$
(80)

The CDF of $\alpha-\mu$ distribution can be expressed as [8, eqs. (66) and
(67)]

$\displaystyle
F_{h_{a}}(h_{a})=1-\sum_{k=0}^{\mu-1}\frac{A_{2}^{k}}{\Gamma(k+1)}h_{a}^{k\alpha}e^{-A_{2}h_{a}^{\alpha}}.$
(81)

Substituting Eqs. (9) and (81) in Eq. (80), we obtain:

$\displaystyle
F_{h}(h)=1-\frac{C_{g}}{G_{0}^{1/\beta_{1}}}\sum_{m=0}^{\mu-1}\sum_{k=0}^{K}\frac{A_{2}^{m}}{\Gamma(m+1)}\frac{\Delta_{k}}{\Gamma(k+2)\beta_{1}^{k+2}}\times$
$\displaystyle\int_{0}^{G_{0}}\left(\frac{h}{h_{L}h_{p}}\right)^{m\alpha}e^{-A_{2}\left(\frac{h}{h_{L}h_{p}}\right)^{\alpha}}h_{p}^{\frac{1}{\beta_{1}}-1}\left(\\!-\ln\left(\frac{h_{p}}{G_{0}}\right)\right)^{k+1}\\!\\!\text{d}h_{p}$
(82)

Applying a change of variable rule $x=-\ln\left(\frac{h_{p}}{G_{0}}\right)$,
Eq. (D) is expressed as

$\displaystyle F_{h}(h)$
$\displaystyle=1-\frac{C_{g}}{1}\sum_{m=0}^{\mu-1}\sum_{k=0}^{K}\frac{A_{2}^{m}}{\Gamma(m+1)}\frac{\Delta_{k}}{\Gamma(k+2)\beta_{1}^{k+2}}$
$\displaystyle~{}~{}~{}\times\int_{0}^{\infty}\left(\frac{h}{h_{L}G_{0}}\right)^{m\alpha}\exp\left(-A_{2}\left(\frac{h}{h_{L}G_{0}}\right)^{\alpha}e^{\alpha
x}\right)$ $\displaystyle~{}~{}~{}\times
e^{(m\alpha-\frac{1}{\beta_{1}})x}x^{k+1}\text{d}x$ (83)

In the following derivation, we use an integral identity [21]

$\displaystyle\int_{0}^{\infty}x^{n}e^{ax^{2}+bx}\text{d}x$
$\displaystyle=\sum_{j=0}^{n}2^{j-n-1}\binom{n}{j}b^{n-j}(-a)^{\frac{j-1}{2}-n}$
$\displaystyle~{}~{}~{}~{}\times\exp\left(-\frac{b^{2}}{4a}\right)\Gamma\left(\frac{j+1}{2},-\frac{b^{2}}{4a}\right)$
(84)

Finally, using Eqs. (D) and (D), and after some manipulations, the closed-form
expression for the CDF of $h$ is derived in Eq. (IV).

## Appendix E

In this appendix, we obtain the PDF and CDF of $h$ for the case where the
intensity of the vibrations along the Yaw and Pitch axes is the same,
approximately. Substituting Eqs. (3) and (15) in Eq. (75), we obtain:

$\displaystyle f_{h}(h)$
$\displaystyle=\frac{A_{1}}{G_{0}(\beta_{t}-\beta_{r})}\int_{0}^{G_{0}}\frac{1}{h_{L}h_{p}}$
$\displaystyle~{}~{}~{}\times\left(\frac{h}{h_{L}h_{p}}\right)^{\alpha\mu-1}\exp\left(-A_{2}\left(\frac{h}{h_{L}h_{p}}\right)^{\alpha}\right)$
$\displaystyle~{}~{}~{}\times\left(\left(\frac{h_{p}}{G_{0}}\right)^{1/\beta_{t}-1}-\left(\frac{h_{p}}{G_{0}}\right)^{1/\beta_{r}-1}\right)\text{d}h_{p}$
(85)

Applying a change of variable rule
$x=\left(\frac{h_{p}}{G_{0}}\right)^{-\alpha}$, Eq. (E) is rewritten as

$\displaystyle
f_{h}(h)=\frac{1}{\alpha(\beta_{t}-\beta_{r})}\frac{1}{h_{L}G_{0}}A_{1}\left(\mathbb{A}_{0}(h)\right)^{\mu-\frac{1}{\alpha}}\times$
$\displaystyle\int_{1}^{\infty}\exp\left(-A_{2}\mathbb{A}_{0}(h)x\right)\left(x^{-1/\alpha\beta_{t}+\mu-1}-x^{-1/\alpha\beta_{r}+\mu-1}\right)\text{d}x$
(86)

In the following derivation, we use an integral identity [21]

$\displaystyle\int_{u}^{\infty}x^{-\nu}{e^{-x}}\text{d}x=u^{-\nu/2}e^{-u/2}\mathbb{W}_{-\frac{\nu}{2},\frac{1-\nu}{2}}(u)$
(87)

Using Eqs. (E) and (87), the closed-form expressions for the PDF of $h$ is
derived in Eq. (IV). Note that the results of Eq. (IV) are based on the
Whittaker function. If
$\big{(}\alpha\mu\min\\{\beta_{t},\beta_{r}\\}\big{)}>1$, we can express Eq.
(IV) based on the upper incomplete gamma function $\Gamma(s,x)$, which has
faster computation than the Whittaker function. Using Eq. (E) and [22], and
after some derivation, the closed-form expressions for the PDF of $h$ is
derived in Eq. (IV).

Substituting Eqs. (IV) and (81) in Eq. (80), we obtain:

$\displaystyle
F_{h}(h)=1-\frac{1}{G_{0}(\beta_{t}-\beta_{r})}\sum_{k=0}^{\mu-1}\frac{A_{2}^{k}}{\Gamma(k+1)}\int_{0}^{G_{0}}\left(\frac{h}{h_{L}h_{p}}\right)^{k\alpha}\times$
$\displaystyle\exp\left(-A_{2}\left(\frac{h}{h_{L}h_{p}}\right)^{\alpha}\right)\\!\left[\\!\left(\frac{h_{p}}{G_{0}}\right)^{1/\beta_{t}-1}-\left(\frac{h_{p}}{G_{0}}\right)^{1/\beta_{r}-1}\right]\text{d}h_{p}$
(88)

Applying a change of variable rule
$x=\left(\frac{h_{p}}{G_{0}}\right)^{-\alpha}$ and using Eqs. (87) and (E),
the closed-form expressions for the CDF of $h$ is derived in Eq. (35).

## Appendix F

In this section, we derive the PDF and CDF of $h$ for the case in which the Tx
and Rx are equipped with ULA antenna. Substituting Eqs. (3) and (IV) in Eq.
(75), we obtain:

$\displaystyle f_{h}(h)$
$\displaystyle=\sum_{n=1}^{N}\frac{A_{1}\Re_{4}\Re_{5}}{\Re_{3}(n)(\beta_{ty}+\beta_{ry})G_{0}}\int_{0}^{G_{0}}\frac{1}{h_{L}h_{p}}\left(\frac{h}{h_{L}h_{p}}\right)^{\alpha\mu-1}$
$\displaystyle\times\exp\left(-A_{2}\left(\frac{h}{h_{L}h_{p}}\right)^{\alpha}\right)\left(\frac{h_{p}}{G_{0}}\right)^{\frac{\Re_{4}}{\beta_{ty}+\beta_{ry}}-1}\text{d}h_{p}$
(89)

Also, substituting Eqs. (IV) and (81) in Eq. (80), we obtain:

$\displaystyle
F_{h}(h)=1-\sum_{k=0}^{\mu-1}\sum_{n=1}^{N}\frac{\Re_{4}\Re_{5}}{\Re_{3}(n)(\beta_{ty}+\beta_{ry})}\frac{A_{2}^{k}\left(\frac{h}{G_{0}h_{L}}\right)^{\alpha
k}}{\Gamma(k+1)}\times$
$\displaystyle\int_{0}^{1}\exp\left(-A_{2}\left(\frac{h}{G_{0}h_{L}}\right)^{\alpha}y^{-\alpha}\right)y^{-k\alpha+\frac{\Re_{4}}{\beta_{ty}+\beta_{ry}}-1}\text{d}y$
(90)

Using Eqs. (F) and (F) and similar to the method used in Appendix E, the
closed-form expressions for the PDF and CDF of $h$ are derived in Eqs. (IV)
and (IV), respectively.

## References

  * [1] D. Moltchanov, E. Sopin, V. Begishev, A. Samuylov, Y. Koucheryavy, and K. Samouylov, “A tutorial on mathematical modeling of 5G/6G millimeter wave and terahertz cellular systems,” _IEEE Communications Surveys & Tutorials_, 2022.
  * [2] V. Petrov, D. Moltchanov, Y. Koucheryavy, and J. M. Jornet, “Capacity and outage of terahertz communications with user micro-mobility and beam misalignment,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 6, pp. 6822–6827, 2020.
  * [3] H. Chen, H. Sarieddeen, T. Ballal, H. Wymeersch, M.-S. Alouini, and T. Y. Al-Naffouri, “A tutorial on terahertz-band localization for 6G communication systems,” _IEEE Communications Surveys & Tutorials_, 2022\.
  * [4] H. Sarieddeen, M.-S. Alouini, and T. Y. Al-Naffouri, “An overview of signal processing techniques for terahertz communications,” _Proceedings of the IEEE_ , 2021.
  * [5] O. S. Badarneh, “Performance analysis of tera-hertz communications in random fog conditions with misalignment,” _IEEE Wireless Commun. Let._ , 2022.
  * [6] P. Bhardwaj and S. M. Zafaruddin, “Performance of dual-hop relaying for THz-RF wireless link over asymmetrical $\alpha$-$\mu$ fading,” _IEEE Transactions on Vehicular Technology_ , vol. 70, no. 10, pp. 10 031–10 047, 2021\.
  * [7] H. Du, J. Zhang, K. Guan, D. Niyato, H. Jiao, Z. Wang, and T. Kürner, “Performance and optimization of reconfigurable intelligent surface aided THz communications,” _IEEE Transactions on Communications_ , 2022.
  * [8] A.-A. A. Boulogeorgos, E. N. Papasotiriou, and A. Alexiou, “Analytical performance assessment of THz wireless systems,” _IEEE Access_ , vol. 7, pp. 11 436–11 453, 2019.
  * [9] X. Tong, B. Chang, Z. Meng, G. Zhao, and Z. Chen, “Calculating terahertz channel capacity under beam misalignment and user mobility,” _IEEE Wireless Communications Letters_ , 2021.
  * [10] A.-A. A. Boulogeorgos, A. Alexiou, and M. Di Renzo, “Outage performance analysis of RIS-assisted UAV wireless systems under disorientation and misalignment,” _arXiv preprint arXiv:2201.12056_ , 2022.
  * [11] A.-A. A. Boulogeorgos, J. M. Riera, and A. Alexiou, “On the joint effect of rain and beam misalignment in terahertz wireless systems,” _IEEE Access_ , 2022.
  * [12] S. Liu, X. Yu, R. Guo, Y. Tang, and Z. Zhao, “THz channel modeling: Consolidating the road to THz communications,” _China Communications_ , vol. 18, no. 5, pp. 33–49, 2021.
  * [13] M. T. Dabiri and M. Hasna, “Pointing error modeling of mmWave to THz high-directional antennas,” _arXiv preprint arXiv:2206.10756_ , 2022.
  * [14] A. A. Farid and S. Hranilovic, “Outage capacity optimization for free-space optical links with pointing errors,” _J. Lightw. technol._ , vol. 25, no. 7, pp. 1702–1710, 2007.
  * [15] M. T. Dabiri, S. M. S. Sadough, and M. A. Khalighi, “Channel modeling and parameter optimization for hovering UAV-based free-space optical links,” _IEEE J. Sel. Areas Commun._ , vol. 36, no. 9, pp. 2104–2113, 2018.
  * [16] J. Kokkoniemi, A.-A. A. Boulogeorgos, M. Aminu, J. Lehtomäki, A. Alexiou, and M. Juntti, “Impact of beam misalignment on THz wireless systems,” _Nano Commun. Netw._ , vol. 24, p. 100302, 2020.
  * [17] C. A. Balanis, _Antenna theory: analysis and design_. John wiley & sons, 2016.
  * [18] S. Tripathi, N. V. Sabu, A. K. Gupta, and H. S. Dhillon, “Millimeter-wave and terahertz spectrum for 6G wireless,” in _6G Mobile Wireless Netw._ Springer, 2021, pp. 83–121.
  * [19] E. N. Papasotiriou, A.-A. A. Boulogeorgos, K. Haneda, M. F. de Guzman, and A. Alexiou, “An experimentally validated fading model for THz wireless systems,” _Scientific Reports_ , vol. 11, no. 1, pp. 1–14, 2021.
  * [20] E. W. Weisstein, “Marcum q-function.” From MathWorld–A Wolfram Web Resource. [Online]. Available: https://mathworld.wolfram.com/MarcumQ-Function.html
  * [21] I. S. Gradshteyn and I. M. Ryzhik, _Table of integrals, series, and products_. 7th ed. Academic press, 2007\.
  * [22] E. W. Weisstein, “Incomplete gamma function.” From MathWorld–A Wolfram Web Resource.
  * [23] A. Papoulis and S. U. Pillai, _Probability, random variables, and stochastic processes_. Tata McGraw-Hill Education, 2002.
  * [24] P. G. Moschopoulos, “The distribution of the sum of independent gamma random variables,” _Annals of the Institute of Statistical Mathematics_ , vol. 37, no. 1, pp. 541–544, 1985.
  * [25] M. Nakagami, “The m-distribution—A general formula of intensity distribution of rapid fading,” in _Statistical methods in radio wave propagation_. Elsevier, 1960, pp. 3–36.
  * [26] (2016) Hoytdistribution. Wolfram Research. [Online]. Available: https://reference.wolfram.com/language/ref/HoytDistribution.html
  * [27] G. Tavares, “Efficient computation of Hoyt cumulative distribution function,” _Electronics letters_ , vol. 46, no. 7, pp. 537–539, 2010.