# Old recurrence formulae for growth series of Coxeter groups

Jan Dymara (Wrocław)  Instytut Matematyczny Uniwersytet Wrocławski pl.
Grunwaldzki 2/4, 50-384 Wrocław Poland<EMAIL_ADDRESS>
###### Abstract.

Several classical formulae for the growth series of a Coxeter group are proved
in a new way, using the structure of the Coxeter complex, the Davis complex,
or the Tits non-complex.

###### Key words and phrases:

Coxeter group, growth series

###### 2010 Mathematics Subject Classification:

Primary 20F55; Secondary 57M07

##  1\. Introduction

Throughout this paper $(W,S)$ is a Coxeter system. This means that $W$ is the
group generated by the finite set $S$, subject to relations of the form
$(st)^{m_{st}}=1$, where $m_{ss}=1$ and for $t\neq s$ we have
$m_{st}=m_{ts}\in\\{2,3,\ldots,\infty\\}$ ($\infty$ means no relation). We
will usually say “Coxeter group $W$” instead of mentioning the whole system.
We denote by $\ell(w)$ the word length of an element $w\in W$ with respect to
the generating set $S$. For any subset $T\subseteq S$ the subgroup $W_{T}$ of
$W$ generated by $T$ is also a Coxeter group, with relations being the
relevant relations of $W$. The word length in $W_{T}$ agrees with the
restriction to $W_{T}$ of the word length of $W$; we denote both by $\ell$.
The growth series of $W$ is the formal power series $W(t)=\sum_{w\in
W}t^{\ell(w)}$. We will prove several well-known formulae for this function.
For infinite $W$:

(1) $\sum_{T\subseteq S}\frac{(-1)^{|T|}}{W_{T}(t)}=0.$

If $W$ is finite, it has a unique element of longest length $m$. Then

(2) $\sum_{T\subseteq S}\frac{(-1)^{|T|}}{W_{T}(t)}=\frac{t^{m}}{W(t)}.$

A subset $T\subseteq S$ is called spherical if $W_{T}$ is finite. We denote by
${\mathcal{S}}$ the set of all spherical subsets of $S$. Let
$\chi_{T}=\sum_{T\subseteq U\in{\mathcal{S}}}(-1)^{|U|}$. Then

(3)
$\sum_{T\in{\mathcal{S}}}\frac{(-1)^{|T|}\chi_{T}}{W_{T}(t)}=\frac{1}{W(t)};$
(4) $\sum_{T\in{\mathcal{S}}}\frac{(-1)^{|T|}}{W_{T}(t)}=\frac{1}{W(t^{-1})}.$

These formulae are classical and quite popular, cf. [Bou, Ex. IV.26], [D,
Chapter Seventeen], [ChD]; variants of Coxeter group growth series are
investigated even nowadays, cf. [OS], [BGM, Prop. 8.3]. The usual proofs use
an inclusion–exclusion principle for appropriate subsets of $W$, or induction
on the cardinality of $S$. Our goal is to use spaces on which $W$ acts to
establish (1), (2), (3) and (4) in an almost uniform way.

##  2\. Some properties of Coxeter groups

The basic reference is [Bou], a newer one [D].

  1. $\bullet$

For any $T\subseteq S$, and any $w\in W$, the coset $wW_{T}$ has a unique
shortest element, say $u$. Moreover, if $w=uw^{\prime}$, then
$\ell(w)=\ell(u)+\ell(w^{\prime})$. ([D, Lemma 4.3.1])

  2. $\bullet$

If $W$ is finite, then it has a unique element $w_{0}$ of longest length
$m=\ell(w_{0})$. Then, for any $w\in W$ we have
$\ell(w)+\ell(w^{-1}w_{0})=\ell(w_{0})$ ([D, Lemma 4.6.1]). Consequently, in
$W$ there are as many elements of length $k$ as of length $m-k$. More
succinctly, $W(t)=t^{m}W(t^{-1})$.

  3. $\bullet$

For any $w\in W$, the set $\mathop{\rm In}(w)=\\{s\in
S\mid\ell(ws)<\ell(w)\\}$ is spherical ([D, Lemma 4.7.2]).

##  3\. Basic spaces

The basic construction ([D], Chapter Five) is performed as follows. Let $Y$ be
a space (topological, or a simplicial complex) with a collection of subspaces
$(Y_{s})_{s\in S}$ (called panels). For $y\in Y$ we put $S(y)=\\{s\in S\mid
y\in Y_{s}\\}$. Then $\mathop{\mathcal{U}}(W,Y)=W\times Y/\sim$, where
$(w,y)\sim(w^{\prime},y^{\prime})\iff y=y^{\prime}$ and $w^{-1}w^{\prime}\in
W_{S(y)}$. The image of $\\{w\\}\times Y$ in $\mathop{\mathcal{U}}(W,Y)$ will
be called a chamber and denoted $wY$. The set of all chambers, $\mathop{\rm
Ch}(\mathop{\mathcal{U}}(W,Y))$, is in bijective correspondence with $W$,
provided $Y\neq\bigcup_{s\in S}Y_{s}$. Via this bijection the length function
is transported to the set of chambers: $\ell(wY)=\ell(w)$. The $W$-action on
$\mathop{\mathcal{U}}(W,Y)$ is defined by $w(u,y)=(wu,y)$. For $T\subseteq S$
we put $Y_{T}=\bigcap_{t\in T}Y_{t}$ (we call it the face of type $T$) and
$Y^{T}=\bigcup_{t\in T}Y_{t}$; this notation and terminology equivariantly
extends to all chambers. We will use the notation $\sigma<X$ to indicate that
$\sigma$ is a simplex of a simplicial complex $X$.

We will use three spaces arising from this construction.

1\. The Coxeter complex $X_{\Delta}=\mathop{\mathcal{U}}(W,\Delta)$, where
$\Delta$ is an $(|S|-1)$-dimensional simplex, and $(\Delta_{s})_{s\in S}$ is
the collection of its codimension-1 faces.

2\. The Davis complex $X_{D}=\mathop{\mathcal{U}}(W,D)$. The chamber model $D$
(the Davis chamber) is the subcomplex of the first barycentric subdivision
$\Delta^{\prime}$ of $\Delta$ spanned by the barycentres of faces of $\Delta$
that are of spherical type (for finite $W$ one should include also the
barycentre of the empty face—we shall not consider this case). Alternatively,
it is the subcomplex $X_{D}$ of the first barycentric subdivision
$X^{\prime}_{\Delta}$ spanned by barycentres of simplices of $X_{\Delta}$ that
are contained in finitely many chambers. It is a locally finite complex. The
panels are $D_{s}=D\cap\Delta_{s}$, and the faces $D_{T}=D\cap\Delta_{T}$. An
important property: for spherical $T$ the face $D_{T}$ is contractible;
indeed, it is a cone with apex at the barycentre of $\Delta_{T}$. For non-
spherical $T$ we have $D_{T}=\emptyset$. For any simplex $\sigma<D$ we define
its type $S(\sigma)=\\{s\in S\mid\sigma<D_{s}\\}$; this definition extends by
invariance to all simplices in $X_{D}$.

3\. The Tits non-complex $X_{\Delta}^{f}=\mathop{\mathcal{U}}(W,\Delta^{f})$.
Here $\Delta^{f}$ is obtained from $\Delta$ by removing the closed simplices
whose types are non-spherical. It is not a simplicial complex. For faces
$\sigma$ of $\Delta$ of spherical type we still use the notation
$\sigma<\Delta^{f}$, even though, in the geometric realisation, some part of
$\sigma$ may not be a subset of $\Delta^{f}$.

For infinite $W$, all of the above spaces are contractible; $X_{D}$ and
$X_{\Delta}^{f}$ even carry good CAT(0) metrics. Though they are generally
very useful, these properties will not be relied upon in this note.

##  4\. Euler characteristic generating function

Let $X$ be either the Coxeter complex or the Davis complex for the Coxeter
system $(W,S)$. (The case of the Tits non-complex is slightly different and
will be dealt with in the last section.) For a simplex $\sigma<X$ we define
its length $\ell(\sigma)$ as the minimum of $\ell(C)$, where $C$ runs through
the chambers of $X$ that contain $\sigma$. There are finitely many chambers
$C$ of $X$ of any given length, and each of them contains finitely many
simplices. Therefore the following formal power series—the Euler
characteristic generating function of $X$—is well defined:

(5) $\chi_{t}(X)=\sum_{\sigma<X}(-1)^{\dim{\sigma}}t^{\ell(\sigma)}.$

We will calculate this series in two ways: grouping simplices according to
type, and by local summation in $X$. Comparing the results we will obtain the
formulae (1),(2) and (3).

##  5\. Summation according to type

We consider $X=\mathop{\mathcal{U}}(W,Y)$ for $Y=\Delta$ or $D$. Let
$\sigma<X$ be a simplex of type $T$. Then the chambers $C$ that contain
$\sigma$ correspond to elements of some coset $uW_{T}$. Assume that $u$ is the
shortest element in that coset. Then $\ell(\sigma)=\ell(u)$, and

(6) $\sum_{C>\sigma}t^{\ell(C)}=\sum_{w\in uW_{T}}t^{\ell(w)}=\sum_{v\in
W_{T}}t^{\ell(u)}t^{\ell(v)}=t^{\ell(u)}W_{T}(t)=t^{\ell(\sigma)}W_{T}(t).$

Multiplying both sides by $(-1)^{\dim{\sigma}}$ and summing over all simplices
of type $T$ we get

(7) $\sum_{C\in\mathop{\rm
Ch}(X)}t^{\ell(C)}\sum_{\genfrac{}{}{0.0pt}{}{\sigma<C}{S(\sigma)=T}}(-1)^{\dim{\sigma}}=W_{T}(t)\sum_{\genfrac{}{}{0.0pt}{}{\sigma<X}{S(\sigma)=T}}(-1)^{\dim\sigma}t^{\ell(\sigma)},$

