# Electromagnetic radiation at extreme angular velocity

Matteo Buzzegoli  Kirill Tuchin  Department of Physics and Astronomy, Iowa
State University, Ames, Iowa, 50011, USA

###### Abstract

We consider a system rotating at extremely high angular velocity, so that its
matter is found mostly at the light-cylinder. We posit that it can be
described by quantum fields confined to the two-dimensional cylindrical
surface rotating about its symmetry axis. We apply this model to study the
electromagnetic radiation. In particular, we compute the photon spectrum
emitted by the quark-gluon plasma.

##  I Introduction

The interest to the rapidly rotating systems has been recently rekindled
thanks to the experimental observation of highly vortical quark-gluon plasma
produced in the relativistic heavy-ion collisions [1, 2, 3, 4, 5, 6]. Previous
studies discussed the thermodynamics and the hydrodynamics of rotating systems
based on the statistical approach [7, 8, 9, 10, 11, 12, 13], developed a
quantum kinetic theory with spin degrees of freedom [14, 15, 16, 17, 18, 19,
20, 21, 22], and made predictions for the spin polarization measured in heavy-
ion collisions [23, 24, 25, 26, 27, 28, 29, 30, 31, 32], see for instance the
reviews [33, 34, 35].

In [36, 37] we initiated a study of the electromagnetic radiation emitted by
rapidly rotating systems. The advantage of the electromagnetic radiation is
that it is only weakly affected by the plasma evolution. The idea is to
observe the impact of rotation on the quantum fields.

The study in [36, 37] discussed “relatively slowly” rotating systems in
magnetic field. Namely, we assumed that the magnetic length $1/\sqrt{eB}$ is
much shorter than the inverse angular velocity $\Omega$***We are using natural
units where $c=\hbar=k_{B}=1$.. Such rotation is slow. On the other hand, the
absolute value of the angular velocity satisfying this condition can be
enormous, hence the qualifying adverb “relatively”. Generally, we can say that
a system is relatively slowly rotating if its transverse size $a$ is much
smaller than $1/\Omega$.

Model simulations show that the vorticity of the quark-gluon plasma can be as
high as its inverse transverse size $a$. This upsets the slow rotation
assumption. A system rotating with the angular velocity $\Omega$ is causally
connected only within the lightcone cylinder of radius $R=1/\Omega$. When
$R<a$ only a part of the rotating plasma is causally connected. This is a
genuine fast rotation. Setting the proper boundary conditions on the quantum
fields at the causal boundary becomes an essential procedure. In [38], in the
spirit of the MIT bag model, we required that the radial current vanishes on
the boundary. However, there may be other possible boundary conditions.

In the regime $1/\Omega<\ell$, where $\ell$ is the mean-free-path, the
rotation is so fast that it overwhelms all inter-particle forces and pushes
the matter towards the light-cylinder wall. Such a medium will break down to a
set of rotating cylindrical regions of radii $R\ll a$. Within each cylindrical
region the matter will be concentrated mostly at the boundary at $R$ due to
the centrifugal force. It seems reasonable therefore, that the dynamics of
such extremely rapidly rotating system can be described by the quantum fields
confined to the cylindrical surface of radius $R$. We will not be concerned
with the statistical properties of the matter within the cylindrical region,
since in view of $\ell>R$, it is simply an ideal rotating gas. Rather we are
interested in the electromagnetic radiation it emits. Since the precise nature
of the particles that make up the rotating cylinder is not very important —
only the fact that they are found at the boundary is — we employ the scalar
QED for calculations. It is reasonable to expect that the qualitative features
of our results should be fairly model-independent.

The three rotation regimes of a system characterized by the radial size $a$,
the light-cylinder radius $R=1/\Omega$, the mean-free-path $\ell$ are

  * •

_Slow rotation_ $\ell\ll a\ll R$. This approximation is used in [36, 37].

  * •

_Fast rotation_ $\ell\ll R\sim a$ [39, 38].

  * •

_Extremely fast rotation_ $R\ll\ell\ll a$. This is the scenario we consider in
the paper.

In summary, we consider model in which the charged scalar particles can freely
move on a thin cylindrical sheet of radius $R=1/\Omega$ rotating with angular
velocity $\Omega$. In the following sections we will compute the
electromagnetic radiation by a single particle and by a system of particles in
thermal equilibrium. We believe that this model describes the universal
properties of extremely rapidly rotating systems.

##  II Radiation by scalar particle on cylindrical sheet

The wave function of a scalar particle of mass $M$ embedded into a cylindrical
surface rotating with the angular velocity $\Omega$ about its symmetry axis
$z$ is

$\displaystyle\psi(t,\phi,z)=\frac{1}{\sqrt{2\pi
L}}\frac{1}{\sqrt{2E}}e^{-iEt+ip_{z}z+im\phi}\,,$ (1)

where energy spectrum is

$\displaystyle E=\sqrt{p_{z}^{2}+\frac{m^{2}}{R^{2}}+M^{2}}+m\Omega\,.$ (2)

The magnetic quantum number $m$ is an integer, while the longitudinal momentum
$p_{z}$ is continuous assuming that the cylinder height $L$ is very large. We
note that $E>0$ for any $m$.

The velocity of the quasi-classical motion along the $z$-direction is
$v=p_{z}/E$. Clearly, only the states with $|v|\leq 1$ are causally connected.
Inspection of (2) reveals that when $m\geq 0$, then $E>|p_{z}|$ for any value
of $p_{z}$. In contrast, when $m<0$ this condition is satisfied only when

$\displaystyle|p_{z}|\leq\frac{M^{2}}{2|m|\Omega}\,,\quad\text{if}\quad
m<0\,.$ (3)

Fig. 1 shows an example of the dispersion relation (2) with $m<0$. Vertical
lines indicated the allowed range of $p_{z}$’s.

Figure 1: (Color online) $m\Omega=-2$, $M=1$. Green line: $v(p_{z})$, red
line: $E(p_{z})$. The allowed values: $-\frac{1}{4}\leq p_{z}\leq\frac{1}{4}$
corresponding to the box in the center.

We are interested to compute the electromagnetic radiation by this particle.
The $S_{fi}$-matrix element reads

$\displaystyle S_{fi}=-ie\int dt\int d\phi\int
dz\,{\bm{j}}_{fi}(t,\phi,z)\cdot{\bm{A}}^{*}(t,\phi,z,R)\,,$ (4)

where the photon wave function in the Coulomb gauge is

$\displaystyle{\bm{A}}(t,\phi,z,R)=\frac{1}{\sqrt{2\omega
V}}{\bm{\epsilon}}_{\lambda}e^{i{\bm{k}}\cdot{\bm{r}}-i\omega
t}\,,\quad{\bm{\epsilon}}_{\lambda}\cdot{\bm{k}}=0\,,$ (5)

and $V$ is the system volume. The transition current is

$\displaystyle{\bm{j}}_{fi}=i(\psi{\bm{\nabla}}\psi^{\prime*}-\psi^{\prime*}{\bm{\nabla}}\psi)\,.$
(6)

In our notation: $\psi_{i}=\psi$, $\psi_{f}=\psi^{\prime}$.

Substituting (1) and (5) into (4) we arrive at

$\displaystyle S_{fi}=$ $\displaystyle\frac{-ie(2\pi)^{2}}{2\pi L\sqrt{2\omega
V}2\sqrt{EE^{\prime}}}\delta(E-E^{\prime}-\omega)\delta(p_{z}-p_{z}^{\prime}-k_{z})$
$\displaystyle\times{\bm{\epsilon}}^{*}_{\lambda}\cdot\left[\frac{1}{R}{\bm{e}}_{\phi}(m+m^{\prime})+{\bm{e}}_{z}(p_{z}+p_{z}^{\prime})\right]2\pi
iJ_{m-m^{\prime}}(k_{\bot}R)\,.$ (7)

Two convenient photon polarizations (following [40]):

$\displaystyle{\bm{\epsilon}}_{1}=-{\bm{e}}_{\phi}\,,\quad{\bm{\epsilon}}_{2}=-\sin\theta{\bm{e}}_{z}+\cos\theta{\bm{e}}_{\bot}\,.$
(8)

where $\theta$ is the polar angle defined with respect to the $z$-axis, e.g.
$k_{z}=\omega\cos\theta$. The photon transverse momentum is then
$k_{\bot}=\omega\sin\theta$.

The photon emission rate can be computed as

$\displaystyle\dot{w}_{fi}$
$\displaystyle=\sum_{\lambda}\sum_{m^{\prime}}\frac{|S_{fi}|^{2}}{T}\frac{dk_{z}L}{2\pi}\frac{dk_{\bot}k_{\bot}\pi
R^{2}}{2\pi}\frac{dp_{z}L}{2\pi}$ (9)

while the radiation intensity is given by $W=\dot{w}_{fi}\omega$:

$\displaystyle W$
$\displaystyle=\sum_{m^{\prime}=-\infty}^{m-1}\frac{e^{2}}{16\pi
EE^{\prime}}\delta(E-E^{\prime}-\omega)\left[\frac{(m+m^{\prime})^{2}}{R^{2}}+\sin^{2}\theta(p_{z}+p_{z}^{\prime})^{2}\right]J_{m-m^{\prime}}^{2}(k_{\bot}R)dk_{z}dk_{\bot}k_{\bot}\,.$
(10)

The sum over the negative $m^{\prime}$ in (10) is constrained by (3). One can
expresses $dk_{z}dk_{\bot}k_{\bot}=\omega^{2}d\omega
d\sin\theta=\omega^{2}d\omega do/2\pi$, where $do$ is the element of the solid
angle in the direction of the emitted photon. In the non-relativistic limit,
the leading term in the multipole expansion of the intensity is the magnetic
dipole one because the electric dipole moment vanishes while the magnetic
moment ${\bm{\mu}}=\frac{1}{2}e{\bm{r}}\times{\bm{v}}$ is finite.

The delta-function in (10) can be re-written as

$\displaystyle\delta(E-E^{\prime}-\omega)=\frac{\delta(\omega-\omega_{0})(E^{\prime}-m^{\prime}\Omega)}{E-m^{\prime}\Omega-\omega\sin^{2}\theta-
p_{z}\cos\theta}\,.$ (11)

where the characteristic frequency is

$\displaystyle\omega_{0}$
$\displaystyle=\frac{1}{\sin^{2}\theta}\bigg{\\{}E-m^{\prime}\Omega-
p_{z}\cos\theta$ $\displaystyle-\sqrt{(E-m^{\prime}\Omega-
p_{z}\cos\theta)^{2}-\sin^{2}\theta\left[(E-m^{\prime}\Omega)^{2}-p_{z}^{2}-m^{\prime
2}/R^{2}-M^{2}\right]}\bigg{\\}}\,.$ (12)

Taking the integral over $\omega$ one is left with the angular spectrum
$dW/do$. Alternatively, we can cast the delta-function in the form

$\displaystyle\delta(E-E^{\prime}-\omega)=\sum_{\pm}\frac{\delta(\cos\theta-\cos\theta_{\pm})(E^{\prime}-m^{\prime}\Omega)}{\omega|\omega\cos\theta-
p_{z}|}\,,$ (13)

where

$\displaystyle\cos\theta_{\pm}=\frac{1}{\omega}\left\\{p_{z}\pm\sqrt{(E-m^{\prime}\Omega-\omega)^{2}-m^{\prime
2}/R^{2}-M^{2}}\right\\}\,.$ (14)

