# Contribution of hadronic light-by-light scattering to the hyperfine
structure of muonium

V. I. Korobov  BLTP JINR, Dubna, Russia  Samara University, Samara, Russia  A.
V. Eskin  Samara University, Samara, Russia  A. P. Martynenko  Samara
University, Samara, Russia  F. A. Martynenko  Samara University, Samara,
Russia

###### Abstract

The contribution of hadronic scattering of light-by-light to the hyperfine
structure of muonium is calculated using experimental data on the transition
form factors of two photons into a hadron. The amplitudes of interaction
between a muon and an electron with horizontal and vertical exchange are
constructed. The contributions due to the exchange of pseudoscalar, axial
vector, scalar and tensor mesons are taken into account.

Muonium hyperfine splitting, one meson exchange interaction, quantum
electrodynamics

###### pacs:

36.10.Dr, 12.20.Ds, 14.40.Aq, 12.40.Vv

##  I Introduction

Exotic atoms such as the muonium atom, positronium atom, positronium ion,
muonic hydrogen, etc. play a very important role in modern physics. Precise
study of their energy levels, decay widths is such a direction of fundamental
research, within which one can look for manifestations of new interactions of
particles. Although such systems do not exist for a long time by the standards
of conventional systems, their creation and experimental study allows one to
look into a field of research that is inaccessible when working with stable
atoms and molecules. We can say that the study of exotic systems, along with
collider physics, is a tool for understanding reality beyond the Standard
Model.

Electromagnetic two-particle bound states make it possible to test one of the
most successful theories of particle interaction - quantum electrodynamics.
Theoretical calculations of the energy levels of the simplest bound states in
quantum electrodynamics have reached a very high accuracy [1, 2, 3, 4]. But
since the accuracy of the experimental study of energy levels has steadily
increased in recent decades, this has led to the need to study not only the
electromagnetic high order contributions but also the contributions of weak
and strong interactions in such systems. For example, the contribution of
hadronic vacuum polarization has already reached the level of experimental
verification for the anomalous magnetic moment (AMM) of the muon, hyperfine
splitting in muonium, in the Lamb shift, and the hyperfine structure of muonic
hydrogen. The most acute situation with the calculation of hadronic
contributions has developed for the AMM muon [5, 6, 7]. But for the other two
problems, the hadronic contributions also become significant, taking into
account the increasing precision of the experiment.

The study of the fine and hyperfine structure (HFS) of muonium has been
central to the study of quantum electrodynamics for decades, since in this
purely lepton system of different leptons there are no nuclear structure
effects, which have always been the main theoretical uncertainty [1, 2, 3]. In
recent years, new more accurate experimental studies related to muonium have
already begun. The Mu-MASS (MuoniuM lAser SpectroScopy) collaboration aims to
measure the $1S-2S$ transition in muonium with a final uncertainty of 10 kHz,
providing a 1000-fold improvement on accuracy [8]. New result of measurement
of the n=2 Lamb shift in muonium comprises an order of magnitude improvement
upon the previous best measurement [9]. The MuSEUM (Muonium Spectroscopy
Experiment Using Microwave) collaboration performed a new precision
measurement of the muonium ground-state hyperfine structure at J-PARC using a
high-intensity pulsed muon beam [10]. The accuracy of the experimental result
in [10] is 4 kHz and is still less than the accuracy of the previous
experiment in 1999 [11]. One can consider experiments with muonium for more
precise determination on the mass ratio $m_{\mu}/m_{e}$, for the test of the
Standard Model with greater accuracy and possibly, for revealing the source of
previously unaccounted interaction between particles forming the bound state
in QED. According to the work [12], the theory predicts
$\nu_{HFS}=4463302872(515)~{}Hz$, $\delta=1.3\times 10^{-7}$, where the most
part of the uncertainty (511 Hz) is dominated by the measurement of the ratio
$m_{\mu}/m_{e}$ (120 ppb). Therefore, from a comparison of the theoretical and
new experimental results for muonium HFS, one can obtain a more accurate value
for the mass ratio $m_{\mu}/m_{e}$. The MuSEUM collaboration aims to precisely
measure the ground-state hyperfine splitting of muonium atoms with the
accuracy 1 ppb [13].

Such a high experimental accuracy of measuring the hyperfine structure of
muonium at a level of 1 Hz requires corresponding theoretical calculations of
various high-order corrections to the fine structure constant. Such
calculations have been carried out over the years by various groups. In this
paper, we study only one of the contributions to the hyperfine structure
connected with the effect of light-by-light scattering, which leads to the
production of various mesons in the intermediate state. In the quark model,
such processes are determined by the production of a pair of light quarks and
antiquarks in the $\gamma^{\ast}\gamma^{\ast}$ interaction, which can then
form a light meson. The corresponding interaction amplitudes are shown in Fig.
1. They can be divided into two parts, which we call vertical and horizontal
exchanges. In our previous work [14], we investigated the contribution
connected with the horizontal exchange of pseudoscalar mesons. A more complete
study of these processes was carried out in Ref. [15], in which, along with
horizontal exchanges, the contribution of vertical exchanges, including axial
vector mesons, was also investigated. In our recent papers, we calculated the
hadronic contributions of light-by-light scattering into the fine and
hyperfine structure of muonic hydrogen [16, 17, 18, 19, 20] (see also Refs.
[21, 22, 23, 24]) and showed that such processes must be taken into account
when obtaining the total value of a specific energy interval, taking into
account the ever-increasing accuracy experiments carried out by the CREMA
collaboration [25, 26, 27, 28], as well as other collaborations are planned
[29, 30, 31]. The purpose of this work is to calculate all possible meson
contributions (pseudoscalar, scalar, axial vector and tensor) to the hyperfine
splitting (HFS) in muonium and to estimate the possible total contribution
from such interactions. The factor determining the order of the contribution,
$m_{e}^{3}\alpha^{7}/\Lambda^{2}h\sim 0.04$ Hz where $\Lambda$ is typical
hadronic mass near 1 GeV, is estimated to be not very large due to recoil
effects and the nature of the hadronic interaction itself. Nevertheless, the
study of such contributions in the hyperfine structure is of interest in
connection with an increase in the accuracy of measurements. Thus, for
example, in the case of muonic hydrogen, hadronic effects of light-by-light
scattering turn out to be rather significant both in the Lamb shift and in the
hyperfine splitting [16, 17, 18, 19, 20].

##  II Contribution of axial vector mesons

We begin the discussion of the contributions of axial vector mesons from the
vertical exchange amplitudes in Fig. 1(c). The diagram has a vertex of the
transition of two virtual photons to an axial vector meson, for which the
following parametrization is used [32, 16]:

$T^{\mu\nu}(k_{1},k_{2})=4\pi
i\alpha\varepsilon_{\mu\nu\alpha\beta}(k_{1}^{\alpha}k_{2}^{2}-k_{2}^{\alpha}k_{1}^{2})\varepsilon_{A}^{\beta}A(t^{2},k_{1}^{2},k_{2}^{2}),$
(1)

where $k_{1}=k$, $k_{2}=t-k$ are four-momenta of virtual photons,
$t=p_{1}-q_{1}=(0,{\bf t})$ is the four-momentum of the meson, $p_{1}$,
$p_{2}$ are four-momenta of electron and muon in initial state, $q_{1}$,
$q_{2}$ are four-momenta of electron and muon in final state, $M_{A}$ is the
mass of axial vector meson. Note that the axial vector decay into two real
photons is forbidden by Landau-Yang theorem but the process with one virtual
photon can already take place. To pick out the electron-muon states with a
certain spin, we use projection operators constructed from the wave functions
of the particles in their rest frame:

$\hat{\Pi}_{S=0}=[u(0)\bar{v}(0)]_{S=0}=\frac{(1+\gamma^{0})}{2\sqrt{2}}\gamma_{5},~{}~{}~{}\hat{\Pi}_{S=1}=[u(0)\bar{v}(0)]_{S=1}=\frac{(1+\gamma^{0})}{2\sqrt{2}}\hat{\varepsilon}.$
(2)

Figure 1: Hadronic light-by-light scattering amplitudes with horizontal and
vertical exchanges. Wavy line corresponds to the virtual photon. The bold dot
denotes the form factor of the transition of two photons into a meson.

As a result, the general expression for the interaction amplitude in Fig.1(c)
can be transformed to the following trace:

$i{\cal
M}^{c}=\frac{\alpha^{2}(Z\alpha)^{2}}{16m_{1}^{2}m_{2}^{2}}\int\frac{d^{4}k}{\pi^{2}}\frac{A(t^{2},k_{1}^{2},k_{2}^{2})}{(k^{2})^{2}}\int\frac{d^{4}r}{\pi^{2}}\frac{A(t^{2},r_{1}^{2},r_{2}^{2})}{(r^{2})^{2}}\frac{\varepsilon_{\mu\nu\alpha\beta}(k_{1}^{\alpha}k_{2}^{2}-k_{2}^{\alpha}k_{1}^{2})}{(k^{2}-2k_{0}m_{1})}\times$
(3)
$\frac{\varepsilon_{\sigma\lambda\rho\omega}(r_{1}^{\rho}r_{2}^{2}-r_{2}^{\rho}r_{1}^{2})}{(r^{2}-2r_{0}m_{2})}D^{\beta\omega}(t)Tr\Bigl{[}(\hat{q}_{1}+m_{1})\gamma^{\nu}(\hat{p}_{1}-\hat{k}+m_{1})\gamma^{\mu})\hat{p}_{1}+m_{1})\hat{\Pi}_{S=1,0}(\hat{p}_{2}-m_{2})\times$
$\gamma^{\sigma}(\hat{r}_{1}-p_{2}+m_{2})\gamma^{\lambda}(\hat{q}_{2}-m_{2})\hat{\Pi}_{S=1,0}\Bigr{]}D^{\beta\omega}(t),$

where $m_{1}$, $m_{2}$ are the masses of electron and muon correspondingly,
$k_{1}=k$, $k_{2}=t-k$ are four-momenta of virtual photons in one loop,
$r_{1}=k$, $r_{2}=t-r$ are four-momenta of virtual photons in other loop.
$D^{\beta\omega}(t)$ is the propagator of axial-vector meson. After taking the
trace in leading order in $\alpha$ and a number of simplifications, the
amplitude numerator for hyperfine splitting can be represented as:

$N^{c}_{AV}=\frac{1}{3}k^{2}r^{2}{\bf k}^{2}{\bf r}^{2},$ (4)

where the index (c) denotes the contribution of the amplitude in Fig. 1(c).
For the purpose of further integration over loop momenta, we pass to the
Euclidean space:

$k^{2}\to-k^{2},~{}~{}~{}r^{2}\to-r^{2},~{}~{}~{}k_{0}^{2}\to-
k_{0}^{2}=-k^{2}\cos^{2}\psi_{1},~{}~{}~{}r_{0}^{2}\to-
r_{0}^{2}=-r^{2}\cos^{2}\psi_{2}.$ (5)

As a result of all transformations, two integrals over k and r are factorized,
and the contribution to the interaction operator in momentum space can be
represented as follows:

$\Delta V^{c}=-\frac{64}{9}\frac{\alpha^{2}(Z\alpha)^{2}}{{\bf
t}^{2}+M_{A}^{2}}\int\frac{d^{4}k}{\pi^{2}}A(t^{2},k^{2},k^{2})\frac{(2k^{2}+k_{0}^{2})}{k^{2}(k^{2}-2m_{1}k_{0})}\int\frac{d^{4}r}{\pi^{2}}A(t^{2},r^{2},r^{2})\frac{(2r^{2}+r_{0}^{2})}{r^{2}(r^{2}-2m_{2}r_{0})}$
(6)

To calculate each of the integrals, it is necessary to know the form of the
transition form factor $A(t^{2},k^{2},k^{2})$ of $1^{++}$ meson to two
photons, which is one of the main structural elements of the formula (6). At
present we have only few experimental data on it [33, 34, 35]. The L3
Collaboration studied the reaction $e^{+}e^{-}\to
e^{+}e^{-}\gamma^{\ast}\gamma^{\ast}\to e^{+}e^{-}f_{1}(1285)\to
e^{+}e^{-}\eta\pi^{+}\pi^{-}$ in [33] and measured the $f_{1}(1285)$
transition form factor for the case when one of the photons is real and
another one is virtual. In [34] the production of $f_{1}(1420)$ was
investigated by the L3 Collaboration in the reaction
$\gamma^{\ast}\gamma^{\ast}\rightarrow K_{S}^{0}K^{\pm}\pi^{\mp}$. By using
these data, we can parameterize the transition form factor for the case of two
photons with equal virtualities as follows:

$A(M_{A}^{2},k^{2},k^{2})=A(M_{A}^{2},0,0)F_{AV}(k^{2}),~{}~{}~{}F_{AV}(k^{2})=\frac{\Lambda_{A}^{4}}{(\Lambda_{A}^{2}-k^{2})^{2}}.$
(7)

The effects of off-shellness for exchange by massive $f_{1}$ mesons might be
important. This effect was investigated in [36, 37], and in [38] a simple
parametrization was proposed. The simplest way to take it into account is the
introduction of the exponential suppression factor [38]:

$\frac{A(t^{2},0,0)}{A(M_{A}^{2},0,0)}\approx
e^{(t^{2}-M_{A}^{2})/M_{A}^{2}},$ (8)

which gives the factor $\sim e^{-1}$ for $t^{2}\approx 0$. The values of the
form factors in (7) for the case of $f_{1}(1285)$ and $f_{1}(1420)$ can be
fixed from L3 data [16]:

$\displaystyle
A_{f_{1}(1285)\gamma^{\ast}\gamma^{\ast}}\left(M_{f_{1}(1285)}^{2},0,0\right)=(0.266\pm
0.043)~{}\mathrm{GeV}^{-2},$ $\displaystyle
A_{f_{1}(1260)\gamma^{\ast}\gamma^{\ast}}\left(M_{f_{1}(1260)}^{2},0,0\right)=(0.160\pm
0.120)~{}\mathrm{GeV}^{-2},$ $\displaystyle
A_{f_{1}(1420)\gamma^{\ast}\gamma^{\ast}}\left(M_{f_{1}(1420)}^{2},0,0\right)=(0.193\pm
0.041)~{}\mathrm{GeV}^{-2}.$ (9)

Table 1:  Hadronic light-by-light contribution to muonium HFS. The top line in
each cell corresponds to a vertical exchange, and the bottom line corresponds
to a horizontal exchange.
Meson | Mass | $I^{G}(J^{PC})$ | $\Lambda$ | $A(M^{2},0,0)$ | $\Delta E^{hfs}(1S)$
---|---|---|---|---|---
| in MeV |  | in MeV |  | in Hz
$f_{1}$ | 1281.9 | $0^{+}(1^{++})$ | 1040 | 0.266 | -0.00028
|  |  |  | $GeV^{-2}$ | -0.00311
$a_{1}$ | 1260 | $1^{-}(1^{++})$ | 1040 | 0.160 | -0.00011
|  |  |  | $GeV^{-2}$ | -0.00115
$f_{1}$ | 1426.3 | $0^{+}(1^{++})$ | 926 | 0.193 | -0.00007
|  |  |  | $GeV^{-2}$ | -0.00096
$\sigma$ | 550 | $0^{+}(0^{++})$ | 2000 | -0.596 | 0
|  |  |  | $GeV^{-1}$ | 0.02701
$f_{0}$ | 980 | $0^{+}(0^{++})$ | 2000 | -0.085 | 0
|  |  |  | $GeV^{-1}$ | 0.00023
$a_{0}$ | 980 | $1^{-}(0^{++})$ | 2000 | -0.086 | 0
|  |  |  | $GeV^{-1}$ | 0.00023
$f_{0}$ | 1370 | $0^{+}(0^{++})$ | 2000 | -0.036 | 0
|  |  |  | $GeV^{-1}$ | 0.00002
$\pi^{0}$ | 134.9768 | $1^{-}(0^{-+})$ | 770 | 0.025 | 0
|  |  |  | $GeV^{-1}$ | -0.00650
$\eta$ | 547.862 | $0^{+}(0^{-+})$ | 774 | 0.024 | 0
|  |  |  | $GeV^{-1}$ | -0.00170
$\eta^{\prime}$ | 957.78 | $0^{+}(0^{-+})$ | 859 | 0.031 | 0
|  |  |  | $GeV^{-1}$ | -0.00165
$f_{2}$ | 1275.4 | $0^{+}(2^{++})$ | 2000 | 0.498 | 0
|  |  |  |  | 0.00006
Total contribution | 0.012 Hz |

Using the dipole parameterization from (7) we can calculate sequentially
analytically the integrals over all variables in the Euclidean space:

$I_{e}=\int
d^{4}k\frac{(2k^{2}+k_{0}^{2})}{k^{2}(k^{2}-2k_{0}m_{1})}\frac{\Lambda^{4}}{(k^{2}-\Lambda^{2})^{2}}=-\int\limits_{0}^{\infty}dk^{2}L_{e}(k^{2})A(0,k^{2},k^{2})=$
(10)
$-\frac{\pi^{2}\Lambda_{A}^{2}}{4(1-a_{e}^{2})^{5/2}}\left[3\sqrt{1-a_{e}^{2}}-a_{e}^{2}(5-2a_{e}^{2})\ln\frac{1+\sqrt{1-a_{e}^{2}}}{a_{e}}\right],$
$L_{e}(k^{2})=\frac{\pi^{2}}{8m_{1}^{4}}\left[k^{2}(k^{2}-6m_{1}^{2})-(k^{2}-8m_{1}^{2})\sqrt{k^{2}(k^{2}+4m_{1}^{2})}\right],~{}~{}~{}a_{e}=\frac{2m_{1}}{\Lambda}.$
(11)

The integral for the muon loop $I_{\mu}$ is obtained by replacing $m_{1}\to
m_{2}$. Thus, the final contribution to the muonium HFS can be represented by
the following analytical formula:

$\Delta
E^{hfs}_{c}(1S)=-\frac{64\alpha^{2}(Z\alpha)^{5}\mu^{3}A(0,0,0)^{2}}{9\pi
M_{A}^{2}\left(1+\frac{2W}{M_{A}}\right)^{2}}I_{e}I_{\mu}.$ (12)

For numerical estimates of this contribution, we take three axial vector
mesons with masses 1285 MeV, 1260 MeV and 1420 MeV. The total numerical value
of the contribution is presented in Table 1. We write out in the Table 1 the
numerical values of the individual contributions to the nearest five digits
after the decimal point, bearing in mind that the smallest contributions are
of this order.

Let us further consider horizontal exchanges with axial vector mesons shown in
Fig.1(a,b). In this case, the use of projection operators (2) also makes it
possible to reduce the product of various factors in the numerator to a common
trace, which can be calculated for the sum of the amplitudes in Fig.1(a,b) as

$N^{(a+b)}_{AV}=(k_{1}^{2}k_{2}^{4}+k_{2}^{2}k_{1}^{4})(2\cos\Omega+\cos\psi_{1}\cos\psi_{2})-k_{1}^{3}k_{2}^{3}(1+3\cos^{2}\Omega+\cos^{2}\psi_{1}+\cos^{2}\psi_{2}),$
(13) $\cos\Omega=\cos\psi_{1}\cos\psi_{2}+\sin\psi_{1}\sin\psi_{2}\cos\theta.$
(14)

To immediately take the sum of the amplitudes in Fig.1(a,b), we multiply the
direct amplitude by the factor $(k_{2}^{2}+2k_{2}^{0}m_{2})$, and the cross
amplitude by the factor $(k_{2}^{2}-2k_{2}^{0}m_{2})$. In addition, we have
passed to the Euclidean space of variables $k_{1}$ and $k_{2}$. After all
transformations, the contribution of horizontal exchanges to the HFS of the
spectrum will be determined by the following integral expression:

$\Delta
E^{hfs}_{AV,(a+b)}=\frac{16\alpha^{2}(Z\alpha)^{5}\mu^{3}\Lambda^{2}}{3\pi}\int_{0}^{\infty}dk_{1}\int\frac{d\Omega_{1}}{\pi^{2}}\int_{0}^{\infty}dk_{2}\int\frac{d\Omega_{2}}{\pi^{2}}\frac{A(M_{A}^{2},k_{1}^{2},k_{2}^{2})}{(k_{1}^{2}+a_{e}^{2}\cos^{2}\psi_{1})}\times$
(15)
$\frac{A(M_{A}^{2},k_{1}^{2},k_{2}^{2})}{(k_{2}^{2}+a_{\mu}^{2}\cos^{2}\psi_{2})}\frac{N^{(a+b)}_{AV}}{(k_{1}^{2}+k_{2}^{2}+2k_{1}k_{2}\cos\Omega+\frac{M_{A}^{2}}{\Lambda^{2}})},$

where $d\Omega_{1}=2\pi\sin^{2}\psi_{1}\sin\theta d\theta d\psi_{1}$,
$d\Omega_{2}=4\pi\sin^{2}\psi_{2}d\psi_{2}$. Further, the calculation of these
integrals is carried out numerically, and the results are presented in Table
1.

##  III Contribution of scalar mesons

Recent results on the properties of light scalar mesons [39] show that they
are being intensively studied, including decays into two photons. But the
accuracy of measuring the decay width $\Gamma_{S\gamma\gamma}$ is currently
not high. Let us consider the contribution of scalar mesons to the interaction
amplitudes and HFS, using the methods formulated in the previous section for
constructing hadronic light-by-light scattering amplitudes. The general
parametrization of scalar meson $\rightarrow\gamma^{\ast}+\gamma^{\ast}$
vertex function takes the form [40, 41, 42, 43]:

$T_{S}^{\mu\nu}(t,k_{1},k_{2})=4\pi\alpha\biggl{\\{}A(t^{2},k_{1}^{2},k_{2}^{2})(g^{\mu\nu}(k_{1}\cdot
k_{2})-k_{1}^{\nu}k_{2}^{\mu})+$ (16)
$B(t^{2},k_{1}^{2},k_{2}^{2})(k_{2}^{\mu}k_{1}^{2}-k_{1}^{\mu}(k_{1}\cdot
k_{2}))(k_{1}^{\nu}k_{2}^{2}-k_{2}^{\nu}(k_{1}\cdot k_{2}))\biggr{\\}},$