which implies

(8)
$\frac{W(t)}{W_{T}(t)}\sum_{\genfrac{}{}{0.0pt}{}{\sigma<Y}{S(\sigma)=T}}(-1)^{\dim{\sigma}}=\sum_{\genfrac{}{}{0.0pt}{}{\sigma<X}{S(\sigma)=T}}(-1)^{\dim\sigma}t^{\ell(\sigma)}.$

Finally, we sum over all $T\subseteq S$:

(9) $W(t)\sum_{T\subseteq
S}\frac{1}{W_{T}(t)}\sum_{\genfrac{}{}{0.0pt}{}{\sigma<Y}{S(\sigma)=T}}(-1)^{\dim{\sigma}}=\chi_{t}(X).$

We will calculate the inner sum separately in each case.

The Coxeter complex. For any proper subset $T\subset S$ there is one
$\sigma<\Delta$ of type $T$, namely $\Delta_{T}$. It has dimension
${|S|-|T|-1}$. Thus

(10) $\chi_{t}(X_{\Delta})=W(t)\sum_{T\subset
S}\frac{(-1)^{|S|-|T|-1}}{W_{T}(t)}.$

The Davis complex. Recall that the face $D_{T}$ is a cone with apex at the
barycentre of $\Delta_{T}$; let $L_{T}^{\prime}$ be its base
$D_{T}\cap\partial\Delta_{T}$. The simplices $\sigma<D$ of type $T$ are the
interior simplices of this cone (including the apex), so that

(11)
$\sum_{\genfrac{}{}{0.0pt}{}{\sigma<D}{S(\sigma)=T}}(-1)^{\dim{\sigma}}=1-\chi(L_{T}^{\prime}).$

Simplices in $L_{T}^{\prime}$ correspond to chains $U_{1}\subset\ldots\subset
U_{k}$, where $U_{i}$ are spherical and properly contain $T$. This means that
$L_{T}^{\prime}$ is, as the notation suggests, the barycentric subdivision of
the following simplicial complex $L_{T}$. The vertices of $L_{T}$ correspond
to spherical sets of the form $T\cup\\{u\\}$ (with $u\in S-T$); vertices
$T\cup\\{u_{1}\\},\ldots,T\cup\\{u_{k}\\}$ span a simplex if
$T\cup\\{u_{1},\ldots,u_{k}\\}$ is spherical. Thus, simplices of $L_{T}$
correspond to spherical $U$ containing $T$, and

(12) $1-\chi(L_{T}^{\prime})=1-\chi(L_{T})=\sum_{T\subseteq
U\in{\mathcal{S}}}(-1)^{|U|-|T|}=(-1)^{|T|}\chi_{T}$

for $\chi_{T}=\sum_{T\subseteq U\in{\mathcal{S}}}(-1)^{|U|}$. Finally,
returning with this to (9) we get

(13)
$\chi_{t}(X_{D})=W(t)\sum_{T\in{\mathcal{S}}}\frac{(-1)^{|T|}\chi_{T}}{W_{T}(t)}.$

##  6\. Local summation

We continue with the cases $Y=\Delta$ or $D$. Now we compute the sum
$\sum_{\sigma<X}(-1)^{\dim{\sigma}}t^{\ell(\sigma)}$ by arranging its summands
according to chambers: we consider each $\sigma$ as part of the chamber
realising $\ell(\sigma)$. We get

(14) $\chi_{t}(X)=\sum_{C\in\mathop{\rm
Ch}(X)}t^{\ell(C)}\sum_{\genfrac{}{}{0.0pt}{}{\sigma<C}{\ell(\sigma)=\ell(C)}}(-1)^{\dim{\sigma}}=\sum_{w\in
W}t^{\ell(w)}\sum_{\genfrac{}{}{0.0pt}{}{\sigma<wY}{\ell(\sigma)=\ell(w)}}(-1)^{\dim{\sigma}}.$

######  Lemma 6.1.

Let $\sigma<wY$. Then $\ell(\sigma)<\ell(w)$ if and only if
$\sigma<(wY)^{\mathop{\rm In}(w)}$.

###### Proof.

Let $T$ be the type of $\sigma$. The chambers containing $\sigma$ correspond
to elements of the coset $wW_{T}$.

Suppose $\ell(\sigma)<\ell(w)$. It follows that $w$ is not the shortest
element of the coset $wW_{T}$. Therefore $\ell(wt)<\ell(w)$ for some $t\in T$.
This $t$ also belongs to $\mathop{\rm In}(w)$. Finally, $\sigma<wY\cap
wtY=(wY)_{t}$.

Conversely, if $\sigma<(wY)^{\mathop{\rm In}(w)}$, then $\sigma<(wY)_{s}$ for
some $s\in\mathop{\rm In}(w)$. But then $\sigma<wsY$ and
$\ell(\sigma)\leq\ell(ws)<\ell(w)$. The Lemma is proved. ∎

By the Lemma, the inner sum in (14) is performed over simplices of $wY$ that
are not in $(wY)^{\mathop{\rm In}(w)}$. Thus

(15)
$\displaystyle\sum_{\genfrac{}{}{0.0pt}{}{\sigma<wY}{\ell(\sigma)=\ell(w)}}(-1)^{\dim{\sigma}}$
$\displaystyle=\chi(wY,(wY)^{\mathop{\rm In}(w)})=\chi(Y,Y^{\mathop{\rm
In}(w)})$ $\displaystyle=\chi(Y)-\chi(Y^{\mathop{\rm
In}(w)})=1-\chi(Y^{\mathop{\rm In}(w)}),$

the last equality due to contractibility of $Y$.

We split the calculation of $\chi(Y^{\mathop{\rm In}(w)})$ into three cases.

1) The generic case: $\mathop{\rm In}(w)$ is a proper non-empty subset of $S$,
of some cardinality $k$. Notice that for any $T\subseteq\mathop{\rm In}(w)$ (a
fortiori, spherical) the face $Y_{T}$ is contractible. Using the
inclusion–exclusion principle we get

(16) $\displaystyle\chi($ $\displaystyle Y^{\mathop{\rm
In}(w)})=\sum_{s\in\mathop{\rm
In}(w)}\chi(Y_{s})-\sum_{\genfrac{}{}{0.0pt}{}{s,t\in\mathop{\rm In}(w)}{s\neq
t}}\chi(Y_{\\{s,t\\}})+\ldots+(-1)^{k-1}\chi(Y_{\mathop{\rm In}(w)})$
$\displaystyle=k-\binom{k}{2}+\binom{k}{3}-\ldots+(-1)^{k-1}\binom{k}{k}=1-(1-1)^{k}=1.$

2) The set $\mathop{\rm In}(w)$ is empty. Then $\chi(Y^{\mathop{\rm
In}(w)})=\chi(\emptyset)=0$. This happens exactly when $w=1$.

3) The case $\mathop{\rm In}(w)=S$. This happens only if $W$ is finite and $w$
is the element of longest length (say length $m$). Then $Y^{S}$ is a
triangulation of the $(|S|-2)$-dimensional sphere,
$\chi(Y^{S})=1-(-1)^{|S|-1}$.

Plugging these results into (15), and then into (14), we get for infinite $W$

(17) $\chi_{t}(X_{\Delta})=1,\qquad\chi_{t}(X_{D})=1,$

while for finite $W$ with element of longest length $m$

(18) $\chi_{t}(X_{\Delta})=1+(-1)^{|S|-1}t^{m}.$

Comparing these results with (10) and (13) we get (1), (2) and (3).

##  7\. The Tits non-complex

Consider now the case $X=X_{\Delta}^{f}$. For $\sigma<X$ we define $L(\sigma)$
as the maximum of $\ell(C)$ over all chambers containing $\sigma$. This is
finite, since the type of $\sigma$ is spherical. Then we define
$\chi^{t}(X)=\sum_{\sigma<X}(-1)^{\dim{\sigma}}t^{L(\sigma)}$. As before, we
will calculate this sum in two ways.

First, we group the summands according to the type of $\sigma$. Let $\sigma<X$
be a simplex of (spherical!) type $T$. Then the chambers $C$ that contain
$\sigma$ correspond to elements of some finite coset $uW_{T}$. Assume that $u$
is the longest element in that coset. Then $L(\sigma)=\ell(u)$, and

(19) $\displaystyle\sum_{C>\sigma}t^{\ell(C)}=\sum_{w\in
uW_{T}}t^{\ell(w)}=\sum_{v\in W_{T}}t^{\ell(u)}t^{-\ell(v)}$
$\displaystyle=t^{\ell(u)}W_{T}(t^{-1})$
$\displaystyle=t^{L(\sigma)}W_{T}(t^{-1}).$

Multiplying both sides by $(-1)^{\dim{\sigma}}$($=(-1)^{|S|-|T|-1}$) and
summing over all simplices of type $T$ (there is just one in each chamber) we
get

(20) $(-1)^{|S|-|T|-1}\sum_{C\in\mathop{\rm
Ch}(X)}t^{\ell(C)}=W_{T}(t^{-1})\sum_{\genfrac{}{}{0.0pt}{}{\sigma<X}{S(\sigma)=T}}(-1)^{\dim\sigma}t^{L(\sigma)},$

which implies

(21)
$\frac{W(t)}{W_{T}(t^{-1})}(-1)^{|S|-|T|-1}=\sum_{\genfrac{}{}{0.0pt}{}{\sigma<X}{S(\sigma)=T}}(-1)^{\dim\sigma}t^{L(\sigma)}.$