We note that in order that $|\cos\theta_{\pm}|\leq 1$, photon energy $\omega$
must not be too small. For the states with $m^{\prime}<0$ one has to take into
account the requirement (3).

The radiation intensity is shown in Fig. 2 as a function of $m$ and $\Omega$.

|
---|---
Figure 2: Left panel: $W$ vs $m$ at $p_{z}=0$, right panel: $W/\Omega^{2}$ vs
$\Omega$. In both panels $M=1$.

##  III Radiation intensity by a single particle at rest

To investigate the radiation intensity analytically, consider a reference
frame in which the incident particle is at rest in the axial direction:
$p_{z}=0$, $p_{z}^{\prime}=-\omega\cos\theta$. In particular, (3) becomes

$\displaystyle\omega_{0}|\cos\theta|\leq\frac{M^{2}}{2|m^{\prime}|\Omega}\,.$
(15)

In addition, we will focus on the limit $\Omega\gg M$. In this case the
spectrum (2) reads

$\displaystyle E\approx\left\\{\begin{array}[]{lc}2m\Omega\,,&m>0\,;\\\
M\,,&m=0\,;\\\ \frac{M^{2}}{2|m|\Omega}\,,&m<0\,.\end{array}\right.\qquad
E^{\prime}\approx\left\\{\begin{array}[]{lc}\sqrt{p_{z}^{\prime 2}+m^{\prime
2}\Omega^{2}}+m^{\prime}\Omega\,,&m^{\prime}>0\,;\\\ \sqrt{p_{z}^{\prime
2}+M^{2}}\,,&m^{\prime}=0\,;\\\ \frac{p_{z}^{\prime
2}+M^{2}}{2|m^{\prime}|\Omega}\,,&m^{\prime}<0\,.\end{array}\right.$ (22)

In order that $E>E^{\prime}$, the range of $m^{\prime}$ must be
$m^{\prime}\leq m-1$.

We distinguish three cases: (A) $m>0$, (B) $m=0$ and (C) $m<0$.

###  III.1 $m>0$

Fig. 3 shows the dependence of a single term in the sum (10) on $m^{\prime}$
for $m>0$. The mode with $m^{\prime}=0$ is enhanced over those with
$m^{\prime}>0$ which in turn are enhanced as compared to the modes
$m^{\prime}<0$.

Figure 3: (Color online) The radiation intensity $W_{m^{\prime}}$ of a single
term in (10) for $m=5$ vs $m^{\prime}$. The total intensity $W$ is the sum of
all $W_{m^{\prime}}$. The blue and orange symbols overlap at $m^{\prime}=-1$.

To understand the dynamics in each case and find an appropriate approximation
we split the summation over $m^{\prime}$ into three parts: (A1)
$m^{\prime}\leq-1$, (A2) $m^{\prime}\geq 1$ and (A3) $m^{\prime}=0$.

####  III.1.1 $m^{\prime}<0$

Using (22),(11),(II) in (10) and keeping the leading terms in $M/\Omega$ we
obtain

$\displaystyle
W^{A1}=\frac{e^{2}}{2\pi}\sum_{m^{\prime}=-\infty}^{-1}\int_{0}^{\pi}\frac{1}{16\pi
EE^{\prime}}\frac{\delta(2m\Omega-\omega)}{\frac{\omega\cos^{2}\theta}{E^{\prime}-m^{\prime}\Omega}+1}\left[\frac{(m+m^{\prime})^{2}}{R^{2}}+\sin^{2}\theta\cos^{2}\theta\omega^{2}\right]$
$\displaystyle\times
J_{m-m^{\prime}}^{2}(k_{\bot}R)\omega^{2}\eta\left(M^{2}+2m^{\prime}\Omega\omega|\cos\theta|\right)d\omega\sin\theta
d\theta$ (23)
$\displaystyle=\frac{e^{2}}{2\pi}\sum_{m^{\prime}=-\infty}^{-1}2\int_{0}^{1}\frac{E^{\prime}-m^{\prime}\Omega}{16\pi
EE^{\prime}(2m\Omega\cos^{2}\theta+E^{\prime}-m^{\prime}\Omega)}\left[\frac{(m+m^{\prime})^{2}}{R^{2}}+(1-x^{2})x^{2}(2m\Omega)^{2}\right]$
$\displaystyle\times
J_{m-m^{\prime}}^{2}\left(2m\sqrt{1-x^{2}}\right)(2m\Omega)^{2}\eta\left(M^{2}+2m^{\prime}\Omega\omega
x\right)dx\,,$ (24)

where $\eta$ is the step function accounting for (15), $x=\cos\theta$ and we
took advantage of the fact that the integrals over $0\leq x\leq 1$ and $-1\leq
x\leq 0$ are equal. In view of (22) and (15), the angular integration is
confined to the small region

$\displaystyle x<\frac{M^{2}}{4\Omega^{2}m|m^{\prime}|}\ll 1\,.$ (25)

Bearing this in mind and using (22) we can then write
$\omega_{0}\cos^{2}\theta+E^{\prime}\approx\frac{M^{2}}{2|m^{\prime}|\Omega}\ll\Omega$.
After expanding the remaining terms in $x$ we obtain

$\displaystyle
W^{A1}=\frac{e^{2}\Omega^{2}}{8\pi}\sum_{j=1}^{\infty}(m-j)^{2}J^{2}_{m+j}(2m)=\frac{e^{2}\Omega^{2}}{8\pi}S_{1}(m)\,,$
(26)

where $j=-m^{\prime}$ and we defined

$\displaystyle S_{1}(m)=\sum_{j=1}^{\infty}(m-j)^{2}J^{2}_{m+j}(2m)\,.$ (27)

This function is shown in Fig. 4.

|
---|---
Figure 4: Left panel: $S_{1}(m)/m^{2}$, right panel: $S_{2}(m)/m^{2}$
appearing in (27) and (30) respectively.

####  III.1.2 $m^{\prime}>0$

Now consider the final states with $m^{\prime}>0$. In this case the constraint
(15) does not apply. It implies that we can safely neglect $M$ in (II),
assuming that the integral over $\theta$ is finite at $\theta=\pi/2$ (which is
indeed the case as will be seen shortly). Thus,

$\displaystyle\omega_{0}\approx\frac{\Omega}{\sin^{2}\theta}\left\\{2m-m^{\prime}-\sqrt{(2m-m^{\prime})^{2}-4m(m-m^{\prime})\sin^{2}\theta}\right\\}\equiv\Omega
z\,,$ (28)

and $E^{\prime}=2\Omega m-\omega_{0}$. A simple calculation yields

$\displaystyle W^{A2}=\frac{e^{2}\Omega^{2}}{32\pi}S_{2}\left(m\right)\,,$
(29)

where we defined the function

$\displaystyle S_{2}(m)=$
$\displaystyle\sum_{m^{\prime}=1}^{m-1}\int_{0}^{\pi}\frac{(2m-m^{\prime}-z)\left[(m+m^{\prime})^{2}+z^{2}\sin^{2}\theta\cos^{2}\theta\right]}{m(2m-z)(2m-m^{\prime}-z\sin^{2}\theta)}J^{2}_{m-m^{\prime}}(z\sin\theta)z^{2}\sin\theta
d\theta\,$ (30)

which depends only on $m$. The numerical calculation of $S_{2}$ is shown in
Fig. 4.

####  III.1.3 $m^{\prime}=0$

In this case the photon frequencies are

$\displaystyle\omega_{0}=\frac{2m\Omega}{\sin^{2}\theta}\left\\{1-\sqrt{\cos^{2}\theta+\sin^{2}\theta\frac{M^{2}}{4m^{2}\Omega^{2}}}\right\\}$
(31)

and the intensity is

$\displaystyle
W^{A3}=\frac{e^{2}}{16\pi(2m\Omega)}\int_{0}^{\pi}\frac{m^{2}\Omega^{2}+\omega_{0}^{2}\sin^{2}\theta\cos^{2}\theta}{\omega_{0}\cos^{2}\theta+\sqrt{\omega_{0}^{2}\cos^{2}\theta+M^{2}}}J_{m}^{2}(\omega_{0}R\sin\theta)\omega_{0}^{2}\sin\theta
d\theta\,.$ (32)

The integrand peaks at small $x=\cos\theta$ which gives rise to the
logarithmically enhanced contribution. In this leading logarithmic
approximation we write

$\displaystyle
W^{A3}\approx\frac{e^{2}}{16\pi}2\int_{0}^{1}\frac{m^{2}\Omega^{2}}{x^{2}+\sqrt{x^{2}+\frac{M^{2}}{4m^{2}\Omega^{2}}}}J_{m}^{2}(2m)dx\approx\frac{e^{2}\Omega^{2}}{8\pi}m^{2}J_{m}^{2}(2m)\log\frac{4m\Omega}{M}\,.$
(33)

###  III.2 $m=0$

In this case all finial states have $m^{\prime}<0$. Using (10)-(14) and
replacing $-m^{\prime}=j$ we get

$\displaystyle
W^{B}=\sum_{j=1}^{\infty}\frac{e^{2}(E^{\prime}-m^{\prime}\Omega)}{16\pi
EE^{\prime}(E^{\prime}-m^{\prime}\Omega+\omega_{0}\cos^{2}\theta)}\left[\frac{j^{2}}{R^{2}}+\omega_{0}^{2}\sin^{2}\theta\cos^{2}\theta\right]\eta\left(M^{2}-2j\Omega\omega|\cos\theta|\right)$
$\displaystyle\times J^{2}_{j}(\omega_{0}R\sin\theta)\omega^{2}_{0}\sin\theta
d\theta\,,$ (34)

where, considering that $x=\cos\theta\ll 1$ (see (25)) we approximate:

$\displaystyle\omega_{0}\approx
M\left(1-\frac{M}{2j\Omega}(1+x^{2})\right)\,.$ (35)

It follows that $E^{\prime}+\omega_{0}\cos^{2}\theta+j\Omega\approx j\Omega$.
We are now left with a trivial integral over $x$ which yields:

$\displaystyle
W^{B}=\frac{e^{2}\Omega^{2}}{8\pi}\sum_{j=1}^{\infty}j^{2}J_{j}^{2}(MR)$ (36)

Expanding the Bessel function and keeping only the leading term with $j=1$ we
finally obtain

$\displaystyle W^{B}=\frac{e^{2}M^{2}}{32\pi}\,.$ (37)

###  III.3 $m<0$

This is the most unusual case because the fermion transitions occur between
the levels with negative $m$ and $m^{\prime}$ which correspond to
$0<E^{\prime}<E<M$. First deal with the delta-function:

$\displaystyle\delta\left(-\frac{M^{2}}{2m\Omega}+\frac{\omega^{2}\cos^{2}\theta+M^{2}}{2m^{\prime}\Omega}-\omega\right)=\frac{1}{\frac{\omega\cos^{2}\theta}{m^{\prime}\Omega}-1}\delta\left(\omega-\frac{M^{2}(|m^{\prime}|-|m|)}{2|m^{\prime}||m|\Omega}\right)\,.$
(38)

were we used (22). We get for the intensity using (10):

$\displaystyle
W^{C}=\frac{e^{2}}{16\pi}\sum_{m^{\prime}=-\infty}^{m-1}\frac{2|m|\Omega}{M^{2}}\frac{2|m^{\prime}|\Omega}{\omega_{0}^{2}\cos^{2}\theta+M^{2}}\frac{|m^{\prime}|\Omega}{\omega_{0}\cos^{2}\theta+|m^{\prime}|\Omega}\left[\frac{(m+m^{\prime})^{2}}{R^{2}}+\omega^{2}_{0}\sin^{2}\theta\cos^{2}\theta\right]$
$\displaystyle\times\eta\left(M^{2}+2m^{\prime}\Omega\omega|\cos\theta|\right)J^{2}_{m-m^{\prime}}(\omega_{0}R\sin\theta)\omega_{0}^{2}\sin\theta
d\theta\,.$ (39)

In view of (38), the step function implies that $|\cos\theta|\leq|m|/\nu$,
where $\nu=|m^{\prime}|-|m|$ is a positive integer. Since $\omega_{0}\sim
M^{2}/\Omega$, we can neglect it in the denominators and in the square
brackets in (III.3). It also allows us to expand the Bessel function
$J_{\nu}(\xi)\approx\frac{1}{\nu!}(\xi/2)^{\nu}$ and retain only the leading
$\nu=1$ term. Since $|m|>1$ in this section, $\cos\theta$ is not restricted at
all in this approximation. The angular integration becomes trivial and yields

$\displaystyle
W^{C}=\frac{e^{2}M^{4}}{192\pi\Omega^{2}}\frac{(2|m|+1)^{2}}{|m|^{3}(|m|+1)^{3}}\,.$
(40)

###  III.4 Summary of $\Omega\gg M$ approximation.

In summary, the radiation intensity at $\Omega\gg M$ is given by

$\displaystyle
W=\frac{e^{2}\Omega^{2}}{8\pi}\left\\{\begin{array}[]{lc}S_{1}(m)+\frac{1}{4}S_{2}(m)+m^{2}J_{m}^{2}(2m)\log\frac{4m\Omega}{M}\,,&m>0\,,\\\
\frac{M^{2}}{4\Omega^{2}}\,,&m=0\,,\\\
\frac{M^{4}}{24\Omega^{4}}\frac{(2|m|+1)^{2}}{|m|^{3}(|m|+1)^{3}}\,,&m<0\,.\end{array}\right.$
(44)

The transitions from $m>0$, with the corresponding energy $E\approx 2\Omega
m$, give the largest contribution. When $m$ is not very large, the leading
channel is the transition to $m^{\prime}=0$ ($E^{\prime}\approx M\ll E$) given
by (33). This is also seen in Fig. 3. For large $m$, we can use the well-known
formula (see e.g. 9.3.15 in [41])

$\displaystyle
J_{m}(2m)\approx\sqrt{\frac{2}{\pi\sqrt{3}m}}\cos\left[m(\sqrt{3}-\pi/3)-\pi/4\right]\,,\quad
m\gg 1\,,$ (45)

to conclude that the contribution from the transitions to $m^{\prime}>0$
becomes dominant since $S_{2}(m)\sim m^{2}$. We verified that (44) is an
accurate approximation of the exact formula at $p_{z}=0$ and $M\ll\Omega$
(plotted in Fig. 2).

##  IV Radiation by spinning ideal gas

###  IV.1 Energy spectrum

If an ideal Maxwell-Boltzmann gas rotates extremely rapidly, then its
radiation intensity is given by

$\displaystyle I$
$\displaystyle=\sum_{m=-\infty}^{\infty}\int\frac{dp_{z}L}{2\pi}\int
do\int_{0}^{E}d\omega\frac{dW}{dod\omega}e^{-\beta E}\,.$ (46)

Using (10),(13) and (14) we obtain

$\displaystyle\frac{dI}{Ld\omega}=\frac{e^{2}}{16\pi}$
$\displaystyle\sum_{m=-\infty}^{\infty}\sum_{\pm}\sum_{m^{\prime}=-\infty}^{m-1}\int\frac{dp_{z}}{2\pi}e^{-\beta
E}\frac{E^{\prime}-m^{\prime}\Omega}{EE^{\prime}\omega|\omega\cos\theta_{\pm}-pz|}$
$\displaystyle\times\left[\frac{(m+m^{\prime})^{2}}{R^{2}}+\sin^{2}\theta_{\pm}(2p_{z}-\omega\cos\theta_{\pm})^{2}\right]J_{m-m^{\prime}}^{2}(\omega
R\sin\theta_{\pm})\omega^{2}\,.$ (47)

where the integral over $p_{z}$ is restricted by (3) for $m<0$. The
corresponding energy spectra are exhibited in Fig. 5. We observe that there is
the threshold frequency below which photon emission is impossible. We also
notice that the spectrum peaks at $\omega\sim\Omega$.

|
---|---
Figure 5: (Color online) The energy spectrum of electromagnetic radiation by a
thermal system. Right panel zooms into the infrared region of the left one.
$M=1$.

###  IV.2 Photon production in heavy-ion collisions

It is tempting to apply the model developed in this article to describe the
photon production by rotating quark-gluon plasma. We certainly realize that
the plasma is not yet in the regime of extremely fast rotation
$\Omega>\ell^{-1}$. On the other hand, it is already in the region of fast
rotation $\Omega\sim a^{-1}$, using the notation introduced in Introduction.
We therefore consider the calculation of the photon spectrum as the no more
than a back-of-the-envelop estimate and the proof of principle that rotation
is a relevant effect.

The photon emission from quark-gluon plasma, is described in terms of the
following variables: $k_{T}$, the photon momentum in the plane perpendicular
to the collision axis (not to be confused with $k_{\bot}$ defined with respect
to the rotation axis), $\phi$ its azimuthal angle in that plane and $y$ its
rapidity. They are related to the photon energy $\omega$ and the emission
angle $\theta$ (see [42] for more details):

$\displaystyle\omega=k_{T}\cosh y,\quad\cos\theta=\sin\phi/\cosh y\,.$ (48)

We now express the photon spectrum as

$\displaystyle\frac{dN(k_{T},y,\phi)}{k_{T}dk_{T}d\phi dy}=g\Delta
t\frac{dI(\omega,\theta)}{\omega^{2}d\omega do}\,,$ (49)

where $\Delta t$ is the time interval. In realistic case photons are emitted
by fermions with two possible polarizations and three colors and three flavors
which make up the degeneracy factor $g=18$. Eq. (49) simplifies at the
midrapidity region $y=0$:

$\displaystyle\frac{dN}{k_{T}dk_{T}dy}\bigg{|}_{y=0}=$
$\displaystyle\frac{g\Delta
tL}{k_{T}^{2}}\int_{0}^{2\pi}d\phi\frac{dI}{Ld\omega
do}\bigg{|}_{\omega=k_{T},\,\theta=\frac{\pi}{2}-\phi}\,.$ (50)

The intensity $I$ is given by (46) and $W$ is given by (10).

It is now advantageous to use the delta-function in (10) to take the integral
over $p_{z}$. To this end we write

$\displaystyle\delta(E-E^{\prime}-\omega)=\delta\left(\sqrt{p_{z}^{2}+m^{2}\Omega^{2}+M^{2}}-\sqrt{(p_{z}-\omega\cos\theta)^{2}+m^{\prime
2}\Omega^{2}+M^{2}}+\Delta\right)$
$\displaystyle=\frac{\delta(p_{z}-p_{z0})(E-m\Omega)(E-\omega-m^{\prime}\Omega)}{|\omega\cos\theta(E-m\Omega)+p_{z}\Delta|}\,,$
(51)

where we introduced a convenient notation

$\displaystyle\Delta=(m-m^{\prime})\Omega-\omega\,.$ (52)

To compute $p_{z0}$ we rewrite the equation in the second delta-function in
(IV.2) as a quadratic equation for $p_{z}$. Of course not all its solutions
necessarily satisfy the original equation. A careful examination of its two
roots reveals that one root satisfies it at $\Delta>0$, while another one at
$\Delta<0$. These can be combined in a single formula:

$\displaystyle p_{z0}=$
$\displaystyle\frac{1}{2(\Delta^{2}-\omega^{2}\cos^{2}\theta)}\Big{\\{}-\omega\cos\theta(\omega^{2}\cos^{2}\theta-\Delta^{2}+m^{2}\Omega^{2}-m^{\prime
2}\Omega^{2})$
$\displaystyle+\Delta\sqrt{(\omega^{2}\cos^{2}\theta-\Delta^{2}+m^{2}\Omega^{2}-m^{\prime
2}\Omega^{2})^{2}+4(m^{2}\Omega^{2}+M^{2})(\omega^{2}\cos^{2}\theta-\Delta^{2})}\Big{\\}}\,,$
(53)

provided that

$\displaystyle\Delta^{2}-\omega^{2}\cos^{2}\theta\leq 0\,.$ (54)

Using (52) in (54) and noting that $|\cos\theta|\leq 1$ we find that the
allowed photon energies are

$\displaystyle\omega\geq\frac{1}{2}(m-m^{\prime})\Omega\,.$ (55)

In particular, the photon spectrum has an infrared threshold at
$\omega_{\text{min}}=\Omega/2$. This is indeed clearly seen in Fig. 5. Eq.
(54) is a constraint on the allowed values of the photon emission angle:

$\displaystyle|\theta|\leq\Theta=\arcsin\sqrt{1-\Delta^{2}/\omega^{2}}\,.$
(56)

Taking all these into account we obtain the final expression for the photon
spectrum:

$\displaystyle\frac{dN}{k_{T}dk_{T}dy}\bigg{|}_{y=0}=$
$\displaystyle\frac{g\Delta
tLe^{2}}{4(2\pi)^{3}}\sum_{m=-\infty}^{\infty}\sum_{m^{\prime}=-\infty}^{m-1}\int_{-\Theta}^{+\Theta}d\theta\frac{1}{EE^{\prime}}e^{-E/T}\frac{(E-m\Omega)(E^{\prime}-m^{\prime}\Omega)}{|k_{T}\cos\theta(E-m\Omega)+p_{z}\Delta|}$
$\displaystyle\times\left[(m+m^{\prime})^{2}\Omega^{2}+\sin^{2}\theta(2p_{z0}-k_{T}\cos\theta)^{2}\right]J^{2}_{m-m^{\prime}}(k_{T}\sin\theta\Omega^{-1})\,,$
(57)

valid for $k_{T}=\omega$ satisfying (55). $E$ and $E^{\prime}$ are the
functions of $p_{z0}$. At negative $m$ and $m^{\prime}$ the angular
integration is further restricted by (3). However, as we have seen, the
negative $m$ and $m^{\prime}$ give a negligible contribution. Therefore in
practice the sums over $m$ and $m^{\prime}$ rum only over the non-negative
values.

Figure 6: (Color online) Photon spectrum at two temperatures and $\Omega=0.1$
fm-1. The data is from [43]. $\Delta t=10$ fm/$c$, $L=10$ fm, $M=150$ MeV.

The results of the calculation are shown in Fig. 6 along with the experimental
data on photon production in relativistic heavy-ion collisions [43]. We
observe a wondrous agreement for a particular plasma temperature. This does
not of course imply that plasma rotation is the only source of photons.
However, it does indicate that it must be accounted for in comprehensive
phenomenological models. And perhaps it will help resolving the long-standing
puzzle of the photon excess at low momenta.

##  V Summary

We posited that systems rotating with extreme angular velocity
$\Omega\sim\ell^{-1}$, where $\ell$ is the mean-free-path, can be described by
matter fields confined to the two dimensional cylindrical surface of radius
$R=1/\Omega$. We developed an application of this idea to the photon
production by the rotating quark-gluon plasma and argued that it is consistent
with the experimental observations.

The statistical properties inside the light-cylinder are equivalent to those
of the two-dimensional ideal gas, as the interactions are screened by
$R<\ell$. However, the statistical properties of the entire plasma are
determined by the interaction of the light-cylinders; the development of these
ideas is left for another study.

###### Acknowledgements.

This work was supported in part by the U.S. Department of Energy Grants No.
DE-SC0023692.

## References

  * Adamczyk  _et al._ [2017] L. Adamczyk _et al._ (STAR), Nature 548, 62 (2017), arXiv:1701.06657 [nucl-ex] .
  * Adam  _et al._ [2018] J. Adam _et al._ (STAR), Phys. Rev. C 98, 014910 (2018), arXiv:1805.04400 [nucl-ex] .
  * Adam  _et al._ [2019] J. Adam _et al._ (STAR), Phys. Rev. Lett. 123, 132301 (2019), arXiv:1905.11917 [nucl-ex] .
  * Acharya  _et al._ [2020] S. Acharya _et al._ (ALICE), Phys. Rev. Lett. 125, 012301 (2020), arXiv:1910.14408 [nucl-ex] .
  * Adam  _et al._ [2021] J. Adam _et al._ (STAR), Phys. Rev. Lett. 126, 162301 (2021), arXiv:2012.13601 [nucl-ex] .
  * Abdallah  _et al._ [2021] M. S. Abdallah _et al._ (STAR), Phys. Rev. C 104, L061901 (2021), arXiv:2108.00044 [nucl-ex] .
  * Zubarev  _et al._ [1979] D. N. Zubarev, A. V. Prozorkevich, and S. A. Smolyanskii, Theoretical and Mathematical Physics 40, 821 (1979).
  * Vilenkin [1980] A. Vilenkin, Phys. Rev. D21, 2260 (1980).
  * Hayata  _et al._ [2015] T. Hayata, Y. Hidaka, T. Noumi, and M. Hongo, Phys. Rev. D 92, 065008 (2015), arXiv:1503.04535 [hep-ph] .
  * Ambrus and Winstanley [2016] V. E. Ambrus and E. Winstanley, Phys. Rev. D 93, 104014 (2016), arXiv:1512.05239 [hep-th] .
  * Buzzegoli  _et al._ [2017] M. Buzzegoli, E. Grossi, and F. Becattini, JHEP 10 (2017), 091, [Erratum: JHEP 07, 119 (2018)], arXiv:1704.02808 [hep-th] .
  * Becattini  _et al._ [2019] F. Becattini, M. Buzzegoli, and E. Grossi, Particles 2, 197 (2019), arXiv:1902.01089 [cond-mat.stat-mech] .
  * Montenegro and Torrieri [2020] D. Montenegro and G. Torrieri, Phys. Rev. D 102, 036007 (2020), arXiv:2004.10195 [hep-th] .
  * Weickgenannt  _et al._ [2019] N. Weickgenannt, X.-L. Sheng, E. Speranza, Q. Wang, and D. H. Rischke, Phys. Rev. D 100, 056018 (2019), arXiv:1902.06513 [hep-ph] .
  * Gao and Liang [2019] J.-H. Gao and Z.-T. Liang, Phys. Rev. D 100, 056021 (2019), arXiv:1902.06510 [hep-ph] .
  * Hattori  _et al._ [2019] K. Hattori, Y. Hidaka, and D.-L. Yang, Phys. Rev. D 100, 096011 (2019), arXiv:1903.01653 [hep-ph] .
  * Bhadury  _et al._ [2021] S. Bhadury, W. Florkowski, A. Jaiswal, A. Kumar, and R. Ryblewski, Phys. Rev. D 103, 014030 (2021), arXiv:2008.10976 [nucl-th] .
  * Liu  _et al._ [2020] Y.-C. Liu, K. Mameda, and X.-G. Huang, Chin. Phys. C 44, 094101 (2020), [Erratum: Chin.Phys.C 45, 089001 (2021)], arXiv:2002.03753 [hep-ph] .
  * Weickgenannt  _et al._ [2021] N. Weickgenannt, E. Speranza, X.-l. Sheng, Q. Wang, and D. H. Rischke, Phys. Rev. Lett. 127, 052301 (2021), arXiv:2005.01506 [hep-ph] .
  * Shi  _et al._ [2021] S. Shi, C. Gale, and S. Jeon, Phys. Rev. C 103, 044906 (2021), arXiv:2008.08618 [nucl-th] .
  * Weickgenannt  _et al._ [2022] N. Weickgenannt, D. Wagner, E. Speranza, and D. H. Rischke, Phys. Rev. D 106, 096014 (2022), arXiv:2203.04766 [nucl-th] .
  * Hidaka  _et al._ [2018] Y. Hidaka, S. Pu, and D.-L. Yang, Phys. Rev. D 97, 016004 (2018), arXiv:1710.00278 [hep-th] .
  * Karpenko and Becattini [2017] I. Karpenko and F. Becattini, Eur. Phys. J. C 77, 213 (2017), arXiv:1610.04717 [nucl-th] .
  * Xie  _et al._ [2017] Y. Xie, D. Wang, and L. P. Csernai, Phys. Rev. C 95, 031901 (2017), arXiv:1703.03770 [nucl-th] .
  * Wu  _et al._ [2019] H.-Z. Wu, L.-G. Pang, X.-G. Huang, and Q. Wang, Phys. Rev. Research. 1, 033058 (2019), arXiv:1906.09385 [nucl-th] .
  * Ivanov [2021] Y. B. Ivanov, Phys. Rev. C 103, L031903 (2021), arXiv:2012.07597 [nucl-th] .
  * Fu  _et al._ [2021] B. Fu, K. Xu, X.-G. Huang, and H. Song, Phys. Rev. C 103, 024903 (2021), arXiv:2011.03740 [nucl-th] .
  * Li  _et al._ [2017] H. Li, L.-G. Pang, Q. Wang, and X.-L. Xia, Phys. Rev. C 96, 054908 (2017), arXiv:1704.01507 [nucl-th] .
  * Wei  _et al._ [2019] D.-X. Wei, W.-T. Deng, and X.-G. Huang, Phys. Rev. C 99, 014905 (2019), arXiv:1810.00151 [nucl-th] .
  * Ivanov  _et al._ [2019] Y. B. Ivanov, V. D. Toneev, and A. A. Soldatov, Phys. Rev. C 100, 014908 (2019), arXiv:1903.05455 [nucl-th] .
  * Guo  _et al._ [2021] Y. Guo, J. Liao, E. Wang, H. Xing, and H. Zhang, Phys. Rev. C 104, L041902 (2021), arXiv:2105.13481 [nucl-th] .
  * Alzhrani  _et al._ [2022] S. Alzhrani, S. Ryu, and C. Shen, Phys. Rev. C 106, 014905 (2022), arXiv:2203.15718 [nucl-th] .
  * Florkowski  _et al._ [2019] W. Florkowski, A. Kumar, and R. Ryblewski, Prog. Part. Nucl. Phys. 108, 103709 (2019), arXiv:1811.04409 [nucl-th] .
  * Becattini  _et al._ [2021] F. Becattini, J. Liao, and M. Lisa, eds.,  _Lecture Notes in Physics_, Vol. 987 (Springer International Publishing, 2021) title: Strongly Interacting Matter under Rotation.
  * Becattini and Lisa [2020] F. Becattini and M. A. Lisa, Annual Review of Nuclear and Particle Science 70, 395 (2020), arXiv:2003.03640 [nucl-ex] .
  * Buzzegoli  _et al._ [2023a] M. Buzzegoli, J. D. Kroth, K. Tuchin, and N. Vijayakumar, Phys. Rev. D 107, L051901 (2023a), arXiv:2209.02597 [hep-ph] .
  * Buzzegoli  _et al._ [2023b] M. Buzzegoli, J. D. Kroth, K. Tuchin, and N. Vijayakumar, arXiv preprint (2023b), arXiv:2306.03863 [hep-ph] .
  * Buzzegoli and Tuchin [2023a] M. Buzzegoli and K. Tuchin, arXiv preprint (2023a), arXiv:2305.13149 [hep-ph] .
  * Buzzegoli and Tuchin [2023b] M. Buzzegoli and K. Tuchin, Nucl. Phys. A 1030, 122577 (2023b), arXiv:2209.03991 [nucl-th] .
  * Sokolov  _et al._ [1986] A. A. Sokolov, I. M. Ternov, and C. W. Kilmister,  _Radiation from relativistic electrons_ (American Institute of Physics, 1986).
  * Abramowitz and Stegun [1968] M. Abramowitz and I. A. Stegun,  _Handbook of mathematical functions with formulas, graphs, and mathematical tables_ , Vol. 55 (US Government printing office, 1968).
  * Tuchin [2015] K. Tuchin, Phys. Rev. C 91, 014902 (2015), arXiv:1406.5097 [nucl-th] .
  * Adare  _et al._ [2015] A. Adare _et al._ (PHENIX), Phys. Rev. C 91, 064904 (2015), arXiv:1405.3940 [nucl-ex] .


# GMM-based Codebook Construction and
Feedback Encoding in FDD Systems††thanks: ©This work has been submitted to the
IEEE for possible publication. Copyright may be transferred without notice,
after which this version may no longer be accessible.

Nurettin Turan∗, Michael Koller∗, Benedikt Fesl∗, Samer Bazzi†, Wen Xu†, and
Wolfgang Utschick∗

∗Professur für Methoden der Signalverarbeitung, Technische Universität
München, 80333 Munich, Germany
†Huawei Technologies Duesseldorf GmbH, 80992 Munich, Germany
Email<EMAIL_ADDRESS>
Email<EMAIL_ADDRESS>
###### Abstract

We propose a precoder codebook construction and feedback encoding scheme which
is based on Gaussian mixture models (GMMs). In an offline phase, the base
station (BS) first fits a GMM to uplink (UL) training samples. Thereafter, it
designs a codebook in an unsupervised manner by exploiting the GMM’s
clustering capability. We design one codebook entry per GMM component. After
offloading the GMM—but not the codebook—to the mobile terminal (MT) in the
online phase, the MT utilizes the GMM to determine the best fitting codebook
entry. To this end, no channel estimation is necessary at the MT. Instead, the
MT’s observed signal is used to evaluate how responsible each component of the
GMM is for the signal. The feedback consists of the index of the GMM component
with the highest responsibility and the BS then employs the corresponding
codebook entry. Simulation results show that the proposed codebook design and
feedback encoding scheme outperforms conventional Lloyd clustering based
codebook design algorithms, especially in configurations with reduced pilot
overhead.

###### Index Terms:

Gaussian mixture models, machine learning, feedback, codebook design,
frequency division duplexing

##  I Introduction

In multiple-input multiple-output (MIMO) communications systems, channel state
information (CSI) has to be acquired at the BS in regular time intervals. In
frequency division duplex (FDD) mode, the BS and the MT transmit in the same
time slot but at different frequencies. This breaks the reciprocity between
the instantaneous UL CSI and downlink (DL) CSI. Accordingly, acquiring DL CSI
in FDD operation is difficult [1]. Possible solutions include to either
extrapolate the DL CSI from the estimate of the UL CSI at the BS, or to
transfer the DL CSI estimated at the MT to the BS directly or in a highly
compressed version. However, the most common solution in practice is to avoid
the direct feedback of the CSI and to use only a small number of feedback
bits. For instance, as it is done in this paper, the feedback can be used as
an index into a predefined codebook of precoders [2].

In recent years, machine learning techniques have been explored in the context
of communications. However, these typically need a large dataset of DL
channels for their training phases. This would require the users to collect
large amounts of DL CSI and either to perform the training themselves or to
share the collected data with the BS. The corresponding computation and/or
signaling overhead involved in such a scheme is generally unaffordable in
practice.

Recently, in [3] it has been shown that DL CSI training data can be replaced
with UL CSI training data even for the design of DL functionalities. This
completely eliminates the aforementioned overhead. UL CSI can, e.g., be
acquired at the BS during the regular UL transmission. In [3], the observation
has been made in the context of training autoencoders. Similar observations
have since been made in [4] for DL channel estimation and in [5, 6] for
codebook design. In this work, we also utilize the idea to centrally learn DL-
related functionalities at the BS using UL training data.

The contributions of this work are summarized as follows. We propose a
codebook construction and feedback encoding scheme which is based on GMMs.
Since GMMs are universal approximators [7], we can use a $K$-components GMM to
approximate the unknown channel probability density function (PDF). In the
_offline_ phase, we propose to fit the GMM centrally at the BS using solely UL
data. Thereafter, we cluster the training data according to the GMM components
and design a codebook entry for every component. In this way, the codebook is
designed for a whole scenario, i.e., for the whole site in which the BS is
located. The GMM is offloaded to every MT in the coverage area of the BS, with
which the MT can then select the best fitting codebook entry in the _online_
phase.

In conventional approaches, first the DL CSI is estimated and subsequently the
best fitting codebook entry is determined. The technique proposed in this
paper allows to bypass explicit CSI estimation and yields the best fitting
codebook entry by using the receive signal to evaluate the GMM
responsibilities. The $k$th responsibility of a GMM corresponds to the
probability that the $k$th GMM component is responsible for the receive
signal. The feedback then consists of the index $k$ of the component with the
highest responsibility. Finally, we make use of a Kronecker approximation to
significantly reduce the number of GMM parameters such that the offloading
overhead is smaller. In simulations, the proposed codebook design and feedback
encoding scheme outperforms conventional Lloyd clustering based codebook
design algorithms [8, 9].

##  II System Model

The DL received signal of a point-to-point MIMO system can be expressed as
$\bm{y}^{\prime}=\bm{H}\bm{x}+\bm{n}^{\prime}$, where
$\bm{y}^{\prime}\in\mathbb{C}^{N_{\mathrm{rx}}}$ is the receive vector,
$\bm{x}\in\mathbb{C}^{N_{\mathrm{tx}}}$ is the transmit vector sent over the
MIMO channel $\bm{H}\in\mathbb{C}^{N_{\mathrm{rx}}\times N_{\mathrm{tx}}}$,
and
$\bm{n}^{\prime}\sim\mathcal{N}_{\mathbb{C}}(\mathbf{0},\sigma_{n}^{2}\bm{I}_{N_{\mathrm{rx}}})$
denotes the additive white Gaussian noise (AWGN). In this paper, we consider
system configurations with $N_{\mathrm{rx}}<N_{\mathrm{tx}}$. The BS is
equipped with a uniform rectangular array (URA) and the MT is equipped with a
uniform linear array (ULA). If perfect CSI is known to both the transmitter
and receiver, and assuming input data with Gaussian distribution, the capacity
of the MIMO channel is [10, 11]:

$C=\max_{\bm{Q}\succeq{\bm{0}},\operatorname{tr}\bm{Q}\leq\rho}\log_{2}\det\left(\bm{I}+\frac{1}{\sigma_{n}^{2}}\bm{H}\bm{Q}\bm{H}^{\operatorname{H}}\right),$
(1)

where $\bm{Q}\in\mathbb{C}^{N_{\mathrm{tx}}\times N_{\mathrm{tx}}}$ is the
transmit covariance matrix and the transmit vector is then given by
$\bm{x}=\bm{Q}^{1/2}\bm{s}$ with
$\operatorname{E}[\bm{s}\bm{s}^{\operatorname{H}}]=\bm{I}_{N_{\mathrm{tx}}}$
[2]. The optimal transmit covariance matrix $\bm{Q}^{\star}$ of the link
between the BS and a MT achieves the capacity and can be obtained by
decomposing the channel into $N_{\mathrm{rx}}$ parallel streams and employing
water-filling [12].

Channel reciprocity can generally not be assumed in FDD systems, e.g., [2].
Therefore, only the MT could compute the optimal transmit covariance matrix
$\bm{Q}^{\star}$ if it estimated the DL CSI. This makes some form of feedback
from the MT to the BS necessary. Ideally, the user would feed the complete DL
CSI back to the BS, which is infeasible in general. Instead, typically limited
feedback is considered where a small number of $B$ bits is sent back to the
BS. The $B$ feedback bits are typically used for encoding an index that
specifies an element from a set of covariance matrices. That is, the MT and BS
share a codebook
$\mathcal{Q}=\\{\bm{Q}_{1},\bm{Q}_{2},\dots,\bm{Q}_{2^{B}}\\}$ of $2^{B}$ pre-
computed transmit covariance matrices, and the MT is assumed to estimate the
DL channel $\bm{H}$ and then uses it to determine the best codebook entry
$\bm{Q}_{k^{\star}}$ via

$k^{\star}=\operatorname*{arg\,max}_{k\in\\{1,\dots,2^{B}\\}}\log_{2}\det\left(\bm{I}+\frac{1}{\sigma_{n}^{2}}\bm{H}\bm{Q}_{k}\bm{H}^{\operatorname{H}}\right).$
(2)

The feedback consists of the index $k^{\star}$ encoded by $B$ bits, and the BS
employs the transmit covariance matrix $\bm{Q}_{k^{\star}}$ for data
transmission.

##  III Channel Model and Data Generation

Version $2.6.1$ of the QuaDRiGa channel simulator [13, 14] is used to generate
CSI for the UL and DL domains in an urban macrocell (UMa) scenario. The
carrier frequencies are $2.53\text{\,}\mathrm{GHz}$ for the UL and
$2.73\text{\,}\mathrm{GHz}$ for the DL such that there is a frequency gap of
$200\text{\,}\mathrm{MHz}$. The BS uses a URA with “3GPP-3D” antennas, and the
MTs use ULAs with “omni-directional” antennas. The BS covers a
$120\text{\,}\mathrm{\SIUnitSymbolDegree}$ sector and is placed at
$25\text{\,}\mathrm{m}$ height. The minimum and maximum distances between MTs
and the BS are $35\text{\,}\mathrm{m}$ and $500\text{\,}\mathrm{m}$,
respectively. In $80\%$ of the cases, the MTs are located indoors at different
floor levels. The outdoor MTs have a height of $1.5\text{\,}\mathrm{m}$.

Many parameters such as path-loss, delay, and angular spreads, path-powers for
each subpath, and antenna patterns are different in the DL and UL domain [13].
However, the following parameters are identical: BS location and the MT
locations, propagation cluster delays and angles for each multi-path component
(MPC), and the spatial consistency of the large scale fading parameters. A
QuaDRiGa MIMO channel is given by
$\bm{H}=\sum_{\ell=1}^{L}\bm{G}_{\ell}e^{-2\pi jf_{c}\tau_{\ell}}$ with $\ell$
the path number, $L$ the number of MPCs, $f_{c}$ the carrier frequency, and
$\tau_{\ell}$ the $\ell$th path delay. The number $L$ depends on whether there
is line of sight (LOS), non-line of sight (NLOS), or outdoor-to-indoor (O2I)
propagation: $L_{\text{LOS}}=37$, $L_{\text{NLOS}}=61$ or $L_{\text{O2I}}=37$.
The coefficients matrix $\bm{G}_{\ell}$ consists of one complex entry for each
antenna pair and comprises the attenuation of a path, the antenna radiation
pattern weighting, and the polarization. As described in the QuaDRiGa manual
[14], the generated channels are post-processed to remove the path gain.

##  IV Channel Estimation

In the pilot transmission phase, the DL received signal is:

$\bm{Y}=\bm{H}\bm{P}+\bm{N}\in\mathbb{C}^{N_{\mathrm{rx}}\times n_{p}},$ (3)

where $n_{p}$ is the number of transmitted pilots and
$\bm{N}=[\bm{n}^{\prime}_{1},\bm{n}^{\prime}_{2},\dots,\bm{n}^{\prime}_{n_{p}}]\in\mathbb{C}^{N_{\mathrm{rx}}\times
n_{p}}$. The pilot matrix $\bm{P}\in\mathbb{C}^{N_{\mathrm{tx}}\times n_{p}}$
is a $2$D-DFT (sub)matrix, constructed by the Kronecker product of two
discrete Fourier transform (DFT) matrices,
$\bm{P}=\bm{P}_{\text{h}}\otimes\bm{P}_{\text{v}}$, where each column
$\bm{p}_{p}$ of $\bm{P}$, for $p\in\\{1,2,\dots,n_{p}\\}$, is normalized such
that $\|\bm{p}_{p}\|^{2}=\rho$, since we employ a URA at the BS [15]. In this
work, we consider $n_{p}\leq N_{\mathrm{tx}}$. For what follows, it is
convenient to vectorize (3):

$\bm{y}=\bm{A}\bm{h}+\bm{n},$ (4)

with the definitions $\bm{h}=\operatorname{vec}(\bm{H})$,
$\bm{y}=\operatorname{vec}(\bm{Y})$, $\bm{n}=\operatorname{vec}(\bm{N})$,
$\bm{A}=\bm{P}^{\operatorname{T}}\otimes\bm{I}_{N_{\mathrm{rx}}}$ and
$\bm{n}\sim\mathcal{N}_{\mathbb{C}}(\mathbf{0},{\bm{\Sigma}}=\sigma_{n}^{2}\bm{I}_{n_{p}\times
N_{\mathrm{rx}}})$.

###  IV-A GMM based Channel Estimation

A GMM is a PDF of the form [16]

$f^{(K)}_{\bm{h}}(\bm{h})=\sum_{k=1}^{K}p(k)\mathcal{N}_{\mathbb{C}}(\bm{h};{\bm{\mu}}_{k},\bm{C}_{k})$
(5)

where every summand is one of its $K$ components. Maximum likelihood estimates
of the parameters of a GMM, i.e, the means ${\bm{\mu}}_{k}$, the covariances
$\bm{C}_{k}$, and the mixing coefficients $p(k)$, can be computed using a
training data set $\mathcal{H}$ and an expectation maximization (EM)
algorithm, see [16]. GMMs allow for the evaluation of responsibilities [16]:

$p(k\mid\bm{h})=\frac{p(k)\mathcal{N}_{\mathbb{C}}(\bm{h};{\bm{\mu}}_{k},\bm{C}_{k})}{\sum_{i=1}^{K}p(i)\mathcal{N}_{\mathbb{C}}(\bm{h};{\bm{\mu}}_{i},\bm{C}_{i})}.$
(6)

These corresponds to the probability that a given $\bm{h}$ was drawn from
component $k$. Interestingly, GMMs can approximate any continuous PDF
arbitrary well [7].

With this GMM background, we now briefly recap the GMM channel estimator from
[17]. Given a training data set of channels
$\mathcal{H}=\\{\bm{h}_{m}=\operatorname{vec}{(\bm{H}_{m})}\\}_{m=1}^{M}$, the
EM algorithm is used to compute a $K$-component GMM $f_{\bm{h}}^{(K)}$ as an
approximation of the true but unknown channel PDF $f_{\bm{h}}$.

The idea in [17, 18] now is to compute the mean squared error (MSE)-optimal
estimator $\hat{\bm{h}}_{\text{GMM}}^{(K)}$ for channels distributed according
to $f_{\bm{h}}^{(K)}$ and to use it to estimate the channels distributed
according to $f_{\bm{h}}$. The motivation for this is that
$\hat{\bm{h}}_{\text{GMM}}^{(K)}$ converges pointwise to the MSE-optimal
estimator for channels distributed according to $f_{\bm{h}}$ as $K\to\infty$
[17, 18].

The estimator $\hat{\bm{h}}_{\text{GMM}}^{(K)}$ can be computed in closed
form:

$\hat{\bm{h}}_{\text{GMM}}^{(K)}(\bm{y})=\sum_{k=1}^{K}p(k\mid\bm{y})\hat{\bm{h}}_{\text{LMMSE},k}(\bm{y})$
(7)

with the responsibilities

$p(k\mid\bm{y})=\frac{p(k)\mathcal{N}_{\mathbb{C}}(\bm{y};\bm{A}{\bm{\mu}}_{k},\bm{A}\bm{C}_{k}\bm{A}^{\operatorname{H}}+{\bm{\Sigma}})}{\sum_{i=1}^{K}p(i)\mathcal{N}_{\mathbb{C}}(\bm{y};\bm{A}{\bm{\mu}}_{i},\bm{A}\bm{C}_{i}\bm{A}^{\operatorname{H}}+{\bm{\Sigma}})}$
(8)

and

$\hat{\bm{h}}_{\text{LMMSE},k}(\bm{y})=\bm{C}_{k}\bm{A}^{\operatorname{H}}(\bm{A}\bm{C}_{k}\bm{A}^{\operatorname{H}}+{\bm{\Sigma}})^{-1}(\bm{y}-\bm{A}{\bm{\mu}}_{k})+{\bm{\mu}}_{k}.$
(9)

The estimator $\hat{\bm{h}}_{\text{GMM}}^{(K)}$ calculates a weighted sum of
$K$ linear minimum mean square error (LMMSE) estimators—one for each
component. The weights $p(k\mid\bm{y})$ are the probabilities that the current
observation $\bm{y}$ corresponds to the $k$th component.

###  IV-B Baseline Channel Estimators

As a first baseline, we consider a sample covariance matrix based channel
estimation approach, where we construct a sample covariance matrix
$\bm{C}_{s}=\frac{1}{M}\sum_{m=1}^{M}\bm{h}_{m}\bm{h}_{m}^{\operatorname{H}}$
given the same set of training samples which is used to fit the GMM and
calculate LMMSE channel estimates:

$\hat{\bm{h}}_{\text{s-cov}}=\bm{C}_{s}\bm{A}^{\operatorname{H}}(\bm{A}\bm{C}_{s}\bm{A}^{\operatorname{H}}+{\bm{\Sigma}})^{-1}\bm{y}.$
(10)

Secondly, compressive sensing approaches commonly assume that the channel
exhibits a certain structure: $\bm{h}\approx\bm{D}\bm{t}$, where
$\bm{D}=\bm{D}_{\text{rx}}\otimes(\bm{D}_{\text{tx,h}}\otimes\bm{D}_{\text{tx,v}})$
is a dictionary with oversampled DFT matrices $\bm{D}_{\text{rx}}$,
$\bm{D}_{\text{tx,h}}$ and $\bm{D}_{\text{tx,v}}$ (cf., e.g., [19]), because
we have a URA at the transmitter and a ULA at the receiver side. A compressive
sensing algorithm like orthogonal matching pursuit (OMP) [20] can now be used
to obtain a sparse vector $\bm{t}$, and the estimated channel is then given by

$\hat{\bm{h}}_{\text{OMP}}=\bm{D}\bm{t}.$ (11)

Since the sparsity order is not known but the algorithm’s performance
crucially depends on it, we use a genie-aided approach to obtain a bound on
the performance of the algorithm. Namely, we use the true channel (perfect CSI
knowledge) to choose the optimal sparsity order.

##  V Codebook Design

###  V-A Proposed Codebook Construction and Encoding Scheme

As explained around (7), the first step in computing channel estimates via
$\hat{\bm{h}}_{\text{GMM}}^{(K)}$ consists of determining how likely it is
that the current observation $\bm{y}$ corresponds to the $k$th component of
the GMM $f_{\bm{h}}^{(K)}$, see the responsibility $p(k\mid\bm{y})$ in (8).
The idea of the proposed method now is to compute a codebook transmit
covariance matrix $\bm{Q}_{k}$ for every component of the GMM and to use the
responsiblities $p(k\mid\bm{y})$ to determine the feedback index.

In detail, in an offline training phase, we take $K=2^{B}$ as the number of
GMM components, use a training data set of channels
$\mathcal{H}=\\{\bm{h}_{m}\\}_{m=1}^{M}$ to fit a $K$-component GMM
$f_{\bm{h}}^{(K)}$, and compute a codebook
$\mathcal{Q}=\\{\bm{Q}_{k}\\}_{k=1}^{K}$ of transmit covariance matrices—one
matrix for every GMM component. We explain the codebook construction in
another paragraph below. During the online phase, we bypass explicit channel
estimation and directly determine a feedback index using the responsibilities
computed via $\bm{y}$:

$k^{\star}=\operatorname*{arg\,max}_{k}{p(k\mid\bm{y})}.$ (12)

Thus, we compute the feedback index $k^{\star}$ without requiring (estimated)
CSI. Note, we thereby also avoid the $\log_{2}\det$ evaluation in (2).
Further, the knowledge of the codebook at the MT is not required. The MT only
requires the GMM to compute (12).

We can think of $p(k\mid\bm{y})$ as an approximation of $p(k\mid\bm{h})$ from
(6), because of the fixed noise covariance of every component. That is, since
there is a true underlying channel $\bm{h}$ leading to the current observation
$\bm{y}=\bm{A}\bm{h}+\bm{n}$, $p(k\mid\bm{y})$ can be seen as an approximation
of the probability $p(k\mid\bm{h})$ that the channel $\bm{h}$ was generated
from the $k$th GMM component. To both gauge the influence of using
$p(k\mid\bm{y})$ instead of $p(k\mid\bm{h})$ and to evaluate the codebook
itself, it is interesting to look at the performance of feedback information
calculated as

$k^{\star}=\operatorname*{arg\,max}_{k}{p(k\mid\bm{h})}.$ (13)

Of course, this approach is not practically feasible because the channel
$\bm{h}$ would have to be known.

_Codebook construction:_ Once the training data set
$\mathcal{H}=\\{\bm{h}_{m}\\}_{m=1}^{M}$ has been used to fit a $K$-component
GMM, we cluster the training data according to their GMM reponsibilities. That
is, we partition $\mathcal{H}$ into $K$ disjoint sets

$\mathcal{V}_{k}=\\{\bm{h}\in\mathcal{H}\mid p(k\mid\bm{h})\geq
p(j\mid\bm{h})\text{ for }k\neq j\\}$ (14)

for $k=1,\dots,K$. For a channel matrix $\bm{H}$ and a covariance matrix
$\bm{Q}$, let

$r(\bm{H},\bm{Q})=\log_{2}\det\left(\bm{I}+\frac{1}{\sigma_{n}^{2}}\bm{H}\bm{Q}\bm{H}^{\operatorname{H}}\right)$
(15)

be the spectral efficiency. We now determine the codebook
$\mathcal{Q}=\\{\bm{Q}_{k}\\}_{k=1}^{K}$ by computing every transmit
covariance matrix $\bm{Q}_{k}$ such that it maximizes the summed rate in
$\mathcal{V}_{k}$:

$\displaystyle\bm{Q}_{k}=\operatorname*{arg\,max}_{\bm{Q}\succeq{\bm{0}}}\frac{1}{|\mathcal{V}_{k}|}\sum_{\operatorname{vec}(\bm{H})\in\mathcal{V}_{k}}r(\bm{H},\bm{Q})$
(16) $\displaystyle\text{subject
to}\quad\operatorname{trace}(\bm{Q})\leq\rho\quad\text{and}\quad\operatorname{rank}\bm{Q}\leq
N_{\mathrm{rx}}.$

This optimization problem is solved via projected gradient descent (PGD), cf.
[5, 21].

In summary, the GMM is used twice: Once for codebook construction (done
offline) and thereafter to determine a feedback index (done online). For the
latter, it is not necessary to estimate the channel and evaluating (2) is
avoided.

###  V-B Conventional Codebook Construction Methods

A standard codebook construction approach makes use of Lloyd’s algorithm [8,
9]. Given a training data set of channels
$\mathcal{H}=\\{\bm{h}_{m}\\}_{m=1}^{M}$, the iterative Lloyd clustering
algorithm alternates between two stages until a convergence criterion is met.
We write $\\{\bm{Q}_{k}^{(i)}\\}_{k=1}^{K}$ for the codebook in iteration $i$.
The two stages in iteration $i$ are:

  1. 1.

Divide the training data set $\mathcal{H}$ into $K$ clusters
$\mathcal{V}_{k}^{(i)}$:

$\mathcal{V}_{k}^{(i)}=\\{\bm{h}\in\mathcal{H}\mid
r(\bm{H},\bm{Q}_{k}^{(i)})\geq r(\bm{H},\bm{Q}_{j}^{(i)})\text{ for }k\neq
j\\}.$ (17)

  2. 2.

Update the codebook:

$\displaystyle\bm{Q}_{k}^{(i+1)}=\operatorname*{arg\,max}_{\bm{Q}\succeq{\bm{0}}}\frac{1}{|\mathcal{V}_{k}^{(i)}|}\sum_{\operatorname{vec}(\bm{H})\in\mathcal{V}_{k}^{(i)}}r(\bm{H},\bm{Q})$
(18) $\displaystyle\text{subject
to}\quad\operatorname{trace}(\bm{Q})\leq\rho\quad\text{and}\quad\operatorname{rank}\bm{Q}\leq
N_{\mathrm{rx}}.$

The optimization problem in stage 2) is again solved via PGD. To initialize
the algorithm, stage 1) is replaced with a random partition of $\mathcal{H}$
in the first iteration.

_Lau’s heuristic:_ In order to avoid solving the costly optimization problem
in stage 2) of every iteration, the authors of [9] provide a heuristic for the
codebook update: A representative matrix
$\bm{S}_{k}^{(i)}=\frac{1}{|\mathcal{V}_{k}^{(i)}|}\sum_{\operatorname{vec}(\bm{H})\in\mathcal{V}_{k}^{(i)}}\bm{H}^{\operatorname{H}}\bm{H}$
is calculated for every cluster $\mathcal{V}_{k}^{(i)}$, and then the matrices
$\bm{S}_{k}^{(i)}$ are decomposed into $N_{\mathrm{rx}}$ parallel streams and
water-filling is employed, yielding the updated codebook entries
$\bm{Q}_{k}^{(i+1)}$, see [9].