where $A(t^{2},k_{1}^{2},k_{2}^{2})$, $B(t^{2},k_{1}^{2},k_{2}^{2})$ are two
scalar functions on three variables, $k_{1,2}$ are four momenta of virtual
photons, t is the four-momentum of scalar meson. The first term in (16)
represents transverse photons interaction, and the second term represents
longitudinal photons interaction. In the leading order, the contribution of
the structure function $A(t^{2},k_{1}^{2},k_{2}^{2})$ is decisive. $t$ is the
four momentum of scalar meson which is equal to $(k_{1}+k_{2})$ for the
horizontal exchanges. The numerator of the sum of the horizontal exchange
amplitudes is equal to

$N^{(a+b)}_{S}=k_{1}^{2}k_{2}^{2}\cos\Omega(\cos\Omega\cos\psi_{1}\cos\psi_{2}-1-\cos^{2}\psi_{1}-\cos^{2}\psi_{2}-\cos^{2}\Omega).$
(17)

The total contribution of scalar mesons to the hyperfine structure is similar
to expression (15) and in euclidean space has the following integral form:

$\Delta
E^{hfs}_{S,(a+b)}=\frac{16\alpha^{2}(Z\alpha)^{5}\mu^{3}}{3\pi}\int_{0}^{\infty}dk_{1}\int\frac{d\Omega_{1}}{\pi^{2}}\int_{0}^{\infty}dk_{2}\int\frac{d\Omega_{2}}{\pi^{2}}\frac{A(M_{S}^{2},k_{1}^{2},k_{2}^{2})}{(k_{1}^{2}+a_{e}^{2}\cos^{2}\psi_{1})}\times$
(18)
$\frac{A(M_{S}^{2},k_{1}^{2},k_{2}^{2})}{(k_{2}^{2}+a_{\mu}^{2}\cos^{2}\psi_{2})}\frac{N^{(a+b)}_{S}}{(k_{1}^{2}+k_{2}^{2}+2k_{1}k_{2}\cos\Omega+\frac{M_{S}^{2}}{\Lambda^{2}})},$

where the parameterization of function $A(M_{S}^{2},k_{1}^{2},k_{2}^{2})$ for
scalar meson chosen in the same form as for axial-vector mesons (monopole form
for variables $k_{1}^{2}$ and $k^{2}_{2}$):

$A(M_{s}^{2},k_{1}^{2},k_{2}^{2})=A_{S}\frac{\Lambda^{4}}{(k_{1}^{2}-\Lambda^{2})(k_{2}^{2}-\Lambda^{2})}.$
(19)

where the $S\gamma\gamma$ coupling $A_{S}$ is related with the
$S\to\gamma\gamma$ partial width [18, 43, 44]:

$A_{S}=\sqrt{\frac{4\Gamma_{S\gamma\gamma}}{\pi\alpha^{2}M_{S}^{3}}}.$ (20)

The vertical exchange amplitudes for scalar mesons are also constructed. The
structure of the interaction vertices is such that the vertical exchanges are
suppressed in comparison with the horizontal ones by the degree of momentum
$|{\bf t}|\sim\mu\alpha$ and therefore give a contribution of a higher order
in $\alpha$, which we omit below.

##  IV Contribution of pseudoscalar mesons

The transition vertex of two virtual photons into pseudoscalar meson is
determined only by one structure function. The effective interaction vertex of
the $\pi^{0}$ meson (or other pseudoscalar mesons $\eta$, $\eta^{\prime}$) and
virtual photons can be expressed in terms of the transition form factor
$F_{\pi^{0}\gamma^{\ast}\gamma^{\ast}}(k_{1}^{2},k_{2}^{2})$ in the form:

$V^{\mu\nu}(k_{1},k_{2})=i\varepsilon^{\mu\nu\alpha\beta}k_{1\alpha}k_{2\beta}A(t^{2},0,0)F_{\pi^{0}\gamma^{\ast}\gamma^{\ast}}(k_{1}^{2},k_{2}^{2}),~{}~{}~{}A(M_{P}^{2},0,0)=\frac{\alpha}{\pi
F_{P}},$ (21)

where the pseudoscalar meson decay constants are $F_{\pi}=0.0924$ GeV,
$F_{\eta}=0.0975$ GeV, $F_{\eta^{\prime}}=0.0744$ GeV. The pseudoscalar decay
constants are related to the two photon partial width of the resonance by the
equation:

$F^{2}_{P}=\frac{\alpha^{2}}{64\pi^{3}}\frac{M_{P}^{3}}{\Gamma(P\to\gamma\gamma)}.$
(22)

The transition form factor
$F_{\pi^{0}\gamma^{\ast}\gamma^{\ast}}(k_{1}^{2},k_{2}^{2})$ is normalized by
the condition: $F_{\pi^{0}\gamma^{\ast}\gamma^{\ast}}(0,0)=1$ and
parametrization (19) is used for it. The form factors of the transition of
pseudoscalar mesons into two photons have been studied experimentally by
various collaborations [46, 47, 48]. Fitting the experimental data using
function (9) gave the following values of the cutoff parameter:
$\Lambda_{\pi}=0.770$ GeV, $\Lambda_{\eta}=0.774$ GeV,
$\Lambda_{\eta^{\prime}}=0.859$ GeV. From the theoretical viewpoint
$F_{\eta}$, $F_{\eta^{\prime}}$ entering in (22) should be considered as
effective decay constants due to $\eta-\eta^{\prime}$ mixing [49].

The general formula that determines the contribution to the ground state HFS
from the horizontal exchange amplitudes can be represented in integral form in
Euclidean space:

$\Delta E^{hfs}_{PS}=-\frac{\alpha^{2}(Z\alpha)^{5}\mu^{3}}{3\pi
F_{\pi}^{2}}\int\frac{d^{4}k_{1}}{k_{1}\pi^{4}}\int\frac{d^{4}k_{2}}{k_{2}\pi^{4}}[F_{\pi^{0}\gamma^{\ast}\gamma^{\ast}}(k_{1}^{2},k_{2}^{2})]^{2}\times$
(23)
$\frac{N^{(a+b)}_{PS}}{(k_{1}^{2}+a_{e}^{2}\cos\psi_{1}^{2})(k_{2}^{2}+a_{\mu}^{2}\cos\psi_{2}^{2})((k_{1}+k_{2})^{2}+\frac{M_{P}^{2}}{\Lambda^{2}})},$

where the function in the numerator

$N^{(a+b)}_{PS}=\left(\cos\Omega+\cos^{2}\Omega-\cos^{3}\Omega+\cos\psi_{1}\cos\psi_{2}-\cos\Omega\cos^{2}\psi_{1}-\cos\Omega\cos^{2}\psi_{2}\right)$
(24)

is obtained by calculating the trace, summing over the Lorentz indices in an
expression like

$\varepsilon^{\mu\nu\alpha\beta}k_{1\alpha}k_{2\beta}\varepsilon^{\lambda\sigma\rho\omega}k_{1\rho}k_{2\omega}Tr[\gamma^{\lambda}(p_{1}-k_{1}+m_{1})\gamma^{\mu}\hat{\Pi}\gamma^{\nu}(-p_{2}+k_{2}+m_{2})\gamma^{\sigma}\hat{\Pi}^{+}].$
(25)

The index $(a+b)$ denotes the contribution of the diagrams $(a)$ and $(b)$ in
Fig. 1.

Subsequently, integrals in (23) are calculated numerically, as in the case of
scalar mesons with horizontal exchanges.

Turning to the vertical exchange amplitudes, it should be noted that they
contain additional powers of the momentum t. So, for example, the numerator of
the amplitude in Fig. 1(c) is equal to

$N^{c}_{PS}=t^{2}k^{2}r^{2}-k^{2}(rt)^{2}-r^{2}(kt)^{2}+(kr)(kt)(rt)-k_{0}r_{0}(kt)(rt).$
(26)

As a result, it turns out that vertical interaction contribution to the
hyperfine splitting is of order $\alpha^{2}(Z\alpha)^{7}$. Therefore, this
contribution can be neglected.

##  V Contribution of tensor mesons

The lowest tensor resonance is the spin 2 $f_{2}(1270)$ dominating in
$\gamma\gamma\to\pi^{+}\pi^{-},\pi^{0}\pi^{0}$ production. The $f_{2}$
parameters extracted are $M_{f_{2}}=1275.4$ MeV, $\Gamma_{f_{2}}=185.8$ MeV
and $\Gamma_{f_{2}\gamma\gamma}/\Gamma_{f_{2}}={1.42\pm 0.24}\times 10^{-5}$.
For tensor mesons consisting from light quarks the experimental analysis of
decay angular distributions for $\gamma\gamma$ cross sections to
$\pi^{+}\pi^{-}$, $\pi^{0}\pi^{0}$, $K^{+}K^{-}$ have shown that the $J=2$
mesons are produced mainly in a state with helicity $\Lambda=2$ [45]. We will
assume further that hadronic light-by-light scattering amplitude for tensor
mesons is dominated be helicity $\Lambda=2$ exchange. Then the amplitude of
the process $\gamma^{\ast}+\gamma^{\ast}\to T$ (see Fig. 1) can be
parametrised as follows [40]:

$T^{T}_{\mu\nu\alpha\beta}(k_{1},k_{2})=4\pi\alpha\frac{k_{1}k_{2}}{M}{\cal
M}_{\mu\nu\alpha\beta}(k_{1},k_{2}){\cal
F}_{T\gamma^{\ast}\gamma^{\ast}}(k_{1}^{2},k_{2}^{2}),$ (27)

where ${\cal F}_{T\gamma^{\ast}\gamma^{\ast}}(k_{1}^{2},k_{2}^{2})$ is a
transition form factor, $k_{1}$, $k_{2}$ are four momenta of virtual photons,

${\cal
M}_{\mu\nu\alpha\beta}(k_{1},k_{2})=\biggl{\\{}R_{\mu\alpha}(k_{1},k_{2})R_{\nu\beta}(k_{1},k_{2})+\frac{1}{8(k_{1}+k_{2})^{2}\left[(k_{1}k_{2})^{2}-k_{1}^{2}k_{2}^{2}\right]}R_{\mu\nu}(k_{1},k_{2})\times$
(28)
$\left[(k_{1}+k_{2})^{2}(k_{1}-k_{2})_{\alpha}-(k_{1}^{2}-k_{2}^{2})(k_{1}+k_{2})_{\alpha}\right]\times\left[(k_{1}+k_{2})^{2}(k_{1}-k_{2})_{\beta}-(k_{1}^{2}-k_{2}^{2})(k_{1}+k_{2})_{\beta}\right]\biggr{\\}},$
$R_{\mu\nu}(k_{1},k_{2})=-g_{\mu\nu}+\frac{1}{X}\left[(k_{1}k_{2})(k_{1}^{\mu}k_{2}^{\nu}+k_{2}^{\mu}k_{1}^{n}u)-k_{1}^{2}k_{2}^{\mu}k_{2}^{\nu}-k_{2}^{2}k_{1}^{\mu}k_{1}^{\nu}\right],~{}X=(k_{1}k_{2})^{2}-k_{1}^{2}k_{2}^{2}.$

Then the electron-muon direct interaction amplitude via horizontal tensor
meson exchange can be presented as follows:

$i{\cal
M}_{T}=\frac{\alpha^{2}(Z\alpha)^{2}}{16m_{1}^{2}m_{2}^{2}}\int\frac{d^{4}k_{1}}{\pi^{2}(k_{1}^{2})^{2}}\int\frac{d^{4}k_{2}}{\pi^{2}(k_{2}^{2})^{2}}\frac{(k_{1}k_{2})^{2}}{M_{T}^{2}}{\cal
F}^{2}_{T\gamma^{\ast}\gamma^{\ast}}(k_{1}^{2},k_{2}^{2}){\cal
M}_{\mu\nu\alpha\beta}(k_{1},k_{2}){\cal
M}_{\sigma\lambda\rho\omega}(k_{1},k_{2})\times$ (29)
$D_{T}^{\alpha\beta\rho\omega}(k_{1}+k_{2})Tr\Bigl{[}\hat{\Pi}(\hat{q}_{1}+m_{1})\gamma^{\sigma}S_{e}(p_{1}-k_{1})\gamma^{\mu}(\hat{p}_{1}+m_{1})\hat{\Pi}(\hat{p}_{2}-m_{2})\gamma^{\nu}S_{\mu}(-p_{2}+k_{2})\gamma^{\lambda}(\hat{q}_{2}-m_{2})\Bigr{]},$

where $S_{e}(p_{1}-k_{1})$ and $S_{\mu}(-p_{2}+k_{2})$ are the propagators of
electron and muon. The massive spin 2 propagator has the form:

$D^{\mu\nu\alpha\beta}_{T}(k)=\frac{f^{\mu\nu\alpha\beta}}{k^{2}-M_{T}^{2}+i0},$
(30)
$f^{\mu\nu\alpha\beta}=\frac{1}{2}\left(g^{\mu\alpha}g^{\nu\beta}+g^{\mu\beta}g^{nu\alpha}-g^{\mu\nu}g^{\alpha\beta}\right)+\frac{1}{2}\left(g^{\mu\alpha}\frac{k^{\nu}k^{\beta}}{M^{2}}+g^{\nu\beta}\frac{k^{\mu}k^{\alpha}}{M^{2}}+g^{\mu\beta}\frac{k^{\nu}k^{\alpha}}{M^{2}}+g^{\nu\alpha}\frac{k^{\mu}k^{\beta}}{M^{2}}\right)+$
(31)
$+\frac{2}{3}\left(\frac{1}{2}g^{\mu\nu}+\frac{k^{\mu}k^{\nu}}{M^{2}}\right)\left(\frac{1}{2}g^{\alpha\beta}+\frac{k^{\alpha}k^{\beta}}{M^{2}}\right).$

The crossed amplitude in Fig. 1(b) has the similar structure. After further
simplifications of the numerator of the expression (29) in the Form package,
it takes the following form:

$N_{T}^{(a,b)}=k_{1}^{0}k_{2}^{0}-k_{1}k_{2}+\frac{1}{k_{1}^{2}k_{2}^{2}-(k_{1}k_{2})^{2}}[k_{1}^{2}(k_{2}^{0})^{2}(k_{1}k_{2})+k_{2}^{2}(k_{1}^{0})^{2}(k_{1}k_{2})-2k_{1}^{2}k_{2}^{2}k_{1}^{0}k_{2}^{0}].$
(32)

For the form factor of the transition of a tensor meson into two virtual
photons, we use the monopole parametrization with respect to each square of
the photon momentum of the form:

${\cal
F}_{T\gamma^{\ast}\gamma^{\ast}}(k_{1}^{2},k_{2}^{2})=\frac{A_{T\gamma^{\ast}\gamma^{\ast}}(M_{T}^{2},0,0)}{(k_{1}^{2}-\Lambda_{T}^{2})(k_{2}^{2}-\Lambda_{T}^{2})},$
(33)

and the value of $A_{T\gamma^{\ast}\gamma^{\ast}}(M_{T}^{2},0,0)$ is
determined using the width $\Gamma_{T~{}\gamma\gamma}$ of the decay of the
tensor meson into two photons:

$A_{T\gamma^{\ast}\gamma^{\ast}}(M_{T}^{2},0,0)=\frac{2\sqrt{5\Gamma_{T~{}\gamma\gamma}}}{\alpha\sqrt{\pi
M_{T}}}.$ (34)

The numerical values of the decay width of one tensor meson $f_{2}(1275)$ is
taken from [39]. After passing to the Euclidean space and a number of
simplifications, we can represent the total contribution in Fig. 1(a+b) to the
hyperfine structure of muonium in the integral form:

$\Delta
E^{hfs}_{T}=\frac{128\pi\alpha^{2}(Z\alpha)^{5}\mu^{3}}{3M_{T}^{2}}\int_{0}^{\infty}k_{1}^{2}dk_{1}\int_{0}^{\pi}\frac{\sin^{2}\psi_{1}}{\pi^{3}}\int_{0}^{\infty}k_{2}^{2}dk_{2}\int_{0}^{\pi}\frac{\sin^{2}\psi_{2}}{\pi^{3}}\int_{0}^{\pi}\sin\theta
d\theta$ (35)
$\frac{A^{2}_{T\gamma^{\ast}\gamma^{\ast}}(M_{T}^{2},0,0)N_{T}^{(a,b)}}{(k_{1}^{2}+1)^{2}(k_{2}^{2}+1)^{2}(k_{1}^{2}+a_{e}^{2}\cos^{2}\psi_{1})(k_{2}^{2}+a_{\mu}^{2}\cos^{2}\psi_{2})[(k_{1}+k_{2})^{2}+\frac{M_{T}^{2}}{\Lambda^{2}}]}.$

The results of the numerical calculation (35) are presented in the Table 1
only for one tensor meson since the contribution of other mesons is negligible
due to the small width $\Gamma_{T~{}\gamma\gamma}$.

##  VI Conclusion

As is known, the last measurement of the hyperfine splitting of the ground
state in muonium was carried out in 1999 [11] with a record-breaking accuracy
for those times up to hundredths of a kHz. In a recent paper [10], the MuSEUM
collaboration announced the start of new measurements of HFS in muonium and
obtained a result that agrees with [11], but is still inferior to it in
accuracy. It can be said that the planned increase in the accuracy of
measuring HFS in muonium to 1 ppb [13] opens a new stage in the theoretical
study of this problem, which is connected with an increase in the accuracy of
calculations of various corrections. It should be emphasized that theoretical
work in this direction did not stop during the last two decades [1, 2, 3, 12,
4]. Various high-order quantum electrodynamic contributions in $\alpha$ were
calculated. A sharp increase in the experimental accuracy leads to the need to
take into account in the theoretical calculations the contributions of other
interactions, as is the case for the anomalous magnetic moment of the muon or
the Lamb shift in muonic hydrogen. This work is devoted to the study of one of
these new contributions, due to the production of hadrons in light-by-light
scattering diagrams.

Compared to our previous work in 2002, this study takes into account the
contributions of light mesons of different spins both in horizontal-type
diagrams and in amplitudes with vertical exchange. The calculated
contributions from various mesons are presented separately in the final Table
1. For all mesons, the parameter $A(M^{2},0,0)$, which is also presented in
the Table 1, plays an important role in the numerical evaluation of the
contribution. The numerical value of this parameter is related to the width of
the meson decay into two photons, which is taken from various experiments. An
analysis of the available experimental data on the decay widths into two
photons shows that the accuracy of their measurement is not high [39].
Therefore, it is more correct to consider the results presented in Table 1 as
possible estimates of contributions of this type.

For mesons for which the $\Gamma_{\gamma\gamma}$ value has not yet been fixed
[39], the average values were taken from the available data. But with
pseudovector and pseudoscalar mesons, which make the main contribution to the
Table 1, the situation with fixing $A(M^{2},0,0)$ is more or less certain, so
that the error of their contribution obtained does not exceed 30 percent.
Nevertheless, there is a significant scatter in the experimental data for the
width of the $\sigma$ meson $\Gamma_{\sigma\gamma\gamma}$. In our
calculations, we used for it a value of 4.5 keV. Since in the end it turns out
that the contribution of this meson is the main one, we estimate the total
error of the calculation in the Table 1 at 50 percent.

It should be noted that the obtained contributions of pseudoscalar mesons
improve our results due to more accurate numerical integration. The
calculation formula (23) is transformed in comparison with [14] in such a way
that the contribution of both direct and cross horizontal amplitudes is taken
into account at once. Numerically, the contributions of $\pi$, $\eta$,
$\eta^{\prime}$ mesons are among the most significant. A close numerical value
is given by the horizontal exchange amplitude with the $\sigma$ meson. As
regards the contribution of axial vector mesons, they contribute from both
types of exchanges (horizontal and vertical). The difference between our
results on vertical exchanges of pseudovector mesons and work [15] is, in our
opinion, that we take into account an additional reducing factor (8), the
square of which just leads to a decrease in our contribution compared to [15]
by an order of magnitude.

The total contribution of all mesons to the hyperfine splitting turned out to
be positive. Although the contributions of axial vector and pseudoscalar
mesons are negative, there is a positive contribution of the $\sigma$ meson,
which exceeds all previous ones in magnitude. The resulting value of 0.012 Hz
can be regarded as an estimate of this small hadronic effect.

###### Acknowledgements.

The authors are grateful to A.E. Radzhabov and A.S. Zhevlakov for useful
discussions. This work is supported by Russian Science Foundation (grant No.
RSF 23-22-00143).