Summing this equality over all spherical types $T$ we get

(22)
$W(t)\sum_{T\in{\mathcal{S}}}\frac{(-1)^{|S|-|T|-1}}{W_{T}(t^{-1})}=\chi^{t}(X^{f}_{\Delta}).$

Second, we group summands in chambers according to $L$:

(23) $\chi^{t}(X)=\sum_{\sigma<X}(-1)^{\dim{\sigma}}t^{L(\sigma)}=\sum_{w\in
W}t^{\ell(w)}\sum_{\genfrac{}{}{0.0pt}{}{\sigma<w\Delta^{f}}{L(\sigma)=\ell(w)}}(-1)^{\dim{\sigma}}.$

Every $\sigma<w\Delta^{f}$ is of the form $w\Delta_{T}$ for some spherical
$T$; moreover, $L(\sigma)=\ell(w)$ if and only if $w$ is the longest element
in the coset $wW_{T}$, which happens exactly when $T\subseteq\mathop{\rm
In}(w)$. Therefore

(24)
$\sum_{\genfrac{}{}{0.0pt}{}{\sigma<w\Delta^{f}}{L(\sigma)=\ell(w)}}(-1)^{\dim{\sigma}}=\sum_{T\subseteq\mathop{\rm
In}(w)}(-1)^{|S|-|T|-1}=(-1)^{|S|-1}(1-1)^{|\mathop{\rm In}(w)|}.$

Thus, the only non-zero summand corresponds to $w$ with empty $\mathop{\rm
In}(w)$, that is to $w=1$. We get

(25) $\chi^{t}(X^{f}_{\Delta})=(-1)^{|S|-1}.$

Comparing (25) with (22), and switching $t$ and $t^{-1}$, we get (4).

### Acknowledgements

The author was partially supported by the Polish National Science Centre (NCN)
grant 2016/23/B/ST1/01556.

## References

  * [1]
  * [Bou] N. Bourbaki. _Groupes et algèbres de Lie, Chapitres IV–VI._ Masson, 2nd edition, 1981.
  * [BGM] M. Bożejko, Ś. Gal, W. Młotkowski. _Positive definite functions on Coxeter groups with applications to operator spaces and noncommutative probability_ , Comm. Math. Phys. 361(2018), no. 2, 583–604.
  * [ChD] R. Charney, M. Davis. _Reciprocity of growth functions of Coxeter groups_. Geom. Dedicata 39(1991), no. 3, 373–378.
  * [D] M. Davis. _The Geometry and Topology of Coxeter Groups_. Princeton University Press, Princeton and Oxford 2008.
  * [OS] B. Okun, R. Scott. _Growth series of CAT(0) cubical complexes_ , Topology Proc. 54(2019), 295–303.


# On the maximum of random assignment process

M.A.Lifshitsa,∗, A.A.Tadevosianb

( aSaint-Petersburg State University, University Emb. 7/9, 199034, St
Petersburg, Russia<EMAIL_ADDRESS>
bSaint-Petersburg State University, University Emb. 7/9, 199034, St
Petersburg, Russia<EMAIL_ADDRESS>
∗Corresponding author
)

###### Abstract

We describe the behavior of the expectation of the maximum for a random
assignment process built upon a square matrix with independent entries. Under
mild assumptions on the underlying distribution, the answer is expressed in
terms of its quantile function.

Keywords: Assignment problem, Random assignment

##  1 Introduction

We consider the following random assignment problem. Let ($X_{ij}$) be an
$n\times n$ random matrix with i.i.d. random entries having a common
distribution $\mathcal{P}$. Let $\SS_{n}$ denote the group of permutations
$\pi:\\{1,2,\dots,n\\}\to\\{1,2,\dots,n\\}$. For every $\pi\in\SS_{n}$ let

$S(\pi)=\sum\limits_{i=1}^{n}X_{i\pi(i)}.$

We call such process $\\{S(\pi),\,\pi\in\SS_{n}\\}$ a random assignment
process. We are interested in the study of asymptotic behaviour of
${\mathbb{E}}\max\limits_{\pi\in\SS_{n}}S(\pi)$ as $n\to\infty$.

We refer to [8, 18] for many applications of assignment problem in various
fields of mathematics.

The setting with $(X_{ij})$ uniformly distributed on $[0,1]$ was studied by
Steele [18] and Mézard and Parisi [13], where the authors proved that

${\mathbb{E}}{\min\limits_{\pi\in\SS_{n}}S(\pi)}=\zeta(2)-\frac{\zeta(2)/2+2\zeta(3)}{n}+O\left(\frac{1}{n^{2}}\right),\text{
as }n\to\infty,$

$\zeta(\cdot)$ being Riemann’s zeta function. Mézard et al. [14] also
conjectured that in the exponential case when $\mathcal{P}=\mathrm{Exp}(1)$ it
is true that

${\mathbb{E}}{\min\limits_{\pi\in\SS_{n}}S(\pi)}\to\zeta(2),\text{ as
}n\to\infty.$

Using replica method from statistical physics [9], they provided an
heuristical argumentation in favor of this conjecture.

Later Parisi [17] conjectured the following explicit expression for every
fixed $n$

${\mathbb{E}}{\min\limits_{\pi\in\SS_{n}}S(\pi)}=\sum\limits_{k=1}^{n}\frac{1}{k^{2}}.$

It is compatible with exact results for $n=1,2$ and for $\zeta(2)$ limit
conjecture as $n\to\infty$. This conjecture was confirmed by numeric
simulations for $n=3,4,5$. Special cases and generalizations are also studied
in [4, 6, 10]. General version of formula above was proved independently by
Nair et. al [16] and by Linusson and Wästlund [12]. Later, Wästlund proposed a
simpler proof [19].

Aldous [2] gave a rigorous proof of $\zeta(2)$ limit in the exponential case.
His approach is based on the assignment analysis of an infinite tree with
edges provided with exponentially distributed weights, see [3]. In the same
work he showed that the distribution $\mathcal{P}$ affects
$\lim_{n}{\mathbb{E}}\,{\min_{\pi}S(\pi)}$ only through the value of its
probability density function at 0, thus results for $\mathcal{P}=U[0,1]$ and
$\mathcal{P}=\mathrm{Exp}(1)$ are equivalent in asymptotic sense.

In [4, 6, 8, 12, 16, 19] similar problems were investigated for rectangular
matrices.

The situation is very different when one deals with the variables having
unbounded distributions ${\mathcal{P}}$. For the study of minima described
above this would mean that ${\mathbb{P}}(X_{ij}<r)>0$ for every real $r$. In
the following we find more natural to switch from minima to maxima and to
assume that ${\mathbb{P}}(X_{ij}>r)>0$ for every real $r$. For such case the
expectation of maxima does not tend anymore to a finite limit but grows to
infinity and the problem is to evaluate the order of this growth.

In this direction, Mordant and Segers [15] showed that in Gaussian case
$\mathcal{P}=\mathcal{N}(0,1)$ it is true that

${\mathbb{E}}{\max\limits_{\pi\in\SS_{n}}S(\pi)}=n\sqrt{2\log{n}}(1+o(1)),$

and proposed a greedy method giving random permutation at which the
permutation process attains the value asymptotically equivalent to the
expectation of maximum. Essentially the same results for Gaussian case were
obtained independently by the authors [11] who also used the greedy method but
the priority of [15] in its preprint version must be of course acknowledged
here. In their preprint, Cheng et al. [7] considered assignment processes on
more general graph structures and more general random variables satisfying
some special assumptions.

In this paper our goal is to provide some mild and natural conditions imposed
on the distribution $\mathcal{P}$ in order to obtain asymptotic expression for
expected maximum of the random assignment process.

For notation convenience let $(X_{i}),X$ be some auxiliary i.i.d. random
variables having the same distribution $\mathcal{P}$ and denote
$M_{n}=\max_{1\leq i\leq n}X_{i}$ the maximum of i.i.d. random variables. Also
let us denote $\mathcal{M}_{n}=\max_{\pi\in\SS_{n}}S(\pi)$ the maximum of the
assignment process built on matrix $(X_{ij})$. Let us define the quantile
function by

$g(p):=\inf\\{r:{\mathbb{P}}(X\geq r)\leq p\\},\qquad 0<p<1.$

Our main result is the following theorem:

###### Theorem 1

Assume that

1) ${\mathbb{E}}|X|<\infty$.

2) At zero the function $g(\cdot)$ tends to infinity and is slowly varying.

Then

${\mathbb{E}}M_{n}=g(1/n)\,(1+o(1)),\qquad\textrm{as }n\to\infty,$ (1)

and

${\mathbb{E}}\mathcal{M}_{n}=n\,g(1/n)\,(1+o(1)),\qquad\textrm{as
}n\to\infty.$ (2)

The structure of the paper is as follows: In Section 2 we study the maxima of
an i.i.d. sequence and prove (1). Then, in Section 3 we establish bilateral
relations between the maximum of the assignment process and that of the
corresponding i.i.d. sequence in order to derive (2) from (1). Finally, we
provide expectation asymptotics of the assignment process maximum for some
standard probability distributions.

##  2 Maximum’s expectation in i.i.d. case

We split the proof of (1) into two parts – a lower bound and an upper bound.

###  2.1 Lower bound

###### Proposition 2

Under theorem’s assumptions, it is true that

${\mathbb{E}}M_{n}\geq g(1/n)\,(1+o(1)),\quad\quad as\,\,n\to\infty.$

Proof: For every $R>0$ we have

$\displaystyle{\mathbb{E}}M_{n}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}{\mathbb{P}}(M_{n}\geq
r)dr-\int_{0}^{\infty}{\mathbb{P}}(M_{n}\leq-r)dr$ (3) $\displaystyle\geq$
$\displaystyle R\,{\mathbb{P}}(M_{n}\geq
R)-\int_{0}^{\infty}{\mathbb{P}}(M_{n}\leq-r)dr.$