Analogously, we can replace the optimization problem in (16) with the
described heuristic from [9] to compute a transmit covariance matrix for every
GMM component. However, as the simulation results in Section VII show, the
performance with PGD is better, especially in the high signal-to-noise ratio
(SNR) regime.

##  VI Complexity Analysis

The responsibilities in (6) or (8) are calculated by evaluating Gaussian
densities. A Gaussian density with mean ${\bm{\mu}}\in\mathbb{C}^{N}$ and
covariance matrix $\bm{C}\in\mathbb{C}^{N\times N}$ can be written as

$\mathcal{N}_{\mathbb{C}}(\bm{h};{\bm{\mu}},\bm{C})=\frac{\exp(-(\bm{h}-{\bm{\mu}})^{\operatorname{H}}\bm{C}^{-1}(\bm{h}-{\bm{\mu}}))}{\pi^{N}\det(\bm{C})}.$
(19)

Since the GMM covariance matrices and mean vectors do not change between
observations, the inverse and the determinant of the densities can be pre-
computed once. Therefore, the online evaluation is dominated by matrix-vector
multiplications and has a complexity of $\mathcal{O}(N^{2})$, where
$N=N_{\mathrm{tx}}N_{\mathrm{rx}}$ to evaluate $p(k\mid\bm{h})$ in (6)
(assuming perfect CSI), or $N=n_{p}N_{\mathrm{rx}}$ to evaluate
$p(k\mid\bm{y})$ in (8) using the observations $\bm{y}$.