## References

  * [1] M. I. Eides, H. Grotch and V. A. Shelyuto, Phys. Rept. 342, 63 (2001).
  * [2] M. I. Eides and V. A. Shelyuto, arXiv:2306.07748 [hep-ph].
  * [3] S. G. Karshenboim and E. Yu. Korzinin, Phys. Rev. A 103, 2, 022805 (2021).
  * [4] M. I. Eides and V. A. Shelyuto, Phys. Rev. D 92, 1, 013010 (2015).
  * [5] F. Jegerlehner, The Anomalous Magnetic Moment of the Muon, Springer Tracts in Modern Physics, V.274, 2nd edition: Springer International Publishing AG 2017.
  * [6] A. E. Radzhabov, A. S. Zhevlakov, A. P. Martynenko, and F. A. Martynenko, Phys. Rev. D 108, 014033 (2023).
  * [7] F. V. Ignatov et al. [CMD-3 Collaboration] e-Print: 2302.08834 [hep-ex].
  * [8] P. Crivelli, Hyperfine Interact. 239 (1), 49 (2018).
  * [9] B. Ohayon, G. Janka, I. Cortinovis et al., Phys. Rev. Lett. 128, 011802 (2022).
  * [10] S. Kanda, Y. Fukao, Y. Ikedo, Phys. Lett. B 815 136154 (2021).
  * [11] W. Liu et al., Phys. Rev. Lett. 82, 711 (1999).
  * [12] M. Eides, Phys. Lett B 795, 113 (2019).
  * [13] R. Iwai, M. Abe, S. Fukumura et al. Jour. of Phys.: Conf. Series 2462, 012019 (2023).
  * [14] R. N. Faustov and A. P. Martynenko, Phys. Lett. B 541, 135 (2002).
  * [15] S. G. Karshenboim, V. A. Shelyuto, and A. I. Vainshtein, Phys. Rev. D 78, 065036 (2008).
  * [16] A. E. Dorokhov, N. I. Kochelev, A. P. Martynenko, F. A. Martynenko, and A. E. Radzhabov, Phys. Lett. B 776, 105 (2018).
  * [17] A. E. Dorokhov, N. I. Kochelev, A. P. Martynenko, F. A. Martynenko, and R. N. Faustov, Phys. Part. Nucl. Lett. 14, 857 (2017).
  * [18] A. E. Dorokhov, A. P. Martynenko, F. A. Martynenko, and A. E. Radzhabov, Phys. Part. Nucl. Lett. 16, 5, 520 (2019).
  * [19] A. E. Dorokhov, A. P. Martynenko, F. A. Martynenko, and A. E. Radzhabov, EPJ Web Conf. 222, 03010 (2019).
  * [20] A. E. Dorokhov, A. P. Martynenko, F. A. Martynenko, and A. E. Radzhabov, J. Phys. Conf. Ser. 1435, 1, 012004 (2020).
  * [21] A. Miranda, P. Roig, and P. Sanchez-Puertas, Phys. Rev. D 105, 016017 (2022).
  * [22] F. Hagelstein and V. Pascalutsa, PoS CD15, 077 (2016).
  * [23] H.-Q. Zhou and H.-R. Pang, Phys. Rev. A 92, 032512 (2015), [Erratum: Phys. Rev. A 93, 069903 (2016)].
  * [24] N. T. Huong, E. Kou, and B. Moussallam, Phys. Rev. D 93, 114005 (2016).
  * [25] A. Antognini, F. Kottmann, F. Biraben et al., Ann. Phys. 331, 127 (2013).
  * [26] J. Krauth, M. Diepold, B. Franke et al., Ann. Phys. 366, 168 (2016).
  * [27] M. Diepold, B. Franke, J. Krauth et al., Ann. Phys. 396, 220 (2018).
  * [28] J. J. Krauth, K. Schuhmann, M. A. Ahmed et al., Nature 589, 527 (2021).
  * [29] A. Antognini, F. Hagelstein, and V. Pascalutsa, Annu. Rev. Nucl. Part. Sci. 72, 389 (2022).
  * [30] C. Pizzolotto, A. Adamczak, D. Bakalov, G. Baldazzi, M. Baruzzo, Eur. Phys. J. A 56, 7, 185 (2020).
  * [31] M. Stoilov, A. Adamczak, D. Bakalov, et al., Phys. Rev. A 107, 032823 (2023).
  * [32] R. N. Cahn, Phys. Rev. D 35, 3342 (1987).
  * [33] P. Achard et al. [L3 Collaboration], Phys. Lett. B 526, 269 (2002).
  * [34] P. Achard et al. [L3 Collaboration], JHEP 03, 018 (2007)
  * [35] H. Aihara et al., Phys. Rev. D 38, 1 (1988).
  * [36] A. E. Dorokhov, A. E. Radzhabov, and A. S. Zhevlakov, Eur. Phys. J. C 71, 1702 (2011).
  * [37] A. E. Dorokhov, A. E. Radzhabov, and A. S. Zhevlakov, Eur. Phys. J. C 72, 2227 (2012).
  * [38] P. Lebiedowicz and A. Szczurek, Phys. Lett. B 772, 330 (2017).
  * [39] R. L. Workman et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2022, 083C01 (2022).
  * [40] V. Pauk and M. Vanderhaeghen, EPJ C 74, 3008 (2014).
  * [41] H. Q. Zhou, Phys. Rev. C 95, 025203 (2017).
  * [42] D. Borisyuk, Phys. Rev. C 96, 055201 (2017).
  * [43] M. K. Volkov, E. A. Kuraev, and Yu. M. Bystritskiy, Phys. Atom. Nucl. 73, 443 (2010).
  * [44] F. Giacosa, Th. Gutsche, and V. E. Lyubovitskij, Phys. Rev. D 77, 034007 (2008).
  * [45] M. R. Pennington et al., Eur. Phys. J. C 56, 1 (2008).
  * [46] J. Gronberg et al. (CLEO Collaboration) Phys. Rev. D 57, 33 (1998).
  * [47] J. P. Lees et al. (BaBar Collaboration) Phys. Rev. D 98, 112002 (2018).
  * [48] S. Uehara et al. (Belle Collaboration) Phys. Rev. D 86, 092007 (2012).
  * [49] Th. Feldmann, P. Kroll, and B. Stech, Phys. Rev. D 58, 114006 (1998).


# OTFS-NOMA System for MIMO Communication Networks with Spatial Diversity

Wafa Hedhly, Leila Musavian, and Nikolaos Thomos The authors are with the
School of Computer Science and Electronic Engineering, University of Essex,
Wivenhoe Park, Colchester CO4 3SQ, United Kingdom (e-mail: {wafa.hedhly,
leila.musavian, nthomos}@essex.ac.uk). This work was funded by the Engineering
and Physical Science Research Council, UK, [grant number EP/X012204/1].

###### Abstract

In this work, we study the use of non-orthogonal multiple access (NOMA) and
orthogonal time frequency space (OTFS) modulation in a multiple-input
multiple-output (MIMO) communication network where mobile users (MUs) with
different mobility profiles are grouped into clusters. We consider a downlink
scenario where a base station (BS) communicates with multiple users that have
diverse mobility profiles. High-mobility (HM) users’ signals are placed in the
delay-Doppler (DD) domain using OTFS modulation in order to transform their
time-varying channel into a sparse static channel, while low-mobility (LM)
users signals are placed in the time-frequency (TF) domain. Precoding is
adopted at the BS to direct focused beams towards each cluster of users.
Moreover, NOMA spectrum sharing is used in each cluster to allow the
coexistence of a single HM user and multiple LM users within the same resource
block. LM users access disjoint subchannels to ensure their orthogonality. All
users within the same cluster first detect the HM user’s signal. Afterward, LM
users suppress the interference from the HM user and detect their own signals.
Closed-form expressions of the detection signal-to-noise ratios (SNRs) are
derived. The numerical results showed that the performance of the proposed
system highly depends on the number of LM users, the number of clusters and
the power allocation factors between HM and LM users.

###### Index Terms:

OTFS modulation, NOMA spectrum sharing, MIMO precoding, delay-Doppler,
performance analysis.

##  I Introduction

Current and future wireless communication networks are characterized by
heterogeneity on many fronts, incorporating diverse technologies and protocols
and serving users with different performance requirements, channel conditions,
and mobility profiles. One way of catering to this diversity is to implement
NOMA, a multiple-access paradigm that promotes user fairness and increases
spectral efficiency by simultaneously sharing the resource block between users
that are distinguished by their channel conditions. On this account, multiple
users can be accommodated in the same resource block while mitigating the
resulting interference through successive interference cancellation at the
receiver [1]. On the other hand, the high mobility of users results in severe
Doppler conditions and doubly dispersive channels, where the orthogonal
frequency division multiplexing (OFDM) modulation fails to achieve acceptable
performance [2]. Along this line, OTFS lately emerged as a new modulation
technique where information symbols are embedded in the DD domain, yielding a
sparse representation of the time-varying wireless channel. Motivated by these
compelling attributes, we aim to study the promise of the parallel
implementation of NOMA and OTFS modulation.

Several research works have investigated the potential of OTFS-NOMA in
improving the performance of communication networks with heterogeneous
mobility profiles [3, 4, 5, 6]. In [3], an OTFS-NOMA transmission protocol is
proposed with different mobility profiles, where a single HM user is served in
the DD domain and modulated using OTFS, and LM users are modulated using OFDM.
In [4], a design of MIMO beamforming is proposed for a downlink OTFS-NOMA
scenario where a HM user and multiple LM users are served by a BS equipped
with multiple antennas. The proposed setups in [3] and [4] are constrained by
the presence of one HM user. An uplink communication scenario was adopted in
[7] where NOMA access is implemented in the DD domain in order to accommodate
stationary and mobile users using OTFS modulation for all users. However, OTFS
modulation is suitable for users with high mobility profiles and high Doppler
spreads. Differently, the work in [6] considered using OTFS modulation at all
terminals in a NOMA setup for low earth orbit (LEO) constellation-enabled IoT.
Similarly, in [5], OTFS-NOMA was implemented in LEO constellations with
massive MIMO IoT terminals and a successive active user identification was
proposed. The above methods highlighted the potential of OTFS-NOMA systems,
but to the best of our knowledge, most studies consider system setups with
particular conditions, such as the use of single antennas or limiting the
number of HM users to one. Hence, the derived conclusions cannot be
generalized. Motivated by the aforementioned lack of systematic study of OTFS-
NOMA in heterogeneous networks, the objective of this work is to leverage the
merits of OTFS-NOMA systems and MIMO precoding to convey reliable
communication to multiple users with heterogeneous mobility profiles.

Thereby, we consider implementing MIMO precoding in the OTFS-NOMA system to
achieve spatial diversity and serve multiple HM and LM users with the required
performance. We propose a suitable system configuration where a BS
communicates with multiple HM and LM users grouped in distinct clusters. The
BS performs spatial MIMO precoding in order to direct focused beams towards
each cluster, where HM and LM users access the resources using NOMA paradigm.
LM users access the spectrum using disjoint subchannels to ensure
orthogonality among each other and, thus, fairness among users. Differently,
HM users’ signals are modulated using OTFS. All users within the same cluster
detect the HM user’s signal in the DD domain and cancel the resulting
interference from their signals. We analyze the proposed framework and derive
closed-form expressions of the signal-to-noise ratios (SNRs) of user
detection. To thoroughly understand the performance of our system, we evaluate
it following outage probability analysis examining different numbers of users
per cluster, different numbers of clusters and various power allocation
factors.

##  II System Model

Figure 1: Proposed system setup.

In this paper, we consider a communication network as depicted in Fig. 1,
where a BS communicates with multiple MUs located within the cell. OTFS-NOMA
is implemented in a downlink (DL) scenario. The BS is equipped with a $W\times
W$ square uniform planar array (UPA) with $A=W^{2}$ antenna elements, while
the users are equipped with a single antenna. The MUs have different mobility
profiles, i.e., low or high mobility. We assume multiple HM users exist in the
cell and can be distinguished based on their angle-of-departures (AoDs).
Moreover, we assume multiple LM users are distributed around each HM user. The
HM user and LM users in each cluster share the spectrum using NOMA principle.
The users are grouped into $Q$ clusters based on their location. Each cluster
$q$ comprises $U_{q}+1,1\leq q\leq Q$ users where one user ($u=0$) is
considered a HM user and the rest $U_{q}$ users are LM. We consider that the
maximum number of LM users is the available subchannels, i.e., $U_{q}\leq
M,\forall\;1\leq q\leq Q$, where $M$ is the number of subchannels. The BS
transmits focused beams towards each cluster of users. The precoding at the BS
is implemented based on the location of the HM user of each cluster. Each user
$u$ of cluster $q$ is located at the elevation angle $\theta_{u,q}$ and
azimuth angle $\phi_{u,q}$ in the spherical coordinate system as depicted in
Fig. 1.