Let us fix some $V>0$ and let $R=R_{n}:=g(V/n)$. Then ${\mathbb{P}}(X\geq
R)\geq\frac{V}{n}$ and

$\displaystyle{\mathbb{P}}(M_{n}\geq R)$ $\displaystyle=$ $\displaystyle
1-{\mathbb{P}}(M_{n}<R)=1-{\mathbb{P}}(X<R)^{n}$ $\displaystyle=$
$\displaystyle 1-(1-{\mathbb{P}}(X\geq R))^{n}\geq
1-\exp\left(-n\,{\mathbb{P}}(X\geq R)\right)$ $\displaystyle\geq$
$\displaystyle 1-\exp\left(-n\,\tfrac{V}{n}\right)=1-\exp\left(-V\right).$

By using slow variation property, we obtain

$R\,{\mathbb{P}}(M_{n}\geq R)\geq
g(V/n)\left(1-\exp\left(-V\right)\right)=g(1/n)\,(1+o(1))\,\left(1-\exp\left(-V\right)\right).$

On the other hand, the integrals in (3) are uniformly bounded because

$\displaystyle\int_{0}^{\infty}{\mathbb{P}}(M_{n}\leq-r)dr\leq\int_{0}^{\infty}{\mathbb{P}}(X_{1}\leq-r)dr$
$\displaystyle\leq$ $\displaystyle\int_{0}^{\infty}{\mathbb{P}}(|X_{1}|\geq
r)dr={\mathbb{E}}|X|<\infty,$

by theorem’s assumption.

Now (3) yields

${\mathbb{E}}M_{n}\geq
g(1/n)\,(1+o(1))\left(1-\exp\left(-V\right)\right)-O(1).$

Using that $g(1/n)\nearrow\infty$ and letting $V\nearrow\infty$ we obtain

${\mathbb{E}}M_{n}\geq g(1/n)\,(1+o(1)).$

$\square$

###  2.2 Upper bound

###### Proposition 3

Under theorem’s assumptions, it is true that

${\mathbb{E}}M_{n}\leq g(1/n)\,(1+o(1)),\quad\quad as\,\,n\to\infty.$

Proof: For every $R>0$ we have

$\displaystyle{\mathbb{E}}M_{n}$ $\displaystyle\leq$
$\displaystyle\int_{0}^{\infty}{\mathbb{P}}(M_{n}\geq
r)\,dr\leq\int_{0}^{R}1\,dr+\int_{R}^{\infty}{\mathbb{P}}(M_{n}\geq r)\,dr$
(4) $\displaystyle=$ $\displaystyle
R\left(1+\int_{1}^{\infty}{\mathbb{P}}(M_{n}\geq Rh)\,dh\right).$

We will now show that under the choice $R=R_{n}:=g(1/n)$ the integral is
asymptotically negligible.

Recall the elementary inequalities

$\displaystyle 1-p$ $\displaystyle\geq$ $\displaystyle e^{-2p},\qquad
0<p<1/2,$ $\displaystyle u$ $\displaystyle\geq$ $\displaystyle 1-e^{-u},\qquad
u>0.$

By combining them, for $p\in(0,1/2)$ we have

$1-(1-p)^{n}\leq 1-e^{-2np}\leq 2np.$

Applying this inequality to $p:={\mathbb{P}}(X\geq r)$ we have, for all $r$
large enough,

${\mathbb{P}}(M_{n}\geq r)=1-(1-{\mathbb{P}}(X\geq r))^{n}\leq
2n\,{\mathbb{P}}(X\geq r).$

Therefore, for all $R$ large enough and all $n\in{\mathbb{N}}$,

$\int_{1}^{\infty}{\mathbb{P}}(M_{n}\geq Rh)\,dh\leq
2n\int_{1}^{\infty}{\mathbb{P}}(X\geq Rh)\,dh.$ (5)

Next, we use the slow variation property of the function $g$. By Karamata
representation theorem [5], for every fixed $\delta>0$ and $V>0$ we have, for
all $\rho$ small enough and all $\ell\leq 1$

$\frac{g(\ell\rho)}{g(\rho)}<(1+\delta)\ell^{-1/V}.$

Let $h\geq 1+\delta$. Then choose
$\ell:=\left(\frac{h}{1+\delta}\right)^{-V}\leq 1$ which satisfies

$(1+\delta)\,\ell^{-1/V}=(1+\delta)\,\frac{h}{1+\delta}=h.$

We obtain $g(\ell\rho)<h\,g(\rho)$. It follows that

${\mathbb{P}}\left(X\geq
g(\rho)h\right)\leq{\mathbb{P}}\left(X>g(\ell\rho)\right)\leq\ell\rho=\left(\frac{h}{1+\delta}\right)^{-V}\rho.$

Integration of this bound yields for $V>1$

$\displaystyle\int_{1+\delta}^{\infty}{\mathbb{P}}\left(X\geq
g(\rho)h\right)\,dh$ $\displaystyle\leq$
$\displaystyle\rho\int_{1+\delta}^{\infty}\left(\frac{h}{1+\delta}\right)^{-V}\,dh$
$\displaystyle=$
$\displaystyle\rho\,(1+\delta)\int_{1}^{\infty}\tilde{h}^{-V}d\tilde{h}=\frac{\rho(1+\delta)}{V-1}.$

We choose here $\rho=\rho_{n}:=1/n$ and $R=R_{n}:=g(\rho_{n})=g(1/n)$ to get

$n\int_{1+\delta}^{\infty}{\mathbb{P}}(X\geq R_{n}h)\,dh\leq
n\cdot\frac{1+\delta}{n(V-1)}=\frac{1+\delta}{V-1}.$

for all $n$ large enough. By letting $V\nearrow\infty$, we have

$n\int_{1+\delta}^{\infty}{\mathbb{P}}(X\geq R_{n}h)\,dh=o(1),\qquad\textrm{as
}n\to\infty,$

for every fixed $\delta$.

Using $R_{n}=g(1/n)$ we also have a trivial bound

$n\int_{1}^{1+\delta}{\mathbb{P}}(X\geq R_{n}h)\,dh\leq
n\,\delta\,{\mathbb{P}}(X>R_{n})\leq n\,\delta\,\frac{1}{n}=\delta.$

By summing up the bounds for two integrals, we obtain

$n\int_{1}^{\infty}{\mathbb{P}}(X\geq R_{n}h)\,dh\leq\delta+o(1),$

and letting $\delta\searrow 0$ we have

$n\int_{1}^{\infty}{\mathbb{P}}(X\geq R_{n}h)\,dh=o(1).$

By combining this estimate with (5) and (4) we finally get

${\mathbb{E}}M_{n}\leq R_{n}(1+o(1))=g(1/n)(1+o(1)).$

$\square$

##  3 Relations between $M_{n}$ and ${\mathcal{M}}_{n}$

The goal of this section is to connect the expectations of $M_{n}$ and
${\mathcal{M}}_{n}$ and to show that under theorem’s assumptions (1) implies
(2).

We start with a trivial estimate

${\mathcal{M}}_{n}\leq\sum_{i=1}^{n}\max_{1\leq j\leq
n}X_{ij}:=\sum_{i=1}^{n}M_{n}^{(i)},$

where $M_{n}^{(i)}$ are independent copies of $M_{n}$. Hence, using (1),

${\mathbb{E}}{\mathcal{M}}_{n}\leq\sum_{i=1}^{n}{\mathbb{E}}M_{n}^{(i)}=n\,{\mathbb{E}}M_{n}=n\,g(1/n)\,(1+o(1)).$
(6)

In order to get an opposite bound, we have to recall the ”greedy method” of
constructing a quasi-optimal permutation $\pi^{*}$ providing sufficiently
large value of the assignment process. In Gaussian context, this method was
introduced in [15] and later in [11].

Let $[i]:=\\{1,2,\dots,i\\}$. Define

$\pi^{*}(1):=\arg\max\limits_{j\in[n]}X_{1j},$

and let for all $i=2,\dots,n$

$\pi^{*}(i):=\arg\max\limits_{j\not\in\pi^{*}([i-1])}X_{ij}.$

It is natural to call this strategy greedy, because at every step we consider
the row $i$, take the maximum of its available elements (without considering
the influence of this choice on subsequent steps) and then forget the row $i$
and the corresponding column $\pi^{*}(i)$.

We stress that at every step we make choice from i.i.d. random variables with
distribution ${\mathcal{P}}$. They are also independent on the variables used
at other steps but we don’t need this fact here. The number of variables used
at consequent steps is decreasing from $n$ to $1$.

By using the greedy method, we have

${\mathcal{M}}_{n}\geq\sum_{i=1}^{n}X_{i\,\pi^{*}(i)}=\sum_{i=1}^{n}\max\limits_{j\not\in\pi^{*}([i-1])}X_{ij}:=\sum_{i=1}^{n}L_{i},$

where $L_{i}$ is equidistributed with $M_{n-i+1}$.

Hence, for every $\delta\in(0,1)$

$\displaystyle{\mathbb{E}}{\mathcal{M}}_{n}$ $\displaystyle\geq$
$\displaystyle\sum_{i=1}^{n}{\mathbb{E}}L_{i}=\sum_{i=1}^{n}{\mathbb{E}}M_{i}$
$\displaystyle\geq$ $\displaystyle(n-[\delta n])\ {\mathbb{E}}M_{[\delta
n]}+\delta n\,{\mathbb{E}}M_{1}$ $\displaystyle=$
$\displaystyle(1-\delta)n\,g(1/[\delta n])\,(1+o(1))+O(n)$ $\displaystyle=$
$\displaystyle(1-\delta)n\,g(1/n)\,(1+o(1)),$