###  VI-A Kronecker Approximation for Saving Complexity

In order for a MT to be able to compute feedback indices, the parameters of
the GMM $f_{\bm{h}}^{(K)}$ need to be offloaded to the MT upon entering the
BS’s coverage area. As demonstrated in a numerical example below, the number
of GMM parameters can be quite large. This is mainly due to the large number
of parameters of the GMM’s covariance matrices. As a remedy, we can constrain
the GMM covariance matrices to a particular form with less parameters.

For spatial correlation scenarios, a well-known assumption is that the
scattering in the vicinity of the transmitter and of the receiver are
independent of each other, cf. [22]. This assumption leads to channel
covariance matrices $\bm{C}$, which can be decomposed into the Kronecker
product of a transmit and receive side spatial covariance matrix:
$\bm{C}=\bm{C}_{\text{tx}}\otimes\bm{C}_{\text{rx}}$. As in [18], we use this
assumption to construct a GMM consisting of Kronecker product covariance
matrices $\bm{C}_{k}=\bm{C}_{\text{tx},k}\otimes\bm{C}_{\text{rx},k}$.

The procedure suggested thus far is to fit a single GMM using the vectorized
channel training data
$\mathcal{H}=\\{\bm{h}_{m}=\operatorname{vec}(\bm{H}_{m})\\}_{m=1}^{M}$ of
dimension $N=N_{\mathrm{tx}}N_{\mathrm{rx}}$. This results in unconstrained
GMM covariance matrices of dimension $N\times N$. To achieve Kronecker product
covariance matrices, a two stage procedure is used in [18]. First, we fit two
independent “transmit and receive GMMs” with respective covariance matrices of
dimensions $N_{\mathrm{tx}}\times N_{\mathrm{tx}}$ and $N_{\mathrm{rx}}\times
N_{\mathrm{rx}}$. To this end, all rows of $\\{\bm{H}_{m}\\}_{m=1}^{M}$ are
used to fit a $K_{\mathrm{tx}}$-component transmit GMM, and all columns of
$\\{\bm{H}_{m}\\}_{m=1}^{M}$ are used to fit a $K_{\mathrm{rx}}$-component
receive GMM. Thereafter, a $K=K_{\mathrm{tx}}K_{\mathrm{rx}}$-component GMM
with Kronecker covariance matrices of dimension $N\times N$ is obtained by
computing all Kronecker products
$\bm{C}_{\text{tx},i}\otimes\bm{C}_{\text{rx},j}$ of the transmit GMM
covariance matrices $\bm{C}_{\text{tx},i}$ and receive GMM covariance matrices
$\bm{C}_{\text{rx},j}$. Please refer to [18] for more details. The advantages
of the Kronecker GMM are a lower offline training complexity, the ability to
parallelize the fitting process, and the need for fewer training samples since
the Kronecker GMM has much fewer parameters.