The HM user receives OTFS-modulated signals in the DD domain, while LM users’
signals are placed in the TF domain (the two domains are non-orthogonal). The
BS implements NOMA access technique and superimposes the HM user signal
modulated in the DD domain with the LM users’ signals modulated in the TF
domain. The HM user accesses the whole resource block while LM users occupy
disjoint subchannels similarly to [3]. This access technique ensures that LM
users do not have frequency-selective channels and are orthogonal amongst each
other. Each MU shares its own channel station information (CSI) with the BS.
All users within the same cluster first detect the HM user’s signal of the
same cluster. Then, the LM users subtract the HM user’s signal from their
received signal. Last, each LM user detects its own signal.

Let us denote the transmit SNR at the BS as $\rho_{\mathrm{T}}$. We assume
that the BS transmits equal power towards all LM users within the same
cluster, while each cluster receives a total power based on the location of
its HM user with respect to the BS.

###  II-A Transmitted Signals

The BS transmits an OTFS frame to each HM user with $NM$ symbols in the DD
domain, where $N$ is the number of Doppler bins and $M$ is the number of delay
bins drawn from a modulation alphabet. Each transmitted packet has a duration
of $NT$ and occupies a bandwidth of $M\Delta f$, where $T$ is the sampling
time interval and $\Delta f=1/T$ is the spacing between subcarriers. The DD
signal is converted to the TF domain using the inverse symplectic finite
Fourier transform (ISFFT) as follows,

$s_{0,q}^{\mathrm{TF}}[n,m]=\frac{1}{NM}\sum_{k=0}^{N-1}\sum_{l=0}^{M-1}s_{0,q}^{\mathrm{DD}}[k,l]e^{j2\pi\left(\frac{kn}{N}-\frac{ml}{M}\right)},$
(1)

where $s_{0,q}^{\mathrm{DD}}[k,l]$ is the baseband transmitted binary phase-
shift keying symbol at the $k^{\mathrm{th}}$ Doppler and $l^{\mathrm{th}}$
delay bin for $1\leq n\leq N$ and $1\leq m\leq M$. For all $0\leq n\leq N-1$,
the TF signal of LM user $u$ of cluster $q$ is expressed as,

$\displaystyle\mathbf{s}_{u,q}^{\mathrm{TF}}[n,m]=\begin{cases}&\mathbf{s}_{u,q}^{\mathrm{TF}}[n,u-1],\quad\mathrm{if}\;m=u-1,\\\
&0,\quad\mathrm{otherwise}.\end{cases}$ (2)

The users in each cluster access the spectrum using NOMA principle. Therefore,
the transmitted signal to cluster $q$ is expressed in the DD domain as,

${s}_{q}[k,l]=\sum_{u=0}^{U_{q}}{s}_{u,q}^{\mathrm{DD}}[k,l],$ (3)

where ${s}_{u,q}^{\mathrm{DD}}$ is the information bearing signal for user $u$
in cluster $q$ expressed in the DD domain. The signal transmitted from each
antenna element $(i,i^{\prime}),0\leq i,i^{\prime}\leq W-1$, is expressed as,

${x}_{i,i^{\prime}}[k,l]=\sum_{q=1}^{Q}p_{i,i^{\prime},q}{s}_{q}[k,l],$ (4)

where $p_{i,i^{\prime},q}$ is the precoding coefficient, $0\leq
i,i^{\prime}\leq W-1,1\leq q\leq Q$. The DD channel impulse response for user
$u$ in cluster $q$ over $P_{u,q}$ channel paths is expressed as
follows111Vectors are defined in boldface letters.,

$\mathbf{h}_{u,q}(\tau,\nu)=\sum_{p=1}^{P_{u,q}}h_{p,u}^{q}\mathbf{v}\left(\theta_{p,u}^{q},\phi_{p,u}^{q}\right)\delta(\tau-\tau_{p,u}^{q})\delta(\nu-\nu_{p,u}^{q}),$
(5)

where $h_{p,u}^{q}$, $\tau_{p,u}^{q}$ and $\nu_{p,u}^{q}$ are the channel
coefficient, delay and Doppler of the $p^{\mathrm{th}}$ path, respectively and
$\mathbf{v}$ is the steering vector depending on the elevation angle
$\theta_{p,u}^{q}$ and azimuth angle $\phi_{p,u}^{q}$. The UPA is placed on
the $(x,z)$ plane of the cartesian coordinates. Therefore, the
$(i,i^{\prime})^{\mathrm{th}}$ element of the vector $\mathrm{v}$, where
$0\leq i,i^{\prime}\leq W-1$, is expressed as [8],

$v_{ii^{\prime}}(\theta,\phi)=\exp{\left(j2\pi\frac{d}{\lambda}\left(i\sin(\theta)\cos(\phi)+i^{\prime}\cos(\theta)\right)\right)},$
(6)

where $d$ is the separation distance between the UPA antenna elements assumed
to be equal in the two dimensions, and $\lambda$ is the wavelength. In the
case of LM users, the channel is not subject to Doppler shifts since they have
low mobility. Therefore, their Doppler shifts are
$\nu_{p,u}^{q}=0,\forall\;1\leq u\leq U_{q}$ and $1\leq q\leq Q$.

###  II-B Received Signals

Assuming bi-orthogonality between the transmit and receive waveforms, the
received signal at the HM user of cluster $q$ from antenna element
$(i,i^{\prime})$, $0\leq i,i^{\prime}\leq W-1$ is expressed in the DD domain
as [2],

$\displaystyle y_{0,q}^{i,i^{\prime}}[k,l]$
$\displaystyle=\sum_{p=1}^{P_{0,q}}h_{p,0}^{q}e^{-j2\pi\tau_{p,0}^{q}\nu_{p,0}^{q}}v_{ii^{\prime}}\left(\theta_{p,0}^{q},\phi_{p,0}^{q}\right)\times$
$\displaystyle
x_{i,i^{\prime}}[\left(k-k_{p,0}^{q}\right)_{N},\left(l-l_{p,0}^{q}\right)_{M}]+w_{0,q}[k,l],$
(7)

where $w_{0,q}$ denotes the complex Gaussian additive noise with zero mean and
variance $\sigma_{\mathrm{w}}^{2}$, $x_{i,i^{\prime}}^{\mathrm{DD}}$ is the
transmitted signal, $k_{p,0}^{q}$ and $l_{p,0}^{q}$ are the indices for the
Doppler and delay taps of path $p$, respectively222Herein, we assume the OTFS
grid resolution is high enough to neglect fractional Doppler. $(x)_{N}$
denotes $x$ modulo $N$.. The Doppler and delay taps are expressed as,

$\nu_{p,0}^{q}=\frac{k_{p,0}^{q}}{NT}\quad\mathrm{and}\quad\tau_{p,0}^{q}=\frac{l_{p,0}^{q}}{M\Delta
f}.$ (8)

Similarly to [9], we assume that the path gains $h_{p,u}^{q},1\leq p\leq
P_{u,q}$ are independent and identically distributed (i.i.d)
$\mathcal{C}(0,\sigma_{\mathrm{p}}^{2})$ with a uniform scattering profile,
i.e., $\sigma_{\mathrm{p}}^{2}=1/P_{u,q}$. We also assume that the path delay
and Doppler taps have the same locations between all transmit antennas and the
user receiver.

The received signal at user $u$ in cluster $q$ is expressed as,

$y_{u,q}=\sum_{i=1}^{W}\sum_{i^{\prime}=1}^{W}y_{u,q}^{i,i^{\prime}},$ (9)

where $y_{u,q}^{i,i^{\prime}}$ is the signal received by user $u$ from antenna
element $(i,i^{\prime})$. To simplify the notations, hereafter, we substitute
the $(i,i^{\prime})$ of each antenna element by the single index $a$, where
$1\leq a\leq A$. As a result, the received signal vector at each user $u$ in
cluster $q$ is expressed as [10],

$\mathbf{y}_{u,q}=\sum_{a=1}^{A}\mathbf{H}_{a,u,q}\mathbf{x}_{a}+\mathbf{w}_{u,q},$
(10)

where $\mathbf{H}_{a,u,q}\in\mathcal{C}^{N\\!M\times N\\!M}$ is a block
circulant channel matrix between user $u$ and antenna $a$ ($M$ circulant
blocks of $N\times N$ circulant matrices) and
$\mathbf{w}_{u,q}\sim\mathcal{CN}(\mathbf{0},{\sigma}_{\mathrm{w}}^{2}\mathbf{I}_{N\\!M})$
is the i.i.d additive noise vector in the DD domain.The signal vector
$\mathbf{x}_{a}$ is constructed from $x_{a}$ expressed in (4) while satisfying
the following condition: the $(k+Nl)^{\mathrm{th}}$ element of
$\mathbf{x}_{a}$ is equal to $x_{a}[k,l]$. The vectors $\mathbf{y}_{u,q}$ and
$\mathbf{w}_{u,q}$ are similarly constructed.

To simplify the notations, we remove the cluster index $q$ and consider the
index $u=0$ for the HM user. Without loss of generality, we consider the users
of cluster 1. Therefore, the received signal in the DD domain at user $u$ of
cluster 1 is expressed as,

$\mathbf{y}_{u}=\sum_{a=1}^{A}\mathbf{H}_{a,u}\mathbf{x}_{a}+\mathbf{w}_{u},0\leq
u\leq U_{1}.$ (11)

The received signal (11) at user $u$ can be re-written as,

$\displaystyle\mathbf{y}_{u}$
$\displaystyle=\sum_{a=1}^{A}p_{a,1}\mathbf{H}_{a,u}\mathbf{s}_{0,1}+\sum_{a=1}^{A}p_{a,1}\mathbf{H}_{a,u}\sum_{i=1}^{U_{1}}\mathbf{s}_{i,1}$
$\displaystyle+\sum_{a=1}^{A}\mathbf{H}_{a,u}\sum_{q=2}^{Q}p_{a,q}\mathbf{s}_{q}+\mathbf{w}_{u},$
(12)

where $p_{a,1}$ is the precoder for cluster 1.

##  III Users Detection

###  III-A High-Mobility User Detection

As explained earlier, all users within each cluster first detect the HM user’s
signal. The system utilizes zero-forcing (ZF) to detect the LM user at each
user of the same cluster. Since the channel matrix $\mathbf{H}_{a,u}$ is a
circulant matrix with $M$ circulant blocks, we express the OTFS channel matrix
between antenna element $a$ and user $u$ using the following decomposition
[11],

$\mathbf{H}_{a,u}=\mathbf{\Psi}^{\mathrm{H}}\mathbf{D}_{a,u}\mathbf{\Psi},$
(13)

where $\mathbf{\Psi}=\mathbf{F}_{N}\otimes\mathbf{F}_{M}$ is an $N\\!M\times
N\\!M$ matrix, $\mathbf{F}_{N}$ and $\mathbf{F}_{M}$ are, respectively,
$N\times N$ and $M\times M$ discrete Fourier transform (DFT)
matrices333$\mathbf{A}\otimes\mathbf{B}$ denotes the Kronecker product of the
matrices $\mathbf{A}$ and $\mathbf{B}$. The superscript H denotes the
Hermitian transform.. The diagonal matrix $\mathbf{D}_{a,u}$ is composed of
the eigenvalues of $\mathbf{H}_{a,u}$ and can be further expressed as,

$\mathbf{D}_{a,u}=\sum_{i=0}^{M-1}\mathbf{O}^{i}\otimes\mathbf{D}_{a,u,i},$
(14)