where we consequently used that the sequence $({\mathbb{E}}M_{n})$ is non-
decreasing, the asymptotics (1), the slow variation property of the function
$g$, and the fact that $g(1/n)\nearrow\infty$ as $n\to\infty$.

Since $\delta$ can be chosen arbitrarily small, we obtain

${\mathbb{E}}{\mathcal{M}}_{n}\geq n\ g(1/n)\,(1+o(1)).$ (7)

Now inequalities (6) and (7) prove that (1) implies (2).

##  4 Examples

In this section we provide the main asymptotic terms of ${\mathbb{E}}M_{n}$
and ${\mathbb{E}}\mathcal{M}_{n}$ as $n\to\infty$ for some standard
distributions satisfying assumptions of Theorem 1.

Distribution $\mathcal{P}$ | Tail quantile function $g(p)$ | ${\mathbb{E}}M_{n}$ | ${\mathbb{E}}\mathcal{M}_{n}$
---|---|---|---
$\mathcal{N}\,(0,1)$ | $\sqrt{2}\,\mathrm{erfc}^{-1}(2p)$ | $\sqrt{2\log{n}}$ | $n\sqrt{2\log{n}}$
$\mathrm{Exp}\,(1)$ | $-\log{p}$ | $\log{n}$ | $n\log{n}$
$\mathrm{Gumbel}\,(0,1)$ | $-\log{\log{\frac{1}{1-p}}}$ | $\log{n}$ | $n\log{n}$
$\mathrm{Laplace}\,(0,1)$ | $-\log{(2p)},\,0<p<1/2$ | $\log{n}$ | $n\log{n}$
$\mathrm{Poisson}\,(\lambda)$ | $-\log{p}/\log{\log{(1/p)}}$ | $\log{n}/\log{\log{n}}$ | $n\log{n}/\log{\log{n}}$

## Acknowledgement

We are grateful to Professor Th. Tkotz for sending us preprint [7] and for
other useful advise.

The work supported by Russian Science Foundation Grant [21-11-00047].

## References

  * [1]
  * [2] Aldous, D. J., 2001. The $\zeta(2)$ limit in the random assignment problem. Random Structures & Algorithms, 18, No. 4, 381–418.
  * [3] Aldous, D. J., 1992\. Asymptotics in the random assignment problem. Probability Theory and Related Fields, 93, No. 4, 507–534.
  * [4] Alm, S. E. and Sorkin, G. B., 2002. Exact expectations and distributions for the random assignment problem. Combinatorics, Probability and Computing, 11, No. 3, 217–248.
  * [5] Bingham, N. H., Goldie, C. M., Teugels, J. L. Regular variation, Encycl. Math. Appl., 27, Cambridge Univ. Press, 1989.
  * [6] Buck, M. W., Chan, C. S., and Robbins, D. P., 2002. On the expected value of the minimum assignment. Random Structures & Algorithms, 21, No. 1, 33–58.
  * [7] Cheng Y., Liu Y., Tkocz T. and Xu A., 2021. Typical values of extremal-weight combinatorial structures with independent symmetric weights. Preprint.
  * [8] Coppersmith, D. and Sorkin, G. B., 1999. Constructive bounds and exact expectations for the random assignment problem. Random Structures & Algorithms, 15, No.2, 113–144.
  * [9] Dotsenko, V. S., 1993. Physics of spin-glass state. Phys. Usp. 36, No.6, 455–485.
  * [10] Dotsenko, V. S., 2000. Exact solution of the random bipartite matching model. Journal of Physics A: Mathematical and General, 33, No.10, 2015–2030.
  * [11] Lifshits, M. and Tadevosian, A., 2021. Gaussian assignment process. Preprint, https://arxiv.org/abs/2107.04883.
  * [12] Linusson, S. and Wästlund, J., 2000. A generalization of the random assignment problem. Preprint, https://arxiv.org/abs/math/0006146.
  * [13] Mézard, M. and Parisi, G., 1987. On the solution of the random link matching problem. Journal de Physique, 48 (9), 1451–1459.
  * [14] Mézard, M., Parisi, G., and Virasoro, M. A., 1987. Spin glass theory and beyond. Ser.: World Scientific Lecture Notes in Physics, Vol.9, World Scientific.
  * [15] Mordant G. and Segers J., 2021, Maxima and near-maxima of a Gaussian random assignment field. Statistics & Probability Letters. 173. 109087. 10.1016/j.spl.2021.109087.
  * [16] Nair, C., Prabhakar, B. and Sharma, M., 2005. Proofs of the Parisi and Coppersmith-Sorkin random assignment conjectures. Random Struct. Alg., 27: 413–444.
  * [17] Parisi, G., 1998. A conjecture on random bipartite matching. Preprint, https://arxiv.org/abs/cond-mat/9801176.
  * [18] Steele, J. M., 1997. Probability theory and combinatorial optimization. Ser.: CBMS-NSF Regional Conference Series in Applied Mathematics, Vol.69.
  * [19] Wästlund J., 2009. An easy proof of the $\zeta(2)$ limit in the random assignment problem. Electron. Commun. Probab. 14 261–269, 2009.


# Terahertz-Band Non-Orthogonal Multiple Access: System- and Link-Level
Considerations

Ahmed Magbool  Hadi Sarieddeen  Nour Kouzayha  Mohamed-Slim Alouini  and Tareq
Y. Al-Naffouri  The authors are with the Computer, Electrical and Mathematical
Sciences and Engineering Division (CEMSE), King Abdullah University of Science
and Technology (KAUST), Thuwal 23955-6900, Makkah Province, Saudi Arabia.
E-mail:{ahmed.magbool, hadi.sarieddeen, nour.kouzayha, slim.alouini,
tareq.alnaffouri }$@$ kaust.edu.sa

###### Abstract

Non-orthogonal multiple access (NOMA) communications promise high spectral
efficiency and massive connectivity, serving multiple users over the same
time-frequency-code resources. Higher data rates and massive connectivity are
also achieved by leveraging wider bandwidths at higher frequencies, especially
in the terahertz (THz) band. This work investigates the prospects and
challenges of combining these algorithmic and spectrum enablers in THz-band
NOMA communications. We consider power-domain NOMA coupled with successive
interference cancellation at the receiver, focusing on multiple-input
multiple-output (MIMO) systems as antenna arrays are crucial for THz
communications. On the system level, we study the scalability of THz-NOMA
beamforming, clustering, and spectrum/power allocation algorithms and motivate
stochastic geometry techniques for performance analysis and system modeling.
On the link level, we highlight the challenges in channel estimation and data
detection and the constraints on computational complexity. We further
illustrate future research directions. When properly configured and given
sufficient densification, THz-band NOMA communications can significantly
improve the performance and capacity of future wireless networks.

##  I Introduction

Following the successful deployment of millimeter-wave (mmWave) communications
in fifth-generation (5G) wireless systems, the terahertz-(THz) band is being
investigated as an enabler that supports larger bandwidths between $0.3$ and
$10\,\mathrm{THz}$ in 6G and beyond. However, in addition to the gap in THz
device technologies, several THz signal processing challenges should be first
overcome [1]. THz signal propagation results in narrow beams, high path
losses, low diffraction, high scattering, high sensitivity to blockages, and
significant differences between line-of-sight (LoS) and non-line-of-sigh
(NLoS) path gains [2]. Infrastructure enablers such as multiple-input-
multiple-output (MIMO) technology can mitigate such problems by packing many
antennas in arrays of tiny footprints at the transmitter and the receiver,
enhancing the beamforming and multiplexing gains despite the ill-conditioned
channels (see [1] and references therein).

Orthogonal multiple access (OMA) techniques have been traditionally utilized
to enhance system throughput and satisfy the required quality of service (QoS)
requirements. OMA systems serve multiple users using orthogonal time,
frequency, or code resources to prevent intra-cell interference between users.
However, such orthogonality limits the number of users served simultaneously
by the same base station (BS). Non-orthogonal multiple access (NOMA)
techniques have been proposed to enhance spectral efficiency, distinguishing
users in the power domain. In NOMA, multiple users are allowed to share the
same time-frequency-code resources, and interference is mitigated using signal
processing techniques [3, 4]. Although spectrum scarcity is not the primary
concern at the THz band, there are several motivations for THz-NOMA. Due to
hardware limitations, the THz band has not been fully conquered yet, and hence
efficient utilization of the available bandwidths remains crucial. Higher
spectral efficiency sustains higher data rates using fewer costly THz
resources. Moreover, diversity in the power domain could help NOMA users to
overcome the highly correlated THz channels of users. Furthermore, THz-NOMA
can enhance user fairness by bridging the severe power difference between the
received signals of users at different distances from the BS through proper
power allocation schemes [4].

THz-band NOMA is challenging because of the narrow THz beamwidths, high-
dimensional ill-conditioned (correlated) THz channels, and processing
complexity (mainly in detection) under Tbps constraints. The main challenges
and their prospective solutions are detailed in Table I. In this article, we
discuss such challenges, starting by describing the THz-NOMA system model in
Sec. II. We address system-level considerations in Sec. III, including
beamforming, user clustering, spectrum and power allocation, and using
stochastic geometry as a performance analysis tool. Afterward, we discuss
link-level considerations of channel estimation and data detection in Sec. IV.
We provide end-to-end simulations of THz-NOMA systems in Sec. V, demonstrating
the performance and complexity tradeoffs. We then compare NOMA with another
broadcasting scheme, namely multi-user linear precoding in Sec. VI. We
illustrate future research directions and conclude the paper in Sec. VII.