_Numerical example:_ To illustrate the difference in the number of GMM
parameters, we plug in the simulation parameters which we consider in Section
VII. There, we have, $N_{\mathrm{tx}}=32$, $N_{\mathrm{rx}}=16$,
$K_{\text{tx}}=16$ and $K_{\text{rx}}=4$, which yields
$N=N_{\mathrm{tx}}N_{\mathrm{rx}}=512$ and
$K=K_{\mathrm{tx}}K_{\mathrm{rx}}=64$. The normal GMM consists of $K=64$
covariance matrices of dimension $N\times N$ which means that it has
$K\frac{N(N+1)}{2}=8404992$ covariance parameters (taking symmetries into
account). By contrast, the Kronecker GMM has only
$K_{\mathrm{rx}}\frac{N_{\mathrm{rx}}(N_{\mathrm{rx}}+1)}{2}+K_{\mathrm{tx}}\frac{N_{\mathrm{tx}}(N_{\mathrm{tx}}+1)}{2}=8992$
covariance parameters. Therefore, with the Kronecker GMM, the number of
parameters which need to be offloaded is drastically reduced. For this reason,
we consider the Kronecker GMM in Section VII.

##  VII Simulation Results

The BS equipped with a URA has in total $N_{\mathrm{tx}}=32$ antenna elements,
with $N_{\mathrm{tx,v}}=4$ vertical and $N_{\mathrm{tx,h}}=8$ horizontal
elements. At the MT we have a ULA with $N_{\mathrm{rx}}=16$. We consider $B=6$
feedback bits and thus $K=2^{6}=64$.