where
$\mathbf{O}=\mathrm{diag}\left[1,e^{j2\pi/M},\ldots,e^{j2\pi(M-1)/M}\right]$
and $\mathbf{D}_{a,u,i}$ is a diagonal matrix consisting of the eigenvalues of
the $i^{\mathrm{th}}$ circulant block of $\mathbf{H}_{a,u}$. This
decomposition significantly reduces the computation complexity of the system.
We consider the following notation,

$\mathbf{H}_{u}=\sum_{a=1}^{A}p_{a,1}\mathbf{H}_{a,u}.$ (15)

Using the diagonalization of $\mathbf{H}_{a,u}$ as shown in (13), the matrix
$\mathbf{H}_{u}$ can be re-written as,

$\mathbf{H}_{u}=\mathbf{\Psi}^{\mathrm{H}}\mathbf{D}_{u}\mathbf{\Psi},$ (16)

where the diagonal matrix $\mathbf{D}_{u}$ is expressed as,

$\mathbf{D}_{u}=\sum_{a=1}^{A}p_{a,1}\mathbf{D}_{a,u}.$ (17)

Therefore, the equalization matrix to detect $\mathbf{s}_{0,1}$ is as follows,

$\mathbf{G}_{u}=\left(\mathbf{H}_{u}^{\mathrm{H}}\mathbf{H}_{u}\right)^{-1}\mathbf{H}_{u}^{\mathrm{H}}.$
(18)

Since,
$\mathbf{\Psi}\mathbf{\Psi}^{\mathrm{H}}=\mathbf{\Psi}^{\mathrm{H}}\mathbf{\Psi}=\mathbf{I}_{MN}$,
the equalization matrix $\mathbf{G}_{u}$ can be re-expressed as,

$\mathbf{G}_{u}=\mathbf{\Psi}^{\mathrm{H}}\left(\mathbf{D}_{u}^{\mathrm{H}}\mathbf{D}_{u}\right)^{-1}\mathbf{D}_{u}^{\mathrm{H}}\mathbf{\Psi}.$
(19)

If we denote by $\lambda_{a,u,i},1\leq i\leq NM$ the eigenvalues of
$\mathbf{H}_{a,u}$, then the matrix $\mathbf{G}_{u}$ can equivalently be
written as,

$\mathbf{G}_{u}=\mathbf{\Psi}^{\mathrm{H}}\mathbf{\Delta}_{u}\mathbf{\Psi},$
(20)

where $\mathbf{\Delta}_{u}$ is a diagonal matrix consisting of elements,

$\delta_{u,i}=\frac{\sum_{a=1}^{A}p_{a,1}^{*}\lambda_{a,u,i}^{*}}{\left|\sum_{a=1}^{A}p_{a,1}\lambda_{a,u,i}\right|^{2}},1\leq
i\leq NM.$ (21)

The HM user detects its own signal while considering the LM users’ signals as
noise. LM users detect and subtract the HM user signal from the received
signal. The proposed equalizer is applied to the received signal at user $u$
as follows,

$\mathbf{G}_{u}\mathbf{y}_{u}=\mathbf{{s}}_{0,1}+\mathbf{\kappa}_{u},\quad
0\leq u\leq U_{1},$ (22)

where the interference-plus-noise term is written as,

$\mathbf{\kappa}_{u}=\sum_{i=1}^{U_{1}}\mathbf{{s}}_{i,1}+\mathbf{G}_{u}\left(\sum_{a=0}^{A}\mathbf{H}_{a,u}\sum_{q=2}^{Q}p_{a,q}\mathbf{s}_{q}+\mathbf{w}_{u}\right).$
(23)

The interference consists of intra-cluster interference from LM users’ signals
and inter-cluster interference from users in other clusters. The interference-
plus-noise term can be equivalently written as,

$\mathbf{\kappa}_{u}=\sum_{i\neq
u}^{U_{1}}\mathbf{{s}}_{i,1}+\sum_{q=2}^{Q}\mathbf{\Omega}_{q,u}\mathbf{s}_{q}+\mathbf{z}_{u},$
(24)

where the noise term is given by $\mathbf{z}_{u}=\mathbf{G}_{u}\mathbf{w}_{u}$
and the matrix $\mathbf{\Omega}_{q,u}$ is given by,

$\mathbf{\Omega}_{q,u}=\mathbf{\Psi}^{\mathrm{H}}\mathbf{\Delta}_{u}\sum_{a=1}^{A}p_{a,q}\mathbf{D}_{a,u}\mathbf{\Psi}.$
(25)

Therefore, the SINR of the $k^{\mathrm{th}}$ symbol of the HM user’s transmit
vector at user $u$ is expressed as,

$\gamma_{u}^{0}=\frac{\rho_{0}}{\sum_{i=1}^{U_{1}}\rho_{u}+\sum_{q=2}^{Q}\rho_{c,q}\omega_{q,u}+\eta_{u}},$
(26)

where $\omega_{q,u}$ and $\eta_{u}$ are written,

$\displaystyle\omega_{q,u}$
$\displaystyle=\frac{1}{NM}\sum_{i=1}^{NM}\left|\frac{\sum_{a=1}^{A}p_{a,q}\lambda_{a,u,i}}{\sum_{a=1}^{A}p_{a,1}\lambda_{a,u,i}}\right|^{2}$
$\displaystyle\eta_{u}$
$\displaystyle=\frac{1}{NM}\sum_{i=1}^{NM}\frac{1}{\left|\sum_{a=1}^{A}p_{a,1}\lambda_{a,u,i}\right|^{2}},$
(27)

and $\rho_{u}$ is the transmit SNR of the $u^{\mathrm{th}}$ user of cluster 1
and $\rho_{c,q}$ is the transmit SNR to cluster $q$. Therefore,
$\rho_{0}+\sum_{i=1}^{U_{1}}\rho_{u}=\rho_{c,1}$. See appendix A for the
proof.

###  III-B Low-Mobility User Detection

After successfully decoding and subtracting the HM user’s signal, each LM user
$u,1\leq u\leq U_{1}$ observes the following received signal in the TF domain,

$\displaystyle y_{u}[n,m]$
$\displaystyle=\sum_{a=1}^{A}p_{a,1}{H}_{a,u}[n,m]\sum_{i=1}^{U_{1}}s_{i,1}[n,m]$
$\displaystyle+\sum_{a=1}^{A}{H}_{a,u}[n,m]\sum_{q=2}^{Q}p_{a,q}s_{q}[n,m]+w_{u}[n,m],$
(28)

where $H_{a,u}[n,m]$ is the channel gain between antenna $a$ and user $u$ and
$w_{u}[n,m]$ is the noise in the TF domain. The samples are conducted at
instant $nT$ and channel $m\Delta f$. The TF domain channel gain can be
expressed as [9],

$H_{a,u}[n,m]=\int_{\tau}\int_{\nu}h_{a,u}(\tau,\nu)e^{j2\pi\nu
nT}e^{-j2\pi(\nu+m\Delta f)\tau}d\tau d\nu,$ (29)

where $h_{a,u}(\tau,\nu)$ is the channel impulse response between antenna $a$
and user $u$. Considering the absence of Doppler shifts in low-mobility
scenarios, the channel gain in (29) can be simplified to a time-independent
expression as follows,

$H_{i,i^{\prime},u}[n,m]=\sum_{p=1}^{P_{u}}h_{p,u}e^{j2\pi m\Delta
f\tau_{p,u}}v_{ii^{\prime}}\left(\theta_{p,u},\phi_{p,u}\right).$ (30)

Since $\tau_{p,u}=l_{p,u}/(M\Delta f)$ where $l_{p,u}$ is the delay tap pf
path $p$, the channel gain can be written as,

$H_{i,i^{\prime},u}[m]=\sum_{p=1}^{P_{u}}h_{p,u}e^{j2\pi\frac{l_{{p,u}}m}{M}}v_{ii^{\prime}}\left(\theta_{p,u},\phi_{p,u}\right).$
(31)

Considering the transmitted signal expression in (2) and thanks to the
orthogonality between LM users within the same cluster, user $u$ disregards
the signals received on subchannels $m\neq u-1$ and the expression in (III-B)
can be re-written as,

$\displaystyle y_{u}[n,u-1]=\sum_{a=1}^{A}p_{a,1}{H}_{a,u}[u-1]s_{u,1}[n,u-1]$
$\displaystyle+\sum_{q=2}^{Q}\sum_{a=1}^{A}{H}_{a,u}[u-1]p_{a,q}s_{q}[n,u-1]+w_{u}[n,u-1].$
(32)

Since the bandwidth of each subcarrier is considered smaller than the
coherence bandwidth of the channel, one-tap equalization can be used to detect
the LM user’s signal [3]. The equalized signal of LM user $u$ at time instant
$n$ is expressed as,

$\displaystyle\frac{y_{u}[n,u-1]}{H_{u}[u-1]}=s_{u}[n,u-1]$
$\displaystyle+\sum_{q=2}^{Q}\frac{\sum_{a=1}^{A}p_{a,q}H_{a,u}[u-1]}{H_{u}[u-1]}s_{q}[n,u-1]+\frac{w_{u}[n,u-1]}{H_{u}[u-1]}.$
(33)

Figure 2: Outage probability of HM and LM users versus the number of LM users
for different modulation schemes.

Since LM users are time-invariant, the SINR expression is the same for all
time instants. Therefore, the SINR for detecting the user $u$ is given by,

$\gamma_{u}=\frac{\rho_{u}\left|{H_{u}[u-1]}\right|^{2}}{\sum_{q=2}^{Q}\rho_{c,q}\left|{\sum_{a=1}^{A}p_{a,q}H_{a,u}[u-1]}\right|^{2}+1}.$
(34)

##  IV Simulation Results

To evaluate the performance of the proposed system, we consider the following
simulation parameters: $A=64,N=16,M=16,\Delta
f=15\;\mathrm{kHz},f_{\mathrm{c}}=60\;GHz$. As in previous studies [2],
Doppler taps are generated according to Jake’s model, considering that the
maximum vehicle speed is 200 km/h. The maximum channel delay tap is
$l_{\max}=4$. The number of paths is 4 for all users. The elevation and
azimuth AoAs of paths are uniformly distributed with mean $\pi/4$ and $0$,
respectively, and with variance $\pi/10$ and $\pi$, respectively. Unless
otherwise specified, we consider the following parameters: $Q=3$,
$\rho_{\mathrm{T}}=30\;\mathrm{dB}$. The SNR of the HM user of cluster 1 is
$\rho_{0}=\alpha\rho_{c,1}$, where $\alpha=3/4$. Therefore, the SNR of each LM
user of cluster 1 is $\rho_{u}={(1-\alpha)\rho_{c,1}}/{U_{1}},1\leq u\leq
U_{1}$. The number of LM users in each cluster is selected uniformly at random
in $[\,1,M]\,$. The LM users’ positions are uniformly placed in a circle
around the HM user of each cluster with a radius
$R_{\mathrm{Q}}=10\;\mathrm{m}$. In the following simulation examples, the
outage probability is the probability that the rate falls below the threshold
$R_{\mathrm{th}}=0.5\;\mathrm{b/s/Hz}$.

Figure 3: Outage probability of HM and LM users for different power allocation
factors. Figure 4: Benchmark comparison of the proposed OTFS-NOMA system with
OMA access.