TABLE I: Main THz-band NOMA challenges and mitigation strategies (Numbers from TeraMIMO channel simulator [2]; carrier frequency of $300\,\mathrm{GHz}$). Challenge |  Example |  Mitigation strategies / opportunities
---|---|---
Narrow THz beamwidths |  For a $25\,\mathrm{dBi}$ sub-array gain, the azimuth and elevation angular spreads do not exceed $11.42^{\circ}$ [2] |  Beamwidth control - multi-beamforming
User fairness |  A user at $5\,\mathrm{m}$ distance from the transmitter experiences 14 dB more path loss than another user at 1 m (according to the free space model) |  Proper user clustering and power and spectrum allocation
High-dimensional channels |  It is typical to transmit 64 data streams over 64 SAs at both the transmitter and the receiver, with 16 antenna elements per SA to achieve a channel gain of $0.0163$ for a user at $5\,\mathrm{m}$ |  Computationally efficient system-level and link-level signal processing algorithms
Low-rank channels |  Channels in the true THz range are of rank $\sim$ 1 |  Leveraging sparsity in the spatial domain to reduce the complexity of baseband algorithms (fast channel updates; exploiting the spherical wave model)
High channel correlation across users |  Users at distances of 1 m and 5 m form the transmitter and $20^{\circ}$ elevation and angular spreads have channels with a correlation coefficient of 0.97 |  Reducing baseband complexity by eliminating redundant computations across users
Terabit per second (Tbps) detection |  At existing circuit clock speeds, $\sim$1000 bits need to be processed per clock cycle to achieve a Tbps [4] |  Computationally efficient SIC data detectors and decoders

##  II System Model

Although MIMO architectures are used in lower frequencies to enhance the
beamforming and multiplexing gains, their use in THz systems is more crucial
to compensate for the severe path loss. Both the transmitter and the receiver
in a THz system employ a large number of antennas, typically organized in an
array-of-subarray (AoSA) architecture; each subarray (SA) transmits one
information symbol, and multiplexing is enabled over different SAs [1]. Thanks
to the tiny wavelengths of THz waves, such multiple-antenna designs can be
realized in small device footprints.

The first step of THz-NOMA is clustering users into groups that share
communication resources, which is very challenging given the narrow width of
the THz beams. Each user cluster is typically served by a single THz beam in
beam-division multiple access (BDMA), as illustrated in Fig. 1 [3]. However,
only a small number of users can be served within a single NOMA group because
user interference severely degrades the quality of decoded signals and causes
a non-negligible delay in successive interference cancellation (SIC).
Alternatively, users per beam can be further divided into sub-clusters
(groups), where orthogonal group resources avoid inter-group interference
(i.e., OMA between groups and NOMA within a group).

The second step of NOMA is power allocation. For simplicity of illustration,
we consider two-user NOMA groups served by a BS over the same time-frequency-
code resources. Superposition coding combines the users’ signals before
transmission, assigning a higher power level to the user with worse channel
conditions. Because THz channels are LoS-dominant and subject to a severe path
loss, the distance between users and the BS dictates to a great extent the
channel condition.

The third step of NOMA is successive interference cancellation (SIC) at the
receiver side. Since the weak user (of worse channel conditions) is assigned
higher power, its component in the received signal is more significant; it can
directly decode its signal, treating the other user’s signal as noise. The
strong user, on the other hand, first decodes the weak user’s signal, treating
its own signal as noise, and then applies SIC to remove the interfering signal
and decode its own signal.

The benefits of THz-NOMA can only be fully realized following optimizations
and efficient signal processing at both the system and link levels. We address
the corresponding system- and link-level prospects and challenges in Sections
III and IV, respectively.

Figure 1: NOMA systems with uniform beamforming, beamwidth control, and multi-
beamforming.

##  III System-Level Considerations

At the system level, we consider multi-user beamforming, NOMA user clustering,
spectrum and power allocation, and motivate the use of stochastic geometry as
a performance analysis technique.

###  III-A Beamforming

THz-NOMA beamforming aims at forming multiple beams of proper direction and
width to serve multiple users. The half-power elevation and azimuth beam
angles are inversely proportional to the SA beamforming gain [2], as
illustrated in Fig. 2. For example, a beam with elevation and azimuth angular
spreads of $27^{\circ}$ can have a SA gain of $17.3\,\mathrm{dBi}$, which is
usually not sufficient for medium-distance communications [2]. The probability
of serving $k$ uniformly distributed users (polar uniform distribution) in an
angular spread of $27^{\circ}$ is small ($0.02$ for $k\\!=\\!2$; so an average
of $100$ users are required to serve two users per beam). Such observations
raise concerns about the practicality of NOMA at high-frequency bands.

Figure 2: Heat map of the subarray gain in dBi as a function of the azimuth
and elevation angular spreads.

Beamwidth control techniques [5] can increase the probability of serving
multiple users by widening the main lobe (e.g., beam $7$ in Fig. 1), for
example. THz-band ultra-massive MIMO (UM-MIMO) beamwidth control should
achieve a trade-off between complexity and performance. Simple solutions such
as conventional beamforming (CBF), where the number of active antennas is
reduced to widen the analog beamwidth with constant modulus phased sifters
(PSs), are favorable due to their low complexity. However, such schemes lack
the flexibility to change beamwidths with user activities and densities. More
computationally complex optimization-based beamwidth control methods can
maximize the worst-user system sum rate (max-min fairness (MMF)). Furthermore,
multi-beamforming methods can steer a beam generated from one radio-frequency
(RF) chain to serve angularly distant users (e.g., beam $6$ in Fig. 1) [6].
However, multi-beamforming reduces the SA gain as the transmit power will be
split into different physical directions.

There are a number of THz-NOMA beamforming solutions in the literature. Simple
approaches such as uniform beamforming minimize system energy consumption but
are not efficient solutions at low user densities as energy can be wasted
forming beams in no-user directions. A more advanced solution in [7] first
clusters users and then forms beams with appropriate directions and widths
that suit the NOMA groups. In particular, the beamforming problem is jointly
optimized with power and bandwidth allocation following user clustering to
maximize system throughput.

###  III-B User Clustering

Clustering is the task of distributing users into small groups in which NOMA
is applied. Although BDMA mitigates inter-beam interference, the interference
between different groups within the same beam should still be avoided, perhaps
using OMA techniques. A basic user-pairing approach is based on the
corresponding user channel conditions, where weak users are typically paired
with strong ones. If the LoS is blocked, the corresponding user is directly
classified as weak as the difference between THz LoS and NLoS components is
significant. Otherwise, users are classified based on their distances from the
BS. So users inside a circle of radius $r$ (smaller than cell radius) are
strong users and the remaining ones are considered weak. Distance-based user
distinction is particularly plausible at THz frequencies due to severe path
losses.

Similar-distance grouping (SDG) [7] is another candidate scheme for THz-NOMA
clustering. In SDG, users are clustered according to their distances from the
BS; the nearest two users form the first group, the second nearest the second
group, etc. As opposed to strong-weak user pairing, SDG pairs a strong user
with another strong user and vice versa. Consequently, spectrum allocation can
be facilitated, where weak-user groups get assigned sub-bands suffering from
less path loss (more details in the next subsection).

Clustering can be updated rather than performed from scratch every time by
utilizing proper clustering update algorithms. By doing so, an optimization-
based clustering technique can be solved to determine the initial association
of the users and then a clustering update technique can be invoked based on
the users’ mobility without the need to solve the optimization problem again.
Unsupervised and supervised machine-learning and reinforcement learning tools
can also be used to cluster users based on features such as geographical
location and channel condition to indicate the best user clustering policy.

###  III-C Spectrum Allocation

Spectrum allocation is the problem of assigning carrier frequencies to NOMA
groups. Frequency-based OMA techniques such as orthogonal frequency-division
multiple access (OFDMA) can still be used to distinguish NOMA groups within a
beam. However, due to THz-band high frequency-selectivity, hardware
constraints and molecular absorption, THz path losses are more severe at
particular sub-bands that need to be avoided. The available bandwidth is thus
divided into smaller transmission windows, the centers of which enjoy less
path loss than the sides [8]. For instance, a total bandwidth of
$45\,\mathrm{GHz}$ is available for transmission at a central frequency of
$300\,\mathrm{GHz}$ (transmission window from $227.5$ to
$322.5\,\mathrm{GHz}$) [8]. Efficient THz carrier frequency assignment in OMA
can leverage this property to allocate window centers to faraway users and
window sides to closer users in a distance-aware multi-carrier (DAMC) approach
[8]. The same concept can be applied to NOMA by allocating each NOMA group a
sub-band at the transmission window center or side, depending on the distance
from the BS. To this end, a group distance metric needs to be defined, the
center of mass of user distances in a group, for example. The distinction
between the group distances is more clear if users are clustered using the SDG
scheme described in Sec. III-B.

###  III-D Power Allocation

Power allocation is the process of determining the power to be assigned for
each NOMA user, with the objective of maximizing sum rate, energy efficiency,
or system fairness. In an AoSA architecture, power allocation can be expressed
in terms of power assigned per transmitting SA, accounting for antenna and
beamforming gains (per SA). Conventional power allocation techniques are
constant or optimization-based. Intuitively, users with worse channel
conditions should be assigned higher power levels for fairness, and low
complexity assignments (without the need for solving complex optimizations)
are favored. For example, users can be ordered based on their channel
conditions; the user with the worst condition gets assigned a certain power,
the second weak user a fraction of that power; and continue this way while
observing the maximum transmit power constraint [9]. In optimization-based
solutions, power allocation can be coupled with beamforming and user pairing.
However, solving such global optimizations is computationally complex,
especially using iterative algorithms. Besides, such optimizations are usually
non-convex, and a global solution is not guaranteed.