We generate datasets with $30\cdot 10^{3}$ channels for both the UL and DL
domain of the scenario: $\mathcal{H}^{\text{UL}}$ and
$\mathcal{H}^{\text{DL}}$. The UL channels have a dimension of $32\times 16$
and the DL channels have a dimension of $16\times 32$. The data samples are
normalized such that
$\operatorname{E}[\|\bm{h}\|^{2}]=N=N_{\mathrm{tx}}N_{\mathrm{rx}}$ holds for
the vectorized channels. We further set $\rho=1$ which allows us to define the
SNR as $\frac{1}{\sigma_{n}^{2}}$. We split the two sets
$\mathcal{H}^{\text{UL}}$ and $\mathcal{H}^{\text{DL}}$ into a training set
with $M=20\cdot 10^{3}$ samples, and the remaining samples constitute an
evaluation set:
$\mathcal{H}_{\text{train}}^{\text{UL}},\mathcal{H}_{\text{eval}}^{\text{UL}},\mathcal{H}_{\text{train}}^{\text{DL}},\
\text{and}\ \mathcal{H}_{\text{eval}}^{\text{DL}}.$ However, the UL evaluation
set $\mathcal{H}_{\text{eval}}^{\text{UL}}$ is not relevant for our
considerations and the following transmit strategies are always evaluated on
$\mathcal{H}_{\text{eval}}^{\text{DL}}$, i.e., in the DL domain. When we fit
the GMM based on $\mathcal{H}_{\text{train}}^{\text{UL}}$, we transpose all
elements of the set to emulate a DL.