In the first simulation example, we investigate the effect of the number of LM
users within each cluster on the system performance. To this end, we plot in
Fig. 2 the outage probability of the HM and LM users versus the number of LM
users for various transmit SNRs. First, we observe that the outage probability
of the HM user is quasi-constant as a function of the number of LM users
within the cluster. This result is valid for both values of the transmit SNR.
This can be justified by the power allocation scheme that is adopted by the
system. The BS allocates a fixed fraction of the power to the HM user and
distributes the remaining power equally among the LM users. Moreover, since
each LM user has access to a single subchannel, some subchannels are unused
when $U_{q}<M$. Hence, we note that the behavior of the HM outage probability
remains constant, implying that the effect of the distribution of interference
coming from LM users between subchannels does not have a significant impact on
the performance of the HM user. In other words, given that the HM OTFS frame
spans the whole resource block, the HM user’s performance remains roughly the
same whether the interference is concentrated within a narrow frequency band
or spread over a larger band. Second, we note from Fig. 2 that increasing
transmit SNR from 20 dB to 30 dB has a significant performance improvement in
terms of quality-of-service (QoS) delivered to the users. Although the
decrease of the outage probability is more significant at the HM user’s side,
LM users also benefit from the resource boost. Since the transmit power
allocated for LM users decreases with their number, increasing the transmit
SNR or tuning the allocated power to LM users helps compensate for the
performance degradation and reach acceptable QoS.

In the second simulation example, we investigate the effect of the power
allocation factor on the system performance. In Fig. 3, the outage probability
is shown versus the transmit SNR at the BS for different allocation factors.
We note that the power distribution between the HM user and LM users has a
significant performance impact, especially for low SNRs. In high SNR regimes,
this effect becomes less significant. Moreover, when $\alpha=4/5$, the
performance gap between the HM and LM users is much more significant than the
performance gap in the case of $\alpha=1/2$. Subsequently, increasing the
transmit power of the HM user will certainly decrease its outage probability,
but at the same time, it can deteriorate the performance of LM users in the
process. As a result, it is important to strategically tune the power
allocation factor in order to keep acceptable QoS for all users.

In the last simulation example, we compare the proposed NOMA system
configuration with an orthogonal multiple access (OMA) scheme. Herein, we
adopt the time division multiple access between HM and LM users. To this end,
in Fig. 4, we present the outage probability of the HM user’s signal detection
at the HM user’s side versus the total DL transmit SNR for different numbers
of clusters. We can observe that the proposed OTFS-NOMA scheme outperforms its
orthogonal counterpart and provides more reliable communication links for all
the examined values of the transmit SNR. Also, the performance gap between the
proposed NOMA setup and OMA benchmark increases with the number of clusters
within the cell. This behavior suggests that the reliability of the proposed
system becomes more significant and impactful in low SNR regimes.

##  V Conclusion

In this paper, we proposed an OTFS-NOMA system where spatial diversity is
adopted to convey data signals to multiple clusters of users simultaneously.
Using NOMA paradigm, the BS superimposes an OTFS signal frame transmitted to
the HM user and TF-modulated signals transmitted to LM users. The LM users
access disjoint subchannels. Subsequently, their signals are, by design,
orthogonal to each other. The results demonstrated that the power distribution
between users and the number of LM users have a remarkable impact on the
outage probability of LM users. However, the performance of the HM user
remains unchanging when the number of LM users increases. Moreover, the
proposed framework outperforms the OMA setup, especially with higher numbers
of clusters.

##  Appendix A SNR Expression

The SNR of $k^{\mathrm{th}}$ symbol detection is expressed as,

$\gamma_{u}^{0}=\frac{\mathbb{E}\left(|\mathbf{s}_{0,1}[k]|^{2}\right)}{\sum_{i=1}^{U_{1}}\\!\mathbb{E}\left(\left|\mathbf{s}_{i,1}[k]\right|^{2}\right)\\!+\\!\sum_{q=2}^{Q}\\!\mathbb{E}\left(\mathbf{\Omega}_{q,u}\mathbf{s}_{q}\mathbf{s}_{q}^{\mathrm{H}}\mathbf{\Omega}_{q,u}^{\mathrm{H}}\right)_{k,k}\\!+\\!\mathbb{E}\left(\mathbf{z}_{u}\mathbf{z}_{u}^{\mathrm{H}}\right)_{k,k}}$
(35)

Taking into account that
$\mathbf{w}_{u,q}\\!\sim\\!\mathcal{CN}(\mathbf{0},{\sigma}_{\mathrm{w}}^{2}\mathbf{I}_{N\\!M})$
and
$\mathbf{G}_{u}\mathbf{G}_{u}^{\mathrm{H}}\\!=\\!\left(\mathbf{H}_{u}^{\mathrm{H}}\mathbf{H}_{u}\right)^{-1}$,
the noise covariance matrix is
$\mathbb{E}\left(\mathbf{z}_{u}\mathbf{z}_{u}^{\mathrm{H}}\right)\\!=\\!\mathbb{E}\left(\mathbf{G}_{u}\mathbf{w}_{u}\mathbf{w}_{u}^{\mathrm{H}}\mathbf{G}_{u}^{\mathrm{H}}\right)\\!=\\!\sigma_{\mathrm{w}}^{2}\left(\mathbf{H}_{u}^{\mathrm{H}}\mathbf{H}_{u}\right)^{{-1}}.$
We assume that the transmit power for cluster $q$ is
$\sum_{u=0}^{U_{q}}P_{u,q}=P_{q},\forall q\in\\{1,\ldots,Q\\}$. The vector
$\mathbf{s}_{q}\sim\mathcal{CN}(\mathbf{0},\sum_{0\leq u\leq
U_{q}}P_{u,q}\mathbf{I}_{N\\!M})$. Thus, the covariance of the inter-cluster
interference term can be re-written as,

$\displaystyle\sum_{q=2}^{Q}\mathbb{E}\left(\mathbf{\Omega}_{q,u}\mathbf{s}_{q}\mathbf{s}_{q}^{\mathrm{H}}\mathbf{\Omega}_{q,u}^{\mathrm{H}}\right)$
$\displaystyle=\sum_{q=2}^{Q}P_{q}\mathbf{\Omega}_{q,u}\mathbf{\Omega}_{q,u}^{\mathrm{H}}.$
(36)

If we denote by $\eta_{u,i}$ the $i^{\mathrm{th}}$ eigenvalue of
$\mathbf{H}_{u}$, then, according to (15), it can be expressed as
$\eta_{u,i}=\sum_{a=1}^{A}p_{a,1}\lambda_{a,u,i}$. According to [11], the term
$\left(\left(\mathbf{H}_{u}^{\mathrm{H}}\mathbf{H}_{u}\right)^{-1}\right)_{k,k}$
can be simplified as follows,

$\left(\left(\mathbf{H}_{u}^{\mathrm{H}}\mathbf{H}_{u}\right)^{-1}\right)_{k,k}=\frac{1}{NM}\sum_{i=1}^{NM}\frac{1}{\left|\eta_{u,i}\right|^{2}}.$
(37)

Similarly, if we denote by $\omega_{q,u,i}$ the $i^{\mathrm{th}}$ eigenvalue
of $\mathbf{\Omega}_{q,u}$, then, according to (25), it can be expressed as,

$\omega_{q,u,i}=\frac{\sum_{a=1}^{A}p_{a,1}^{*}\lambda_{a,u,i}^{*}}{\left|\sum_{a=1}^{A}p_{a,1}\lambda_{a,u,i}\right|^{2}}\sum_{a=1}^{A}p_{a,q}\lambda_{a,u,i}.$
(38)

The $(k,k)^{\mathrm{th}}$ element of
$\mathbf{\Omega}_{q,u}\mathbf{\Omega}_{q,u}^{\mathrm{H}}$ can be simplified as
follows,

$\left(\mathbf{\Omega}_{q,u}\mathbf{\Omega}_{q,u}^{\mathrm{H}}\right)_{k,k}=\frac{1}{NM}\sum_{i=1}^{NM}{\left|\omega_{q,u,i}\right|^{2}}.$
(39)

Normalizing the SNR by $\sigma_{\mathrm{w}}^{2}$, we get the expression in
(26).

## References

  * [1] S. R. Islam, M. Zeng, O. A. Dobre, and K.-S. Kwak, “Resource allocation for downlink NOMA systems: Key techniques and open issues,” _IEEE Wirel. Commun._ , vol. 25, no. 2, pp. 40–47, 2018.
  * [2] P. Raviteja, K. T. Phan, Y. Hong, and E. Viterbo, “Interference cancellation and iterative detection for orthogonal time frequency space modulation,” _IEEE Trans. Wirel. Commun.s_ , vol. 17, no. 10, pp. 6501–6515, 2018.
  * [3] Z. Ding, R. Schober, P. Fan, and H. V. Poor, “OTFS-NOMA: An efficient approach for exploiting heterogenous user mobility profiles,” _IEEE Trans. Commun._ , vol. 67, no. 11, pp. 7950–7965, 2019.
  * [4] Z. Ding, “Robust beamforming design for OTFS-NOMA,” _IEEE Open. J. Commun. Soc._ , vol. 1, pp. 33–40, 2019.
  * [5] X. Zhou, K. Ying, Z. Gao, Y. Wu, Z. Xiao, S. Chatzinotas, J. Yuan, and B. Ottersten, “Active terminal identification, channel estimation, and signal detection for grant-free NOMA-OTFS in LEO satellite Internet-of-Things,” _IEEE Trans. Wirel. Commun._ , vol. 22, no. 4, pp. 2847–2866, 2022.
  * [6] X. Zhou and Z. Gao, “Joint active user detection and channel estimation for grant-free NOMA-OTFS in LEO constellation Internet-of-Things,” in _IEEE/CIC Inter. Conf. Commun. China. (ICCC)_ , 2021, pp. 735–740.
  * [7] Y. Ge, Q. Deng, P. Ching, and Z. Ding, “OTFS signaling for uplink NOMA of heterogeneous mobility users,” _IEEE Trans. Commun._ , vol. 69, no. 5, pp. 3147–3161, 2021.
  * [8] Z. D. Zaharis, I. P. Gravas, P. I. Lazaridis, T. V. Yioultsis, and T. D. Xenos, “Improved beamforming in 3D space applied to realistic planar antenna arrays by using the embedded element patterns,” _IEEE Trans. Veh. Technol._ , vol. 71, no. 6, pp. 6145–6157, 2022.
  * [9] P. Singh, A. Gupta, H. B. Mishra, and R. Budhiraja, “Low-complexity ZF/MMSE MIMO-OTFS receivers for high-speed vehicular communication,” _IEEE Open. J. Commun. Soc._ , vol. 3, pp. 209–227, 2022.
  * [10] B. C. Pandey, S. K. Mohammed, P. Raviteja, Y. Hong, and E. Viterbo, “Low complexity precoding and detection in multi-user massive MIMO OTFS downlink,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 5, pp. 4389–4405, 2021.
  * [11] P. Singh, K. Yadav, H. B. Mishra, and R. Budhiraja, “BER analysis for OTFS zero forcing receiver,” _IEEE Trans. Commun._ , vol. 70, no. 4, pp. 2281–2297, 2022.