As an illustration, consider a NOMA group with three users at distances $1$,
$1.2$, and $7\,\mathrm{m}$ from the BS; assume a transmit power budget of
$20\,\mathrm{mW}$ and a power fraction of $0.5$. The power allocated to the
furthest, middle, and nearest user is $11.43$, $5.71$, and $2.86\,\mathrm{W}$,
respectively, which results in a receive power of $-88$, $-76$ and
$-77.4\,\mathrm{dBm}$ according to the free space model. This approach clearly
lacks fairness, which can be enhanced by linking power allocation to user
channel gains, as in fractional transmit power allocation (FTPA) [9], for
example. Applying FTPA to our example, with the assumption that the channel
gain is inversely proportional to the square of the distance from a BS, the
power distribution becomes $19$, $0.6$, and $0.4\,\mathrm{mW}$ for the
furthest, middle, and nearest user, respectively. This will result in a
received power of around $-86\,\mathrm{dBm}$ for all users.

###  III-E Stochastic Geometry for THz-NOMA

Stochastic geometry (SG) can be used as a powerful model-driven tool that
provides a mathematical framework to capture the properties of THz networks
and quantify the benefits of THz-NOMA. Several research works have leveraged
SG for the performance analysis of power-domain NOMA networks operating on
sub-$6$ GHz and mmWave frequency bands [10]. However, there is still a gap in
the literature in extending the available models to account for specific
characteristics of THz-NOMA. In the SG analysis of multi-cell NOMA systems,
clustering algorithms that capture the users to be simultaneously served in
the same resource block should be clearly defined. This is important because
the distances between the users and the serving BS within a NOMA group are
dependent on the adopted clustering technique [11]. A more general setup is
also considered by modeling the users as a Poisson cluster process (PCP). The
beamforming technique and the narrow beamwidth will affect user clustering and
hence the distance distributions of the serving and the interfering BSs. The
accuracy of the SIC process at the receiver should also be taken into
consideration, which can be modeled by introducing a fraction between $0$ and
$1$ in the SINR expression (for more details, see [10] and the references
therein).

In addition to previous modeling choices, several modifications should be
applied to conventional SG frameworks to account for THz-specific
characteristics in THz-NOMA networks. Typically: 1) THz antennas should be
modeled as directional antennas. 2) 3D modeling of THz networks is mandatory
as it accounts for the effect of vertical heights. While it may be acceptable
to ignore the heights in sub-$6$ GHz networks as they enjoy large
communication distances, such an assumption is inaccurate for dense THz
networks with limited communication ranges. This will require the use of 3D
directional antenna models such as the 3D pyramidal-plus-sphere-sectored
antenna model [10, 12] at THz systems. 3) The impact of blockages can be
captured by a LoS/NLoS probabilistic model that identifies each link’s
condition based on the link length and the propagation environment. In THz
systems, blocking LoS paths is mainly due to attributes such as dynamic self-
blockage and human blockage and static blockages like buildings in outdoor
environments and walls or furniture in indoor environments [10, 12] 4) Since
THz beams are usually narrow and as THz communication is more sensitive to
blockages compared to sub-$6$ GHz, the small scale fading can be ignored in
the analysis of THz networks. 5) As opposed to lower frequencies, THz
communication is affected by water molecules in the environment. To account
for this, the THz propagation model is a frequency-dependent exponential
function that highlights the effect of both molecular absorption loss and
spreading loss [10, 12]. 6) Finally, SG needs to account for thermal noise as
well as molecular absorption noise caused by the re-radiation of the absorbed
signal energy.

##  IV Link-Level Considerations

At the link level, we need to consider the problems of efficient channel
estimation and data detection.

###  IV-A Channel Estimation

Channel state information (CSI) is required at the transmitter for
beamforming, user pairing, and power allocation, and at the receiver for SIC
and symbol decoding. Imperfect channel estimation causes power allocation
ambiguity which reduces the NOMA sum rate. It is essential to investigate the
effect of channel mismatch caused by THz propagation characteristics such as
beam squint and spherical wave propagation. Due to the large number of
antennas at both the transmitter and receiver, computationally efficient
channel estimation solutions are crucial [2].

THz signal propagation models should also account for some widely ignored
characteristics at lower frequency bands. For instance, spherical wave
propagation (SWP) models should replace plane wave approximations that are
inaccurate at THz frequencies and may lead to channel estimation errors in
massive MIMO systems. In particular, hybrid spherical- and planar-wave channel
modeling solutions are found to be sufficiently accurate and less complex than
fully spherical wave models [13]. Similarly, the beam squint effect, caused by
frequency-independent steering vectors in wide-band massive MIMO THz
beamforming, should be accounted for [14]; ignoring such frequency dependency
results in signals being steered into unintended physical directions.

When accounted for, beam squint and SWP can result in enhanced THz-NOMA system
performance. For instance, the beam squint effect can be exploited to serve
more users as beams are split in different directions (a concept similar to
multi-beamforming). Moreover, the SWP model can better distinguish the
channels of different users as it reduces the correlation between these
channels (i.e., SWP introduces more degrees of freedom). Such characteristics
are useful when modeled carefully and integrated into the system model.

In addition, despite being problematic for NOMA, the high correlation between
user channels can be leveraged to reduce channel estimation complexity
significantly. For example, instead of estimating the CSI for every user
independently, minimal channel estimation updates on a specific CSI instance
based on geometry information can construct the channels of the neighboring
users. The computational complexity of the channel estimation process can be
reduced using multiple methods. In general, the low-rank LoS-dominated THz
channels can convert the channel estimation problem into a sparse recovery
problem in the angle domain. Since THz channels have usually low ranks, the
single-input-single-output (SISO) channel can be estimated between one
transmitting and one receiving SA and appropriate tuning is performed to find
the other entries of the channel matrix. Furthermore, channel update
algorithms can be utilized to expand the channel estimated in one spatial,
temporal or spectral domain to others.

###  IV-B Data Detection

While channel coding and decoding account for most of the computational cost
in the THz baseband, a proper design of data detectors can significantly
reduce the complexity of THz UM-MIMO systems. When deciding on the detector of
choice, it is important to note the constraints on the processing capabilities
towards achieving the promised Tbps data rates. With a few GHz clock speeds,
baseband processing should be executed at a rate of at least 1000 bits per
clock cycle to achieve Tbps. Existing THz transceivers are not yet capable of
supporting such rates. On top of detection constraints in conventional THz
systems, SIC in NOMA at least doubles (in two-user NOMA) the resources
required to decode signals; strong users have to detect both their symbols and
those of the weak users, all while satisfying the Tbps constraint of THz
communications. Due to the corresponding SIC delay, having more than two users
in a THz-NOMA group is impractical.

The maximum-likelihood (ML) detector is optimal but is incapable of meeting
the Tbps constraint because of its computational complexity (exponential
increase with MIMO dimension). Furthermore, unlike conventional low-frequency
massive MIMO systems, the high-dimensional (doubly-massive) ill-conditioned
THz channels make linear detectors such as zero-forcing inefficient (suffering
from performance loss caused by inverting low-rank channel matrices). THz MIMO
systems thus require more sophisticated, yet low-complexity non-linear
detectors such as nulling-and-cancellation (NC), chase detectors, and the
layered orthogonal lattice detector (LORD), all of which are realizations of
subspace detectors [4]. These detectors can be employed to puncture the
channel and force certain structures to help in the data detection. In
particular, subspace decompositions can provide sufficient parallelizability
in the detector and SIC NOMA stages to achieve the Tbps complexity and memory
requirements [4]. In near-static THz scenarios where the channel is flat,
preprocessing matrix decompositions can be retained across multiple frames,
further reducing complexity. Moreover, even under slight channel variations,
matrix decompositions can be updated without recomputing them from scratch.

##  V Simulations and Discussion

In this section, we provide end-to-end simulations of THz-NOMA systems
addressing system-level and link-level considerations.

###  V-A System-Level Simulation

We consider $n$ users randomly and uniformly distributed inside a 2-D circular
cell of radius $10\,\mathrm{m}$. The SA beamforming gain is set to
$16.61\mathrm{dBi}$, which results in beams with half-power azimuth and
elevation angles of $30^{\circ}$. The transmission window is centered at
$300\,\mathrm{GHz}$ with a total bandwidth of $10\,\mathrm{GHz}$. We adopt
uniform beamforming, so $12$ beams are needed to cover the cell.

At the system level, we compare three THz-NOMA baselines: (1) distance-based
clustering and power allocation [4], (2) exhaustive-search-based user
clustering and power allocation that maximizes user fairness per beam, and (3)
THz-OMA with frequency-division-multiple-access (FDMA). Each NOMA group is
assigned $1$ GHz; each OMA user is allocated $1$ GHz. The channels are
generated using the TeraMIMO channel simulator [2]; four transmitting and
receiving SAs; $64$ antennas per SA; carrier frequency of $300\,\mathrm{GHz}$.
Although the channel dimension is $4\times 4$, this is still considered an UM-
MIMO system as the transmitter and the receiver consist of $256$ antenna
elements each. We assume $20$ or 50 users per cell. Since $12$ beams are
required to cover the cell, on average, $24$ users should exist in a cell for
allocating two users per NOMA beam.

Figure 3: Fairness factor for NOMA and OMA systems with 20 and 50 users.

Given that spectral efficiency is not the main bottleneck of THz systems due
to the huge bandwidth availability, we examine user fairness with NOMA. We
introduce a fairness metric, computed as the ratio of the lowest and highest
rates in the beam. The fairness metric is 1 for equal maximum and minimum user
rates (maximum fairness) and 0 for a lowest user rate of zero; a fairness
metric closer to $1$ indicates even distribution of user rates. A weighted
average by the number of users per beam gives the overall fairness factor of
the system.