In the following, we depict the normalized spectral efficiency (nSE) as
performance measure. The spectral efficiencies achieved with a given transmit
covariance matrix are normalized by the spectral efficiency achieved with the
optimal transmit covariance matrix which is given by decomposing the channel
into $N_{\mathrm{rx}}$ parallel streams and employing water-filling [12]. The
empirical complementary cumulative distribution function (cCDF)
$P(\text{nSE}>s)$ of the normalized spectral efficiency denoted by the
variable $s$ (the corresponding random variable is simply denoted by nSE), is
used to depict the empirical probability that the nSE exceeds a specific value
$s$.

We consider the following baseline transmit strategies:

_i_) The curves labeled “uni pow cov” represent uniform power allocation where
the transmit covariance matrix is given by
$\bm{Q}=\frac{\rho}{N_{\mathrm{tx}}}\bm{I}$. In this case, no CSI knowledge or
codebook is used.
_ii_) Moreover, “uni pow eigsp” depicts the transmit strategy where a transmit
covariance matrix is calculated by allocating equal power on the eigenvectors
of the channel. That is, the channel is decomposed into $N_{\mathrm{rx}}$
parallel streams and $\frac{\rho}{N_{\mathrm{rx}}}$ power is allocated to each
stream. Note, this approach is infeasible because the BS would require full
knowledge of the DL channel (or its eigenvectors).

$0.3$$0.35$$0.4$$0.45$$0.5$$0.55$$0.6$$0.65$$0.7$$0.75$$0.8$$0.85$$0.9$$0.95$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$normalized
spectral efficiency, $s$$P(\text{nSE}>s)$(a) Codebooks @
$0\text{\,}\mathrm{d}\mathrm{B}$ SNR, evaluated with perfect CSIuni pow covuni
pow eigspLloyd PGD DLLloyd PGD ULLloyd Lau DLLloyd Lau ULGMM PGD DLGMM PGD
ULGMM Lau DLGMM Lau
UL$0.5$$0.55$$0.6$$0.65$$0.7$$0.75$$0.8$$0.85$$0.9$$0.95$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$normalized
spectral efficiency, $s$$P(\text{nSE}>s)$(b) Codebooks @
$10\text{\,}\mathrm{d}\mathrm{B}$ SNR, evaluated with perfect CSI Figure 1:
Empirical cCDFs of the normalized (by the optimal transmit strategy) spectral
efficiencies achieved with different codebooks and transmit strategies
evaluated with perfect CSI.

In Fig. 1(a), we set the $\text{SNR}=$0\text{\,}\mathrm{d}\mathrm{B}$$. The
conventional codebook construction approaches are denoted by “Lloyd PGD UL/DL”
and “Lloyd Lau UL/DL”, depending on whether PGD or Lau’s heuristic is used to
update the codebook in the second stage of the iterative Lloyd clustering
algorithm, and depending on whether $\mathcal{H}_{\text{train}}^{\text{UL}}$
or $\mathcal{H}_{\text{train}}^{\text{DL}}$ is used as training data to
construct the codebooks. With these approaches, the codebook is known to the
BS and the MT and additionally perfect CSI is assumed at the MT. Each user
then selects the best possible codebook entry by evaluating (2). PGD seems to
be slightly better than Lau’s heuristic. Further, using DL or UL training data
results in approximately the same performance.

The proposed codebook construction and encoding scheme is denoted by “GMM PGD
UL/DL” and “GMM Lau UL/DL”, again depending on whether PGD or Lau’s heuristic
is used to construct the codebook. We either use
$\mathcal{H}_{\text{train}}^{\text{UL}}$ or
$\mathcal{H}_{\text{train}}^{\text{DL}}$ as training data to fit the GMM and
to construct the codebook as described in V-A. With our proposed approach, the
knowledge of the codebook at the MT is not required. After offloading the GMM
to the MT and given perfect CSI knowledge, the MT can then simply determine
the feedback index by evaluating (13). Again, PGD is slightly better than
Lau’s heuristic, and using DL or UL training data results in approximately the
same performance. The proposed GMM approach performs slightly worse in
comparison to the conventional Lloyd clustering approach. In Fig. 1(b), we set
$\text{SNR}=$10\text{\,}\mathrm{d}\mathrm{B}$$, and observe similar results.
Interestingly, “GMM Lau UL/DL” performs better than “Lloyd Lau UL/DL”.