Fig. 3 depicts the fairness versus the signal-to-noise ratio (SNR) defined as
the ratio of total transmit power by the noise power. The fairness factor
increases with the number of users, which is intuitive and further highlights
the importance of user densification in THz systems. Furthermore, distance-
based NOMA outperforms OMA, and exhaustive-search-based NOMA forms an upper
bound on NOMA system fairness. The results can be partitioned into three
regions based on SNR. The first region (lower SNR regime) is noise-limited and
does not achieve good fairness. The second region introduces fairness
improvements to all schemes. However, in the third region (high SNR regime),
OMA fairness starts to decay. Therefore, despite the slight increase in the
fairness factor of NOMA systems at high SNR, NOMA algorithms can still
maintain user fairness.

###  V-B Link-Level Simulation

For link-level simulations, we consider a cell of $10\,\mathrm{m}$ radius and
100 users. We compare two detectors, nulling and cancelling (NC) and LORD,
assuming distance-based clustering and power allocation (see [4] for more
details). We study the corresponding bit-error-rate (BER) performance in Fig.
4 assuming"LoS + multi-path" THz channels, and accounting for CSI errors
caused by ignoring beam squint and SWP. LORD generally outperforms NC at a
marginal complexity cost. It can be noted that the weak user has a
significantly higher BER compared to the strong user. This can be attributed
to two main reasons. First, as opposed to the strong user, the weak user does
not apply any form of interference cancellation. Second, the weak user has a
much lower channel gain than the strong one due to the severe path loss. The
error of detecting the bits for the weak user can be propagated to degrade the
performance of the BER of the strong user as well. Furthermore, the system
performance is degraded by ignoring beam squint, more significantly and more
by ignoring SWP. The BER can increase by more than an order of magnitude when
both beam squint and SWP are neglected compared to the case of perfect CSI.
This observation further highlights the importance of designing joint THz-NOMA
clustering, power allocation, channel estimation, and detection algorithms.

Figure 4: BER versus SNR when the detectors have full knowledge of CSI (red),
do not account for beam squint (green), do not account for SWP (blue) and do
not account for both beam squint and SWP (black) for LORD detector (solid) and
NC detector (dashed). All curves are plotted for the strong user, whereas, for
clarity, only the case of perfect CSI is illustrated for the weak user (doted-
dashed for LORD detector and dotted for NC detector).

##  VI NOMA versus MU-LP

The achievable gains of MIMO NOMA at lower frequencies are disputed, citing
equivalent or superior performance provided by the much simpler multi-user
linear-precoding (MU-LP) scheme [15]. MU-LP is a special class of broadcasting
techniques that, unlike NOMA, do not employ SIC at any receiver. That is, a
superposition signal is transmitted, and each user recovers its own data by
treating other users signals’ as noise. The multiplexing gain, which measures
how fast the rate changes with SNR (system performance at high SNR), is
studied for NOMA and MU-LP in [15]. NOMA achieves lower multiplexing gains
than MU-LP under full-rank channels due to the SIC rate constraint in the
latter. In particular, the NOMA transmission rate is restricted to support
both paired users while preventing SIC error, which does not apply to MU-LP.
This argument assumes a zero-forcing beamformer (ZFBF) to mitigate intra-
cluster interference.

To check the validity of this conclusion in THz scenarios, the system capacity
versus SNR is simulated in Fig. 5 for three channels types: Gaussian (zero
mean and unit variance; typical for lower frequencies), correlated, and
orthogonal THz channels ($4\\!\times\\!4$ SA-channel at $300\,\mathrm{GHz}$
using TeraMIMO [2]). THz channels can be made orthogonal by careful spatial
tuning of SA positions as described in [4]. Due to high path loss, inter-cell
interference has a negligible effect in the THz band compared to the lower
frequencies and can be safely ignored. For MU-LP, we apply ZFBF to mitigate
intra-cluster interference in both users. ZFBF is accompanied with rate-
maximization to obtain the strong-user power matrix in NOMA (no interference
following SIC). Fig. 5 verifies that MU-LP achieves a higher multiplexing gain
(slope of sum-rate) in Gaussian and orthogonal THz channels. However, this
observation does not apply to naturally occurring THz correlated channels.
Applying ZFBF on ill-conditioned channels to cancel the interfering-user
signal corrupts the main signal of the intended user, reducing the achievable
rate. Therefore, ZFBF-based MU-LP is not a valid solution at the THz-band,
despite its favorable low complexity.

Figure 5: Sum-rate per Hz versus SNR using ZFBF-based power allocation for
Gaussian and THz channels (correlated and orthogonal) for NOMA and MU-LP
systems two-user systems.

##  VII Conclusions

In this work, we discuss the prospects and challenges of applying NOMA
techniques to THz communications. We investigate system- and link-level
considerations of THz-NOMA, highlighting the effects of massive antenna
arrays, high channel correlation, narrow beams, etc. We demonstrate that THz-
NOMA can indeed enhance system throughput and user fairness in cells with
sufficient densifications, which is a promising research direction towards 6G.
Therefore, this work provides solid grounds for future THz-NOMA research that
considers realistic THz-specific constraints.

Beside the future THz-NOMA research directions mentioned previously in the
paper, several others can be considered to achieve ultimate network promises
of maximizing area traffic capacity, spectral efficiency, user fairness, and
bandwidth utility. For example, on the algorithmic level, learning-based
solutions can address system-level THz constraints in beamforming and beam
steering, forming beams by learning user activity and dynamicity. Furthermore,
on the infrastructure level, intelligent reflecting surfaces (IRSs) can be
utilized to control the propagation environment and enhance THz channel
conditions. IRSs can boost the THz multiplexing gain and reduce the huge gap
between the LoS and NLoS components.

## Acknowledgment

Research reported in this publication was supported by King Abdullah
University of Science and Technology (KAUST).

## References

  * [1] H. Sarieddeen, M.-S. Alouini, and T. Y. Al-Naffouri, “An overview of signal processing techniques for terahertz communications,” _Proceedings of the IEEE_ , vol. 109, no. 10, pp. 1628–1665, 2021.
  * [2] S. Tarboush, H. Sarieddeen, H. Chen, M. H. Loukil, H. Jemaa, M. S. Alouini, and T. Y. Al-Naffouri, “TeraMIMO: A channel simulator for wideband ultra-massive MIMO terahertz communications,” _arXiv preprint arXiv:2104.11054_ , 2021.
  * [3] B. Makki, K. Chitti, A. Behravan, and M.-S. Alouini, “A survey of NOMA: Current status and open research challenges,” _IEEE Open Journal of the Communications Society_ , vol. 1, pp. 179–189, 2020.
  * [4] H. Sarieddeen, A. Abdallah, M. M. Mansour, M.-S. Alouini, and T. Y. Al-Naffouri, “Terahertz-band MIMO-NOMA: Adaptive superposition coding and subspace detection,” _arXiv_ , 2021.
  * [5] Z. Wei, D. W. K. Ng, and J. Yuan, “NOMA for hybrid mmWave communication systems with beamwidth control,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 13, no. 3, pp. 567–583, 2019.
  * [6] Z. Wei, L. Zhao, J. Guo, D. W. K. Ng, and J. Yuan, “Multi-beam NOMA for hybrid mmWave systems,” _IEEE Transactions on Communications_ , vol. 67, no. 2, pp. 1705–1719, 2019.
  * [7] X. Zhang, C. Han, and X. Wang, “Joint beamforming-power-bandwidth allocation in terahertz NOMA networks,” in _2019 16th Annual IEEE International Conference on Sensing, Communication, and Networking (SECON)_ , 2019, pp. 1–9.
  * [8] H. Yuan, N. Yang, K. Yang, C. Han, and J. An, “Enabling massive connections using hybrid beamforming in terahertz micro-scale networks,” in _2020 IEEE Wireless Communications and Networking Conference (WCNC)_ , 2020, pp. 1–7.
  * [9] H. Zhang, D.-K. Zhang, W.-X. Meng, and C. Li, “User pairing algorithm with SIC in non-orthogonal multiple access system,” in _2016 IEEE International Conference on Communications (ICC)_ , 2016, pp. 1–6.
  * [10] Y. Hmamouche, M. Benjillali, S. Saoudi, H. Yanikomeroglu, and M. D. Renzo, “New trends in stochastic geometry for wireless networks: A tutorial and survey,” _Proceedings of the IEEE_ , vol. 109, no. 7, pp. 1200–1252, 2021\.
  * [11] W. Yi, Y. Liu, A. Nallanathan, and M. Elkashlan, “Clustered millimeter-wave networks with non-orthogonal multiple access,” _IEEE Transactions on Communications_ , vol. 67, no. 6, pp. 4350–4364, 2019.
  * [12] K. M. S. Huq, J. Rodriguez, and I. E. Otung, “3D network modeling for THz-enabled ultra-fast dense networks: A 6G perspective,” _IEEE Communications Standards Magazine_ , vol. 5, no. 2, pp. 84–90, 2021.
  * [13] Y. Chen, L. Yan, , and C. Han, “Hybrid spherical- and planar-wave channel modeling and DCNN-powered estimation for terahertz ultra-massive MIMO systems,” 2021.
  * [14] B. Wang, M. Jian, F. Gao, G. Y. Li, and H. Lin, “Beam squint and channel estimation for wideband mmwave massive MIMO-OFDM systems,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 23, pp. 5893–5908, 2019.
  * [15] B. Clerckx, Y. Mao, R. Schober, E. Jorswieck, D. J. Love, L. H. Jinhong Yuan, G. Y. Li, E. G. Larsson, and G. Caire, “Is NOMA efficient in multi-antenna networks? A critical look at next generation multiple access techniques,” _arXiv_ , 2021.