However, assuming perfect CSI at the MT is not feasible. In fact, it is
desired to obtain relatively good system performances with estimated CSI,
where typically only a fraction of the number of transmit antennas
$N_{\mathrm{tx}}$ is used as the number of pilots $n_{p}$, i.e., when
considering systems with reduced pilot overhead. In the following, we consider
the proposed codebook construction and encoding scheme and the conventional
Lloyd clustering algorithm exclusively with PGD due to its superior
performance. Additionally, we only consider UL training data in the remainder.

In Fig. 2(a), the $\text{SNR}=$0\text{\,}\mathrm{d}\mathrm{B}$$ and we have
$n_{p}=8$. We depict results for the conventional Lloyd clustering approach,
where we first estimate the channel either via OMP (11), or the sample
covariance approach (10), or via the GMM estimator (7), and then select a
transmit covariance matrix by evaluating (2) given the estimated channel:
“Lloyd PGD, $\hat{\bm{h}}_{\text{OMP}}$”, “Lloyd PGD,
$\hat{\bm{h}}_{\text{s-cov}}$”, and “Lloyd PGD,
$\hat{\bm{h}}_{\text{GMM}}^{(K)}$”. As can be seen, estimating the channel via
the GMM estimator gives the best performance when considering the conventional
approach.

By contrast, with our proposed approach denoted by “GMM PGD, $\bm{y}$”, where
we bypass channel estimation and directly evaluate (12) for determining a
feedback index, we achieve an even better performance as compared to the
conventional approach. With the curves “Lloyd PGD, $\bm{h}$”, and “GMM PGD,
$\bm{h}$” we depict the case of assuming perfect CSI knowledge (this is a
performance bound). A similar observation can also be made in Fig. 2(b), where
the $\text{SNR}=$15\text{\,}\mathrm{d}\mathrm{B}$$ and we only have $n_{p}=4$
pilots.

In Fig. 3(a), we set $\text{SNR}=$0\text{\,}\mathrm{d}\mathrm{B}$$ and in Fig.
3(b) we have $\text{SNR}=$5\text{\,}\mathrm{d}\mathrm{B}$$, we fix $s=0.8$ and
consider $P(\text{nSE}>0.8)$ for a varying number of pilots $n_{p}$. We see,
that our proposed approach is especially beneficial in the low number of
pilots regime and outperforms the conventional approach, which requires both
channel estimation and the evaluation of (2).

$0.3$$0.35$$0.4$$0.45$$0.5$$0.55$$0.6$$0.65$$0.7$$0.75$$0.8$$0.85$$0.9$$0.95$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$normalized
spectral efficiency, $s$$P(\text{nSE}>s)$(a) Codebooks @
$0\text{\,}\mathrm{d}\mathrm{B}$ SNR, evaluated with imperfect CSI @
$n_{p}=8$uni pow covuni pow eigspLloyd PGD, $\bm{h}$Lloyd PGD,
$\hat{\bm{h}}_{\text{s-cov}}$Lloyd PGD, $\hat{\bm{h}}_{\text{OMP}}$Lloyd PGD,
$\hat{\bm{h}}_{\text{GMM}}^{(K)}$GMM PGD, $\bm{h}$GMM PGD,
$\bm{y}$$0.5$$0.55$$0.6$$0.65$$0.7$$0.75$$0.8$$0.85$$0.9$$0.95$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$normalized
spectral efficiency, $s$$P(\text{nSE}>s)$(b) Codebooks @
$15\text{\,}\mathrm{d}\mathrm{B}$ SNR, evaluated with imperfect CSI @
$n_{p}=4$ Figure 2: Empirical cCDFs of the normalized (by the optimal transmit
strategy) spectral efficiencies achieved with different codebooks and transmit
strategies evaluated with imperfect CSI.
$2$$4$$8$$16$$32$$0$$0.2$$0.4$$0.6$$0.8$$1$$n_{p}$$P(\text{nSE}>0.8)$(a) SNR
$=$0\text{\,}\mathrm{d}\mathrm{B}$$Lloyd PGD, $\bm{h}$Lloyd PGD,
$\hat{\bm{h}}_{\text{s-cov}}$Lloyd PGD, $\hat{\bm{h}}_{\text{OMP}}$Lloyd PGD,
$\hat{\bm{h}}_{\text{GMM}}^{(K)}$GMM PGD, $\bm{h}$GMM PGD,
$\bm{y}$$2$$4$$8$$16$$32$$0$$0.2$$0.4$$0.6$$0.8$$1$$n_{p}$$P(\text{nSE}>0.8)$(b)
SNR $=$5\text{\,}\mathrm{d}\mathrm{B}$$ Figure 3: The probability that the nSE
of a certain transmit strategy exceeds $80\%$ of the optimal transmit
strategy’s spectral efficiency for a varying number of pilots.

##  VIII Conclusion and Outlook

We proposed a codebook construction and feedback encoding scheme which is
based on GMMs. The proposed approach involves an offline phase where a GMM is
fitted and a codebook is constructed at the BS using solely UL data. In the
online phase, the same GMM, which is offloaded to a MT upon entering the
coverage area of the BS, was used for feedback encoding. Simulation results
confirmed the validity of this approach, especially in configurations with a
reduced pilot overhead. In future work, we will investigate the GMM based
feedback encoding principle for systems with multiple users.

## References

  * [1] E. Björnson, L. Sanguinetti, H. Wymeersch, J. Hoydis, and T. L. Marzetta, “Massive MIMO is a reality–What is next? five promising research directions for antenna arrays,” _Digital Signal Processing_ , vol. 94, pp. 3 – 20, 2019, special Issue on Source Localization in Massive MIMO.
  * [2] D. J. Love, R. W. Heath, V. K. N. Lau, D. Gesbert, B. D. Rao, and M. Andrews, “An overview of limited feedback in wireless communication systems,” _IEEE J. Sel. Areas Commun._ , vol. 26, no. 8, pp. 1341–1365, 2008.
  * [3] W. Utschick, V. Rizzello, M. Joham, Z. Ma, and L. Piazzi, “Learning the csi recovery in fdd systems,” _IEEE Trans. on Wireless Commun._ , 2022.
  * [4] B. Fesl, N. Turan, M. Koller, M. Joham, and W. Utschick, “Centralized learning of the distributed downlink channel estimators in FDD systems using uplink data,” in _WSA 2021; 25th Int. ITG Workshop on Smart Antennas_ , 2021, pp. 1–6.
  * [5] N. Turan, M. Koller, S. Bazzi, W. Xu, and W. Utschick, “Unsupervised learning of adaptive codebooks for deep feedback encoding in fdd systems,” in _2021 55th Asilomar Conference on Signals, Systems, and Computers_ , 2021, pp. 1464–1469.
  * [6] N. Turan, M. Koller, V. Rizzello, B. Fesl, S. Bazzi, W. Xu, and W. Utschick, “On distributional invariances between downlink and uplink MIMO channels,” in _WSA 2021; 25th Int. ITG Workshop on Smart Antennas_ , 2021, pp. 1–6.
  * [7] T. T. Nguyen, H. D. Nguyen, F. Chamroukhi, and G. J. McLachlan, “Approximation by finite mixtures of continuous density functions that vanish at infinity,” _Cogent Math. Statist._ , vol. 7, no. 1, p. 1750861, 2020.
  * [8] Y. Linde, A. Buzo, and R. Gray, “An algorithm for vector quantizer design,” _IEEE Trans. Commun._ , vol. 28, no. 1, pp. 84–95, Jan. 1980.
  * [9] V. Lau, Y. Liu, and T.-A. Chen, “On the design of MIMO block-fading channels with feedback-link capacity constraint,” _IEEE Trans. on Commun._ , vol. 52, no. 1, pp. 62–70, 2004.
  * [10] A. Goldsmith, _Wireless Communications_. Cambridge Univ. Press, 2005.
  * [11] A. Goldsmith, S. Jafar, N. Jindal, and S. Vishwanath, “Capacity limits of MIMO channels,” _IEEE J. Sel. Areas Commun._ , vol. 21, no. 5, pp. 684–702, 2003.
  * [12] I. E. Telatar, “Capacity of multi-antenna gaussian channels,” _Eur. Trans. on Telecommun._ , vol. 10, pp. 585–595, 1999.
  * [13] S. Jaeckel, L. Raschkowski, K. Börner, and L. Thiele, “Quadriga: A 3-d multi-cell channel model with time evolution for enabling virtual field trials,” _IEEE Trans. Antennas Propag._ , vol. 62, no. 6, pp. 3242–3256, 2014.
  * [14] S. Jaeckel, L. Raschkowski, K. Börner, L. Thiele, F. Burkhardt, and E. Eberlein, “Quadriga: Quasi deterministic radio channel generator, user manual and documentation,” Fraunhofer Heinrich Hertz Institute, Tech. Rep., v2.2.0, 2019.
  * [15] Y. Tsai, L. Zheng, and X. Wang, “Millimeter-wave beamformed full-dimensional MIMO channel estimation based on atomic norm minimization,” _IEEE Trans. on Commun._ , vol. 66, no. 12, pp. 6150–6163, 2018.
  * [16] C. M. Bishop, _Pattern Recognition and Machine Learning (Information Science and Statistics)_. Berlin, Heidelberg: Springer-Verlag, 2006.
  * [17] M. Koller, B. Fesl, N. Turan, and W. Utschick, “An Asymptotically Optimal Approximation of the Conditional Mean Channel Estimator based on Gaussian Mixture Models,” in _IEEE Int. Conf. on Acoust., Speech and Signal Process._ , 2022, pp. 5268–5272.
  * [18] ——, “An Asymptotically Optimal Approximation of the Conditional Mean Channel Estimator based on Gaussian Mixture Models,” 2021, arXiv preprint: 2112.12499.
  * [19] A. Alkhateeb, G. Leus, and R. W. Heath, “Compressed sensing based multi-user millimeter wave systems: How many measurements are needed?” in _2015 IEEE Int. Conf. on Acoust., Speech and Signal Process. (ICASSP)_ , 2015, pp. 2909–2913.
  * [20] M. Gharavi-Alkhansari and T. Huang, “A fast orthogonal matching pursuit algorithm,” in _Proceedings of the 1998 IEEE Int. Conf. on Acoust., Speech and Signal Process., ICASSP ’98 (Cat. No.98CH36181)_ , vol. 3, 1998, pp. 1389–1392.
  * [21] R. Hunger, D. A. Schmidt, M. Joham, and W. Utschick, “A general covariance-based optimization framework using orthogonal projections,” in _2008 IEEE 9th Workshop on Signal Process. Advances in Wireless Commun._ , 2008, pp. 76–80.
  * [22] J. Kermoal, L. Schumacher, K. Pedersen, P. Mogensen, and F. Frederiksen, “A stochastic MIMO radio channel model with experimental validation,” _IEEE J. Sel. Areas Commun._ , vol. 20, no. 6, pp. 1211–1226, 2002.