# Exact quantum revivals for the Dirac equation

Fernando Chamizo  Osvaldo P. Santillán

###### Abstract

In the present work, the results obtained in [25] about the revivals of a
relativistic fermion wave function on a torus are considerably enlarged. In
fact, all the possible quantum states exhibiting revivals are fully
characterized. The revivals are exact, that is, are true revivals without
taking any particular limit such as the non relativistic one. The present
results are of interest since they generalize the Talbot effect and the
revivals already found by studying the Schrödinger equation to a relativistic
situation with nonzero mass. This makes the problem nontrivial, as the
dispersion relation is modified and is not linear. The results are obtained by
the use of arithmetic tools which are described in certain detail along the
text.

##  1 Introduction

In quantum mechanics, the solutions of the Schrödinger equation composed by
waves of different wavelengths posses varying phase velocities. It is said
that the equation is dispersive and it implies that when time evolves the
solution is expected to randomize or spread out. This can be considered as a
manifestation of the uncertainty principle. A basic example is the evolution
of an initial Gaussian packet, which at a time $t$ is given by

$\Psi(\mathbf{r},t)=\frac{\exp\big{(}-\frac{1}{2}(a+i\hbar
t/m)^{-1}\mathbf{r}\cdot\mathbf{r}\big{)}}{(\pi a)^{3/4}(1+i\hbar
t/ma)^{3/2}}.$

Despite this dispersive behaviour, in some situations there is a coherence
between the phase velocities causing _quantum revivals_. This is an exact or
approximate periodicity in time. For instance, if the energy spectrum is
discrete, say $\\{E_{n}\\}_{n=1}^{\infty}$, and there exists
$s\in{\mathbb{R}}^{+}$ such that
$\\{sE_{n}\\}_{n=1}^{\infty}\subset{\mathbb{Z}}$ then the solutions of the
Schrödinger equation are of the form

$\Psi(\mathbf{r},t)=\sum_{n=1}^{\infty}a_{n}\psi_{n}(\mathbf{r})\,e^{-iE_{n}t/\hbar},$

and an exact revival $\Psi(\mathbf{r},t)=\Psi(\mathbf{r},t+T_{\text{\rm
rev}})$ appears with $T_{\text{\rm rev}}=2\pi s\hbar$ the _revival time_. A
simple example of this situation is the infinite quantum well. The changes
suffered by the initial state at a fractional multiple of the revival time and
its fractalization for irrational multiples have been studied for several
quantum systems in many works [1], [22], [23], [12], [15], [7], establishing
some interesting connections with the classic Talbot effect in optics and some
topics in number theory [3].

Some simulations related to quantum particles in boxes [20], [21] may suggest
that taking into account relativistic effects, the coherence of the phase
velocities is lost and one can only expect approximate quantum revivals (as it
occurs in many models cf. [26]). However, in [25] it is shown a different
truth. In this reference, it is introduced a model associated to a fermionic
particle confined to a circle of radius $R$ under Dirac equation and it is
shown that if the dimensionless quantity $q=McR/\hbar$, where $M$ is the mass,
is a positive even integer one can construct a finite number of plane wave
spinor solutions such that a linear combination of them presents an exact
quantum revival. An elegant part of the argument is that these plane waves are
associated to Pythagorean triples. These are integral solutions of the
algebraic equation $x^{2}+y^{2}=z^{2}$ [16, §13.2]) which constitute a very
old topic in pure mathematics. This observation introduces a new connection
between quantum revivals and number theory, to be explored further in the
present paper.

With the methods developed in [25], it is seen and it will be explained along
the present text that, for each $q\in 2{\mathbb{Z}}^{+}$ there are only few
possible plane wave solutions. These solutions are related to the divisors of
$q$ (on average this number goes like $\log q$ [16, §18.2]). This fact
suggests that there is little freedom to construct examples following those
lines. In this paper this model is analyzed with more advanced arithmetic
tools. In these terms, it turns out that for each $q$ satisfying a less
restrictive condition, apart from some trivial cases, there are infinitely
many plane wave solutions such that any linear combination of them presents
exact revivals. In fact, for each valid $q$ all the states showing these
revivals are characterized in an algorithmic way. The exponential growth of
the energies explains some fractal like _quantum carpets_ 111Which are density
plots of the state along a period [4]. and connects with some topics in
fractal geometry. The case of particles confined in a square flat torus is
also studied. The number of plane waves serving as building blocks to obtain
the states possessing exact quantum revivals is still infinite for suitable
fixed values of $q$ in a dense set, but in this 2-dimensional case the
possible energies present a mild growth allowing high degeneracies and more
freedom to choose the initial conditions. Here it appears another new
connection with arithmetical topics through $\text{\rm SL}_{2}({\mathbb{Z}})$
and some allied discrete groups.

The present work is organized as follows. In section 2 the defining equation
and the representation employed along the text are clarified. Also, the main
problem of characterization is stated. In section 3 the possible revivals in a
one dimensional torus are characterized in terms of solutions of Pell’s
equations, and the presence of an infinite number of states exhibiting quantum
revivals is pointed out. In section 4 the two dimensional case is analyzed.
This situation is more complex and its analysis relies on arithmetic topics
related to Fuchsian groups. Section 5 contains the discussion of the presented
results.

##  2 The basic model

The plane wave solutions of the Dirac equation under the Dirac representation
are given by the spinors [18]

$\Psi=\begin{pmatrix}(k^{0}+Mc)\varphi\\\
(\vec{\sigma}\cdot\vec{k})\varphi\end{pmatrix}e^{-ik\cdot
x/\hbar}\qquad{and}\qquad\widetilde{\Psi}=-\begin{pmatrix}(\vec{\sigma}\cdot\vec{k})\eta\\\
(k^{0}+Mc)\eta\end{pmatrix}e^{ik\cdot x/\hbar}$

where $\varphi,\eta\in{\mathbb{C}}^{2}$ with $|\varphi|=|\eta|=1$ and $M$ is
the mass. At this point, these spinors are not normalized. In the following,
the attention will focused on $\Psi$ since $\widetilde{\Psi}$ corresponds to
antiparticles, and could be treated in a symmetric way.

Consider particles confined in a flat torus obtained identifying the opposite
sides of

$[2\pi R_{1},0]\times[0,2\pi R_{2}]\times\\{0\\}.$

In this situation the cartesian coordinates may be parameterized as
$x=R_{1}\phi_{1}$, $y=R_{2}\phi_{2}$, $z=0$ with $\phi_{1}$ and $\phi_{2}$ two
$2\pi$-periodic angles. The formula for $\Psi$ imposes the following
quantization of the momentum

$\vec{k}=\Big{(}\frac{\hbar n_{1}}{R_{1}},\frac{\hbar
n_{2}}{R_{2}},0\Big{)}\qquad\text{with}\quad\vec{n}=(n_{1},n_{2})\in{\mathbb{Z}}^{2}.$

Let $\omega_{\vec{n}}=\hbar^{-1}E_{\vec{n}}>0$ be the angular frequency for
each $\vec{n}$ and consider the complex quatity

$z_{\vec{n}}=\frac{\hbar n_{1}}{R_{1}}+\frac{i\hbar n_{2}}{R_{2}}$

associated to the momentum. The generic energy equation
$E^{2}=M^{2}c^{4}+c^{2}\vec{k}^{2}$ implies

$\omega_{\vec{n}}=\frac{c}{\hbar}\sqrt{M^{2}c^{2}+|z_{\vec{n}}|^{2}}$

which in particular gives the following useful relation

$\Big{(}\frac{\hbar\omega_{\vec{n}}}{c}+Mc\Big{)}^{2}+|z_{\vec{n}}|^{2}=\frac{2\hbar\omega_{\vec{n}}}{c}\Big{(}\frac{\hbar\omega_{\vec{n}}}{c}+Mc\Big{)}.$

The state $\Psi$ is not normalized. The normalization (up to a $1/2\pi$
factor) is achieved by dividing by $\sqrt{(\frac{\hbar
E}{c})^{2}+|z_{\vec{n}}|^{2}}$, this leads to the following normalized spinor

$\Psi_{\vec{n}}(\phi_{1},\phi_{2},t)=\sqrt{\frac{\hbar\omega_{\vec{n}}+Mc^{2}}{2\hbar\omega_{\vec{n}}}}\begin{pmatrix}\varphi\\\\[5.0pt]
\frac{cz_{\vec{n}}^{*}}{\hbar\omega_{\vec{n}}+Mc^{2}}\varphi_{2}\\\\[5.0pt]
\frac{cz_{\vec{n}}}{\hbar\omega_{\vec{n}}+Mc^{2}}\varphi_{1}\end{pmatrix}e^{i(n_{1}\phi_{1}+n_{2}\phi_{2}-\omega_{\vec{n}}t)}.$

By use of appropriate spinorial transformation laws, with a suitable rotation
of the reference frame $\varphi$ is mapped into $(1,0)^{\dagger}$. Assume in
the sequel $\varphi=(1,0)^{\dagger}$ for the sake of simplicity, as this
choice allows to avoid spurious multiplicities in the solutions. Different
choices of $\varphi$ only affect to the amplitude of the components of the
Dirac spinor while the frequency, which is the important issue for quantum
revivals, remains invariant.

Some particular limits are in order. The 1-dimensional case corresponds to let
$R_{1}\to\infty$, so the periodicity in the $x$ coordinate is killed,
$k^{1}=0$ and the particle is confined to move in a (flat) circle of radius
$R=R_{2}$ given by a $2\pi R$-periodic $y$ coordinate. The 2-dimensional case
corresponds to the model above choosing $R=R_{1}=R_{2}$ (a square flat torus).
Therefore, $z_{\vec{n}}=i\hbar n_{2}/R$ in the first case and
$z_{\vec{n}}=\hbar(n_{1}+in_{2})/R$ in the second one.

Both cases can be written in a unified way by introducing the adimensional
positive quantity

$q=\frac{McR}{\hbar},$

from where

$E_{\vec{n}}=\frac{c\hbar}{R}\sqrt{\vec{n}^{2}+q^{2}},\qquad\omega_{\vec{n}}=\frac{c}{R}\sqrt{\vec{n}^{2}+q^{2}},$

and the formula

$\Psi_{\vec{n}}=\frac{1}{\sqrt{2}}\bigg{(}1+\frac{q}{\sqrt{\vec{n}^{2}+q^{2}}}\bigg{)}^{1/2}\begin{pmatrix}1\\\
0\\\ 0\\\
\frac{n_{1}+in_{2}}{q+\sqrt{\vec{n}^{2}+q^{2}}}\end{pmatrix}e^{i(n_{1}\phi_{1}+n_{2}\phi_{2}-\sqrt{\vec{n}^{2}+q^{2}}\,ct/R)},$
(2.1)

follows. The 1-dimensional case corresponds to set $n_{1}=0$ (so, the state
does not depend on $\phi_{1}$) and the $2$-dimensional case to leave it free.
We took apart a $1/2\pi$ factor from the normalization to avoid conflicts with
the choice of the dimension.

The spinor given above is an energy eigenfunction. A generic state with
$\varphi=(1,0)^{\dagger}$, having a discrete, finite and nonnegative energy
spectrum is a superposition of the form

$\Phi=\sum_{\vec{n}\in\mathcal{N}}\lambda_{\vec{n}}\Psi_{\vec{n}}\qquad\text{with}\quad\lambda_{\vec{n}}\in{\mathbb{C}}-\\{0\\}$
(2.2)

where
$\mathcal{N}=\big{\\{}\vec{n}_{0},\vec{n}_{1},\dots,\vec{n}_{N}\big{\\}}$,
$\vec{n}_{j}=(k_{j},l_{j})\in{\mathbb{Z}}^{2}$ and $k_{j}=0$ in the one
dimensional model. The energy spectrum is related to the squared norm of the
$\vec{n}_{j}$. It is the set

$\mathcal{E}=\Big{\\{}\frac{c\hbar}{R}\sqrt{k_{j}^{2}+\ell_{j}^{2}+q^{2}}\,:\,(k_{j},\ell_{j})\in\mathcal{N}\Big{\\}}.$

Note that cardinality $\\#\mathcal{E}$ of this set may be different from $N$,
as there is a possible degeneracy of states with the same energy value. At
this point, it is important to state the following affirmation, which is a
direct and simple consequence of the spinor formulas given above.

> _Affirmation:_ The state described by formulas (2.1)-(2.2) is periodic in
> time, and thus exhibits quantum revival, if and only if for each $0\leq
> j\leq N$ there exists an irreducible fraction $a_{j}/b_{j}>0$ such that
>
>
> $k_{j}^{2}+\ell_{j}^{2}+q^{2}=\frac{a_{j}^{2}}{b_{j}^{2}}\big{(}k_{0}^{2}+\ell_{0}^{2}+q^{2}\big{)}.$
> (2.3)
>
> Here $N$ is the number of elements of the state (2.2). The revival time is
> given by the formula
>
> $T_{\rm rev}=\frac{2\pi L}{\omega_{\vec{n}_{0}}}\qquad\text{with}\quad
> L=\text{\rm lcm}\,\big{(}\\{b_{j}\\}_{j=0}^{N}\big{)}$
>
> where, as usual, lcm denotes the least common multiple.

This affirmation is clear from the formulas just stated, and it is plain that
the time periodicity in (2.2) is equivalent to the commensurability of the
energies in $\mathcal{E}$. However, the solutions (2.3) are not necessarily
trivial from the arithmetic point of view. In fact, as it will be shown below,
they are described in terms of generalized Pell’s equations in the one
dimensional case and they relate to some Fuchsian groups in the 2D case.

The philosophy that will guide the present paper is perhaps different than the
standard one. Instead of looking for particular states which exhibit quantum
revivals, the goal is to characterize all of them and it is achieved fixing an
specific vector state $\vec{n}_{0}$ and considering all possible $\Phi$ in
(2.2) containing $\Psi_{\vec{n}_{0}}$.

> _Philosophy:_ To characterize the possible sets $\mathcal{N}$ containing
> $\vec{n}_{0}$ such that quantum revivals appear. More concretely, given
> $\vec{n}_{0}$ we want to determine the maximal set of quantum numbers
> $\mathcal{N}_{0}$ including $\vec{n}_{0}$,
>
> $\boxed{\Phi\text{ is periodic in
> }t\quad\Leftrightarrow\quad\mathcal{N}\text{ is a finite subset of
> }\mathcal{N}_{0}.}$ (2.4)
>
> It is also addressed the computation of the period for a given
> $\mathcal{N}$.

In principle, $\mathcal{N}_{0}$ may be finite or infinite although it will be
shown that it is infinite if somewhat trivial revivals are excluded. Let us
elucidate this pointing out some remarks and advancing some results. The
reader that finds the following explanation too technical may skip this
paragraph and read it later on. First, note that if $\\#\mathcal{E}>1$ i.e.,
if $\Phi$ in (2.2) has not a definite energy, then there is some
$a_{j}/b_{j}\neq 1$ in (2.3) and it implies $q^{2}\in{\mathbb{Q}}^{+}$.
Assuming this, it will be shown in the subsequent sections that
$\mathcal{N}_{0}$ is infinite except in the one dimensional case when
$\vec{n}_{0}^{2}+q^{2}$ is the square of an irreducible fraction, which leads
to rather obvious divisibility conditions for the numerator and denominator of
$q^{2}$ that allow to construct $\mathcal{N}_{0}$.

If $\\#\mathcal{E}=1$ there is only an energy, thus $a_{j}/b_{j}=1$ and no
restriction is needed on $q$. Assume $q^{2}\not\in{\mathbb{Q}}$ because
otherwise the previous considerations apply. In the one dimensional situation
$\Psi_{(0,\ell_{0})}$ can be only accompanied with $\Psi_{(0,-\ell_{0})}$
giving obvious revivals in (2.2) and $\mathcal{N}_{0}=2$. In the 2D case the
situation is less trivial because $\\#\mathcal{N}$ can be arbitrarily large
keeping $\\#\mathcal{E}=1$ because the function giving the number of
representations as a sum of two squares is unbounded222Although it is $\pi$ on
average, as it is easily checked approximating the number of integer points in
a circle by its area [16, Th. 339].. So, $\mathcal{N}_{0}$ is finite but not
uniformly bounded in terms of $\vec{n}_{0}$. For instance,
$\vec{n}_{0}=(178,19)$ can be complemented with other $63$ integral pairs to
get the same
$R\omega_{\vec{n}_{j}}/c=\sqrt{\vec{n}_{0}^{2}+q^{2}}=\sqrt{32045+q^{2}}$

$178^{2}+19^{2}=166^{2}+67^{2}=157^{2}+86^{2}=179^{2}+2^{2}=142^{2}+109^{2}=163^{2}+74^{2}=\dots$

This phenomenon of unbounded degeneracy of the energy leads typically to
intricate density probabilities for a fixed time (see the section of numerical
results) and it is linked to some unsolved problems about the structure of the
nodal lines [5], [2].

##  3 The one dimensional case

Recall that in the one dimensional case the first coordinate of $\vec{n}$ is
set to $0$, and the necessary and sufficient condition (2.3) for having
quantum revivals becomes

$\ell_{j}^{2}+q^{2}=\frac{a_{j}^{2}}{b_{j}^{2}}\big{(}\ell_{0}^{2}+q^{2}\big{)}\qquad\text{for}\quad
0\leq j\leq N.$ (3.1)

If $q^{2}\not\in{\mathbb{Q}}^{+}$ then $a_{j}/b_{j}=1$ and $\\#\mathcal{E}=1$
leading to a trivial case, as shown in the previous section. Then it can be
assumed $q^{2}\in{\mathbb{Q}}^{+}$. In fact, it is convenient to focus on
$q^{2}\in{\mathbb{Z}}^{+}$ because it is somewhat simpler and contains all the
ingredients appearing in the general situation. Some concepts of arithmetic
are in order. A positive integer is said to be _squarefree_ if it is not
divisible by a square different from $1^{2}$. Each positive integer admits a
unique factorization as a product of a squarefree and a square. Consider such
decomposition for $\ell_{0}^{2}+q^{2}$:

$\ell_{0}^{2}+q^{2}=Ds^{2}\qquad\text{with $D$ squarefree}.$ (3.2)

Then (3.1) can be arranged as

$\ell_{j}-D\Big{(}\frac{a_{j}s}{b_{j}}\Big{)}^{2}=-q^{2}.$ (3.3)

As $D$ is squarefree and $q^{2}$ is integer, $a_{j}s/b_{j}$ is integer too.
So, all the possibilities for $\ell_{j}$ and the corresponding $a_{j}/b_{j}$
are in one to one correspondence with the integer solutions of the
_generalized Pell’s equation_

$x^{2}-Dy^{2}=-q^{2}.$ (3.4)

In these terms, the main affirmation of the previous section, adapted to the
one dimensional case, is the following.

> _Affirmation:_ Given $\ell_{0}\in{\mathbb{Z}}$ and
> $q^{2}\in{\mathbb{Z}}^{+}$, defining $D$ and $s$ as in (3.2), the set
> $\mathcal{N}_{0}$ of quantum numbers to have (2.4) when
> $(0,\ell_{0})\in\mathcal{N}$ is
>
> $\mathcal{N}_{0}=\big{\\{}(0,x)\,:\,x^{2}-Dy^{2}=-q^{2}\text{ with
> }(x,y)\in{\mathbb{Z}}^{2}\big{\\}}.$
>
> Moreover the revival time is given by
>
> $T_{\rm rev}=\frac{2\pi RL}{cs\sqrt{D}}=\frac{2\pi
> L}{\omega_{\vec{n}_{0}}},$
>
> where $L$ is the lcm of the denominators of
> $a_{j}/b_{j}=s^{-1}\sqrt{(\ell_{j}^{2}+q^{2})/D}$ for the
> $(0,\ell_{j})\in\mathcal{N}$.

In view of the last affirmation, it is mandatory to review the solutions of
the generalized Pell’s equation (3.4) that allows to describe the structure of
$\mathcal{N}_{0}$ for any choice of $\ell_{0}\in{\mathbb{Z}}$ and
$q^{2}\in{\mathbb{Z}}^{+}$.

The case $D=1$ is somewhat trivial because $x^{2}-Dy^{2}=(x-y)(x+y)=-q^{2}$
implies that $x-y$ and $x+y$ are divisors of $q^{2}$ leaving only a finite
number of possibilities. In other words, if $\ell_{0}^{2}+q^{2}$ is a square,
$\mathcal{N}_{0}$ is finite. This includes the case covered in [25] using
Pythagorean triples. By introducing the quantity $\mp d=x-y$ and, consequently
$\pm q^{2}/d=x+y$, it is obtained the explicit description

$\mathcal{N}_{0}=\Big{\\{}\Big{(}0,\frac{q^{2}-d^{2}}{2d}\Big{)}\,:\,d\mid
q^{2}\text{ with }q^{2}/d-d\text{ even}\Big{\\}}.$

The parity condition is to assure $(q^{2}-d^{2})/2d$ to be integer. It is
superfluous if $q^{2}$ is odd. Note that the $\pm$ is unnecessary by the anti-
symmetry $d\leftrightarrow q^{2}/d$.

If $D\neq 1$, $\mathcal{N}_{0}$ is instead an infinite set, which can be seen
as follows. For the standard _Pell’s equation_

$x^{2}-Dy^{2}=1,$

a classic result due to Lagrange asserts that if $D\in{\mathbb{Z}}^{+}$ is not
a square (this is assured in our case because $D$ is squarefree and $D\neq 1$)
then this equation always has infinitely many solutions
$(x,y)\in{\mathbb{Z}}^{+}\times{\mathbb{Z}}^{+}$. Lagrange’s theorem also
characterizes these solutions. They can be expressed in the synthetic way as
follows

$x+y\sqrt{D}=\big{(}x_{p}+y_{p}\sqrt{D}\big{)}^{n}\qquad\text{with}\quad
n\in{\mathbb{Z}}^{+}$

where $(x_{p},y_{p})$ is the minimal positive solution. Namely, the one having
$x_{p}$ and hence $y_{p}$ as small as possible. It is easy to check that the
previous formula gives solutions by simply multiplying by the conjugate
because

$\big{(}x_{p}+y_{p}\sqrt{D}\big{)}\big{(}x_{p}-y_{p}\sqrt{D}\big{)}=x_{p}^{2}-Dy_{p}^{2}=1.$

In the present case, it is deduced from (3.2) that $(\ell_{0},s)$ is an
integer solution of (3.4). The same argument just employed, using the
conjugate, shows that

$x+y\sqrt{D}=\pm\big{(}\ell_{0}+s\sqrt{D}\big{)}\big{(}x_{p}+y_{p}\sqrt{D}\big{)}^{n}\qquad\text{with}\quad
n\in{\mathbb{Z}}$ (3.5)

are infinitely many integer solutions of (3.4), proving that $\mathcal{N}_{0}$
is infinite. The $\pm$ and the extension of the range of $n$ from
${\mathbb{Z}}^{+}$ to ${\mathbb{Z}}$ are only to take into consideration
different combinations of signs. Note that
$x_{p}-y_{p}\sqrt{D}=\big{(}x_{p}+y_{p}\sqrt{D}\big{)}^{-1}$.

Let us illustrate the situation with an example. If $\ell_{0}=5$ and
$q=\sqrt{2}$ then (3.2) produces $D=3$ and $s=3$. A direct search shows that
$(x_{p},y_{p})=(2,1)$ is the minimal positive solution of $x^{2}-3y^{2}=1$.
The values $n=-1,0,1,2$ in (3.5) give

$x+y\sqrt{3}=\pm(1+\sqrt{3}),\,\pm(5+3\sqrt{3}),\,\pm(19+11\sqrt{3}),\,\pm(71+41\sqrt{3}).$

Hence $(0,\pm 1)$, $(0,-5)$, $(0,\pm 19)$ and $(0,\pm 71)$ are in the
$\mathcal{N}_{0}$ corresponding to $(0,\ell_{0})=(0,5)$.

The formula (3.5) establishes a linear recurrence of the first order between
the solutions corresponding to $n$ and $n+1$ (cf. [10]). Eliminating a
variable, it is obtained a second order recurrence. A complication of the
theory (related to deep topics in number theory as unique factorization and
the class number [8]) is that a finite number of similar families of solutions
may exists. Namely, in our case they may appear families of integer solutions
replacing $\ell_{0}$ and $s$ by $a$ and $b$ with $(x,y)=(a,b)$ satisfying
(3.4) but it can be proved that $0<a<q(1+\sqrt{2x_{p}})/2$ [9] then there is
only a finite number of possibilities to explore. All the possible families
obey to the same recurrence law because they always come from the
multiplication by $x_{p}+y_{p}\sqrt{D}$.

In the example given above, $(x_{p},y_{p})=(2,1)$ and the only valid pairs
$(a,b)$ are $(1,\pm 1)$ which belong to the same family as
$(\ell_{0},s)=(5,3)$ because $1+\sqrt{3}=(5+3\sqrt{3})(2+\sqrt{3})^{-1}$ and
$1-\sqrt{3}=-(5+3\sqrt{3})(2+\sqrt{3})^{-2}$.

Note that always (3.5) can be translated in one sided recurrence formulas
dividing into $n\geq n_{0}$ and $n<n_{0}$. With this remark and the previous
considerations, it is deduced

> There exists a finite number of sequences $c_{n}^{(1)}$,…, $c_{n}^{(J)}$
> satisfying a second order linear recurrence such that
>
> $\mathcal{N}_{0}=\big{\\{}(0,\pm c_{n}^{(j)})\text{ for {$1\leq j\leq J$}
> and $n\geq 0$}\big{\\}}.$

For instance, in the previous example it can be proved (see the section of
numerical results)

$\mathcal{N}_{0}=\big{\\{}(0,\pm c_{n})\text{ where
}c_{n+2}=4c_{n+1}-c_{n}\text{ with }c_{0}=1,\ c_{1}=5\big{\\}}.$ (3.6)

The sequence has an exponential growth, due to (3.5),

$\\{c_{n}\\}=\\{1,5,19,71,265,989,3691,13775,51409,191861,716035,2672279,9973081,\dots\\}.$

Depending on the elements chosen to compose $\mathcal{N}_{0}$, the revival
time can be at most $6\pi/\omega_{\vec{n}_{0}}$, because
$a_{j}s/b_{j}=3a_{j}/b_{j}$ must be an integer and hence $L$ divides $3$.

The previous paragraphs constitute a quick review of the theory. For the
reader interested in deepening in the theory and practice of the generalized
Pell’s equation, some interesting references are [24], [9], [8] and [17]. It
may be remarked that the classical continued fractions (nowadays mostly
excluded from the standard undergraduate syllabi) are of great help for the
efficient computation of the families composing the solution.

The next task is to mention how the previous ideas can be modified to cover
the remaining cases $q^{2}\in{\mathbb{Q}}^{+}-{\mathbb{Z}}$ preserving the
same result. This will be done briefly. In this situation, (3.2) is
substituted by the squarefree-square decomposition of the numerator and the
denominator of $\ell_{0}^{2}+q^{2}$ to get

$\ell_{0}^{2}+q^{2}=\frac{Ds^{2}}{D_{*}s_{*}^{2}}\qquad\text{with $D$ and
$D_{*}$ squarefree}.$

Clearing denominators in (3.3), it is obtained that $(\ell_{j},a_{j}s/b_{j})$
is a solution of

$D_{*}s_{*}^{2}x^{2}-Dy^{2}=-q^{2}D_{*}s_{*}^{2}.$

In fact it is an integral solution because
$-q^{2}D_{*}s_{*}^{2}\in{\mathbb{Z}}$ and the tentative denominator of
$(a_{j}s/b_{j})^{2}$ could not be canceled with $D$ because the latter is
squarefree. So, to cover the full rational case in the affirmation above,
$\mathcal{N}_{0}$ must be generalized to

$\mathcal{N}_{0}=\big{\\{}(0,x)\,:\,D_{*}s_{*}^{2}x^{2}-Dy^{2}=-q^{2}D_{*}s_{*}^{2}\text{
with }(x,y)\in{\mathbb{Z}}^{2}\big{\\}}.$ (3.7)

Note that this is actually a generalization because $D_{*}s_{*}^{2}=1$ if
$q^{2}\in{\mathbb{Z}}$.

The case $D=D_{*}=1$ leads, as before, to a finite number of possibilities
because $s_{*}x+y$ and $s_{*}x-y$ are divisors of $-q^{2}s_{*}^{2}$. The
explicit description of $\mathcal{N}_{0}$ is

$\mathcal{N}_{0}=\Big{\\{}\Big{(}0,\frac{q^{2}s_{*}^{2}-d^{2}}{2s_{*}d}\Big{)}\,:\,d\mid
q^{2}s_{*}^{2}\text{ and }q^{2}s_{*}^{2}/d-d\text{ multiple of
}2s_{*}\Big{\\}}.$

When $s_{*}$ grows, being a multiple of $2s_{*}$ imposes a strong condition
and $\mathcal{N}_{0}$ might reduce to the trivial set $\\{(0,\pm\ell_{0})\\}$
even if $q^{2}s_{*}^{2}$ has nontrivial divisors. An example of this situation
is $\ell_{0}=2$, $q=\frac{2}{5}\sqrt{299}$.

In the rest of the cases ($D$ and $D_{*}$ not simultaneously $1$), a
modification of the arguments for $q^{2}\in{\mathbb{Z}}$ applies. Recall that
$(\ell_{0},s)$ is a solution of the equation in (3.7). Reasoning as in (3.5)
multiplying by the conjugate, it is deduced that

$xs_{*}\sqrt{D_{*}}+ys\sqrt{D}=\pm\big{(}\ell_{0}s_{*}\sqrt{D_{*}}+s\sqrt{D}\big{)}\big{(}x_{p}+y_{p}ss_{*}\sqrt{DD_{*}}\big{)}^{n}\qquad\text{with}\quad
n\in{\mathbb{Z}}$

gives infinitely many integer solutions where $(x_{p},y_{p})$ is the minimal
positive solution of the Pell equation

$x^{2}-s^{2}s_{*}^{2}DD_{*}y^{2}=1.$

Its existence is assured because, once $D=D_{*}=1$ is excluded,
$s^{2}s_{*}^{2}DD_{*}$ cannot be a square. Again, a finite number of other
families may appear.

Summing up, for $q^{2}\in{\mathbb{Q}}$ the set $\mathcal{N}_{0}$ is infinite
except in the case in which $\ell_{0}^{2}+q^{2}$ is the square of an
irreducible fraction.

##  4 The two dimensional case

The two dimensional case is much more involved because, as it will be shown
below, the families of eigenstates presenting exact quantum revivals are
parametrized Fuchsian groups i.e., discrete subgroups of $\text{\rm
PSL}_{2}({\mathbb{R}})$. Recall that (2.3) implies $q^{2}\in{\mathbb{Q}}$
except in the case $\\#\mathcal{E}=1$ in which $a_{j}/b_{j}=1$. As in the one
dimensional case, the focus will be on the case $q^{2}\in{\mathbb{Z}}$ and
some closing comments about the fractional case will be made later on. Some
considerations about $\\#\mathcal{E}=1$ were included above and it will be
expanded in the next section.

Given $q^{2}\in{\mathbb{Z}}$, consider the decomposition

$k_{0}^{2}+\ell_{0}^{2}+q^{2}=Ds^{2}\qquad\text{with $D$ squarefree}.$ (4.1)

The condition (2.3) reads

$k_{j}^{2}+\ell_{j}^{2}+q^{2}=D(a_{j}s/b_{j})^{2},$

where $a_{j}s/b_{j}$ must be integer333This is because $D$ is squarefree and
$k_{j}^{2},\ell_{j}^{2},q^{2}\in{\mathbb{Z}}$. and it follows a natural
analogue of the one dimensional result.

> _Affirmation:_ Given $(k_{0},\ell_{0})\in{\mathbb{Z}}^{2}$ and
> $q^{2}\in{\mathbb{Z}}^{+}$, then by defining $D$ and $s$ as in (4.1), the
> set $\mathcal{N}_{0}$ of quantum numbers needed for (2.4) with
> $(k_{0},\ell_{0})\in\mathcal{N}$ to hold is
>
> $\mathcal{N}_{0}=\big{\\{}(y,z)\,:\,Dx^{2}-y^{2}-z^{2}=q^{2}\text{ with
> }(x,y,z)\in{\mathbb{Z}}^{3}\big{\\}}.$
>
> Moreover the revival time is given by
>
> $T_{\rm rev}=\frac{2\pi RL}{cs\sqrt{D}}=\frac{2\pi L}{\omega_{\vec{n}_{0}}}$
>
> where $L$ is the lcm of the denominators of
> $s^{-1}\sqrt{(k_{j}^{2}+\ell_{j}^{2}+q^{2})/D}$ for the
> $(k_{j},\ell_{j})\in\mathcal{N}$.

This statement suggests that the study of the group that leaves invariant the
quadratic form $Dx^{2}-y^{2}-z^{2}$ and which applies integer solutions into
integer solutions is of importance, as it leads to an algorithm generating new
solutions. The resulting group is related to the Lorentz group, but it is
rather clear that it does not contain all the elements. The task is to
determine which elements are to be included. This group is not present in the
one dimensional case.

It may be convenient to explain the difference between the one and two
dimensional case from an abstract point of view. As reviewed above, the
structure of the solutions in the one dimensional case can be summarized
saying that there were finitely many families of solutions and each family was
given by the action of the discrete group
$\big{\\{}\big{(}x_{p}+y_{p}\sqrt{D}\big{)}^{n}\big{\\}}_{n\in{\mathbb{Z}}}$,
which is isomorphic to ${\mathbb{Z}}$, producing sparse solutions because of
the exponential growth. In the two dimensional case instead, the group is non-
abelian and gives a denser set of solutions.

In order to study the above mentioned group, let us start with the case $D=1$,
which is no longer trivial. Following [11], we note that the map

$\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\longmapsto\begin{pmatrix}(a^{2}+b^{2}+c^{2}+d^{2})/2&(-a^{2}+b^{2}-c^{2}+d^{2})/2&-ab-
cd\\\ (-a^{2}-b^{2}+c^{2}+d^{2})/2&(a^{2}-b^{2}-c^{2}+d^{2})/2&ab-cd\\\ -ac-
bd&ac-bd&ad+bc\end{pmatrix}$ (4.2)

establishes an isomorphism between $\text{\rm PSL}_{2}({\mathbb{R}})=\text{\rm
SL}_{2}({\mathbb{R}})/\\{\pm I\\}$ and the proper Lorentz group $\text{\rm
SO}^{+}(1,2)$ of linear transformations leaving invariant the quadratic form
$x^{2}-y^{2}-z^{2}$, except for a global change of sign. In general, only the
matrices in $\text{\rm SO}^{+}(1,2)\cap\mathcal{M}_{3}({\mathbb{Z}})$ with
$\mathcal{M}_{3}({\mathbb{Z}})$ the set of integral $3\times 3$ matrices,
apply integral solutions into integral solutions. In the jargon of the
quadratic form theory, these matrices are called _integral automorphs_.

Clearly, the fundamental problem is to characterize the preimage of $\text{\rm
SO}^{+}(1,2)\cap\mathcal{M}_{3}({\mathbb{Z}})$ by the isomorphism (4.2) in a
simple way. Note that it is unclear whether any integral matrix in the right
side is obtained from an integral matrix, an element of $\text{\rm
PSL}_{2}({\mathbb{Z}})$, in the left side444Note that the image of some
matrices in $\text{\rm PSL}_{2}({\mathbb{Z}})$ has semi-integral entries. For
instance, this is the case for $a=2$, $b=c=d=1$.. The nature of the allowed
values of these four parameters has to be established as follows. First, let
$A$ be the matrix in the image of (4.2). If
$A\in\mathcal{M}_{3}({\mathbb{Z}})$, as its entries are integers, it is seen
that

$2a^{2}=a_{11}+a_{22}-a_{12}-a_{21}\in{\mathbb{Z}}.$

Other choices of the signs prove that in general

$2a^{2},2b^{2},2c^{2},2d^{2}\in{\mathbb{Z}}.$

On the other hand, using the rest of the entries and the determinant equation
$ad-bc=1$, it is seen that

$2ab=a_{23}-a_{13}\in{\mathbb{Z}}$

and, similarly,

$2ac,\,2ad,\,2bc,\,2bd,\,2cd\in{\mathbb{Z}}.$

Considering $a_{ij}\pm a_{i^{\prime}j^{\prime}}$ with
$i,j,i^{\prime},j^{\prime}\in\\{1,2\\}$, it is deduced that
$2a^{2},\,2b^{2},\,2c^{2},\,2d^{2}$ are all even or all odd. Consider first
the even case. If one selects a nonzero variable, say $a$, then $a=s\sqrt{R}$
with $R$ squarefree. This, together with $2ab,2ac,2ad\in{\mathbb{Z}}$ implies
that $b$, $c$ and $d$ are integers when divided by $\sqrt{R}$. Finally, the
determinant equation $ad-bc=1$ implies $R=1$ and it is concluded
$a,b,c,d\in{\mathbb{Z}}$. On the other hand, the extra condition $2\mid
a+b+c+d$ is needed to assure $a_{11},a_{12},a_{21},a_{22}\in{\mathbb{Z}}$. In
these terms, the elements of the sometimes called $\theta$-group

$\Gamma_{\theta}=\left\\{\begin{pmatrix}a&b\\\ c&d\end{pmatrix}\in\text{\rm
PSL}_{2}({\mathbb{Z}})\,:\,2\mid a+b+c+d\right\\}$

are obtained. For the odd case instead, one may consider the quantities
$2(a\sqrt{2})^{2}$, $2(b\sqrt{2})^{2}$, $2(c\sqrt{2})^{2}$, $2(d\sqrt{2})^{2}$
which are clearly even. By repeating the argument it is deduced that
$a\sqrt{2}$, $b\sqrt{2}$, $c\sqrt{2}$, $d\sqrt{2}$ are odd integers. It is not
difficult to see that the matrices in $\text{\rm PSL}_{2}({\mathbb{R}})$
having odd entries when multiplied by $\sqrt{2}$ belong to the coset

$\frac{1}{\sqrt{2}}\begin{pmatrix}1&-1\\\ 1&1\end{pmatrix}\Gamma_{\theta}$

and, conversely, any matrix in it has that property. Therefore the isomorphism

$F:\,\Gamma\longrightarrow\text{\rm
SO}^{+}(1,2)\cap\mathcal{M}_{3}({\mathbb{Z}})$

holds where

$\Gamma=\Gamma_{\theta}\cup C_{\theta}\qquad\text{with}\quad
C_{\theta}=\frac{1}{\sqrt{2}}\begin{pmatrix}1&-1\\\
1&1\end{pmatrix}\Gamma_{\theta}.$

By (4.1), and by making the corresponding variable identification, it follows
that $x=s$, $y=k_{0}$, $z=\ell_{0}$ is a valid solution of
$Dx^{2}-y^{2}-z^{2}=q^{2}$ and

$F(\gamma)\begin{pmatrix}s\\\ k_{0}\\\
\ell_{0}\end{pmatrix}\qquad\text{with}\quad\gamma\in\Gamma,$ (4.3)

gives a family of infinitely many of them for $D=1$. Generating elements of
$\Gamma$ is a simple task because the extended Euclidean algorithm produces
easily integer solutions of $aX-bY=1$ (even under a parity condition) to
obtain matrices in $\Gamma_{\theta}$ composed by $a$, $b$, $c=Y$, $d=X$.
Another possibility is to use group generators of $\Gamma_{\theta}$ to get
generators of $\Gamma$.

To give an example of the above procedure, consider the choice $q=\sqrt{6}$,
$k_{0}=1$, $\ell_{0}=3$ gives $D=1$ and $s=4$ in (4.1). By applying
$F(\gamma)$ with $\gamma$ each of the 12 matrices in $\Gamma$ with
$|a|,|b|,|c|,|d|\leq 2$, they are obtained the following valid values
$(k_{j},\ell_{j})$ satisfying (2.3)

$\pm(1,3),\quad\pm(1,-3),\quad\pm(3,-1),\quad\pm(3,7),\quad\pm(13,9),\quad\pm(15,-13).$

The case $D>1$ is much more difficult although it runs in similar lines.
Fortunately, there exists literature about the subject that dates from more
than a century ago [14]. The initial idea is that the change of variables
$x\mapsto x/\sqrt{D}$ passes the quadratic form $Dx^{2}-y^{2}-z^{2}$ into
$x^{2}-y^{2}-z^{2}$. It induces a change in the image of (4.2) and now the
isomorphism between $\text{\rm PSL}_{2}({\mathbb{R}})=\text{\rm
SL}_{2}({\mathbb{R}})/\\{\pm I\\}$ and the proper Lorentz group preserving
$Dx^{2}-y^{2}-z^{2}$, except for a global change of sign, becomes

$\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\longmapsto\begin{pmatrix}\frac{1}{2}\big{(}a^{2}+b^{2}+c^{2}+d^{2}\big{)}&\frac{1}{2\sqrt{D}}\big{(}-a^{2}+b^{2}-c^{2}+d^{2}\big{)}&-\frac{1}{\sqrt{D}}(ab+cd)\\\
\frac{\sqrt{D}}{2}\big{(}-a^{2}-b^{2}+c^{2}+d^{2}\big{)}&\frac{1}{2}\big{(}a^{2}-b^{2}-c^{2}+d^{2}\big{)}&ab-
cd\\\ -(ac+bd)\sqrt{D}&ac-bd&bc+ad\end{pmatrix}.$

The crucial question is again how to characterize $a,\ b,\,c,\,d$ in the
initial matrix having an integral matrix as image. A simplification is
reached, re-parametrizing the initial matrix as

$M=\begin{pmatrix}-a-c\sqrt{D}&-b-d\sqrt{D}\\\
b-d\sqrt{D}&-a+c\sqrt{D}\end{pmatrix}$ (4.4)

because it rules out the square roots in the image, namely, the previous
isomorphism is

$F\,:\,M\longmapsto\begin{pmatrix}a^{2}+b^{2}+D(c^{2}+d^{2})&-2ac+2bd&-2bc-2ad\\\
-2D(ac+bd)&a^{2}-b^{2}+D(c^{2}-d^{2})&2ab+2Dcd\\\ 2D(bc-
ad)&-2ab+2Dcd&a^{2}-b^{2}-D(c^{2}-d^{2})\end{pmatrix}.$ (4.5)

In these terms it is plain to see that for $a,b,c,d\in{\mathbb{Z}}$ it is
obtained an integral matrix. In fact, as all the entries are quadratic, one
can expect that it is possible to introduce $\sqrt{2}$ denominators in $M$ if
the parity conditions cancel the $1/2$ in the diagonal entries of the image,
as above. Elaborating this idea, it can be proved [14] (cf. [19]) that the
case $4\mid D-1$ parallels the case $D=1$. Namely, the group is composed by
two parts that reduce to $\Gamma_{\theta}$ and the companion coset
$C_{\theta}$ for $D=1$. Specifically,

$\Gamma=\Gamma^{\prime}_{\theta}\cup C^{\prime}_{\theta}$

where

$\Gamma^{\prime}_{\theta}=\big{\\{}M\in\text{\rm
PSL}_{2}({\mathbb{R}})\,:\,a,b,c,d\in{\mathbb{Z}}\big{\\}}\qquad\text{and}\qquad
C^{\prime}_{\theta}=\Big{\\{}\frac{1}{\sqrt{2}}M\in\text{\rm
PSL}_{2}({\mathbb{R}})\,:\,a,b,c,d\in{\mathbb{Z}}\Big{\\}}.$

Here $M$ is like in (4.4).

If $4\nmid D-1$ then it is necessary to add a new set, which is a kind of
semi-integral version of $\Gamma^{\prime}_{\theta}$ namely

$\Gamma=\Gamma^{\prime}_{\theta}\cup C^{\prime}_{\theta}\cup
C^{\prime\prime}_{\theta}$

where $\Gamma^{\prime}_{\theta}$ and $C^{\prime}_{\theta}$ are as before and

$C^{\prime\prime}_{\theta}=\Big{\\{}\frac{1}{2}M\in\text{\rm
PSL}_{2}({\mathbb{R}})\,:\,a,b,c,d\text{ odd integers}\Big{\\}}.$

In fact, the distinction between $4\mid D-1$ and $4\nmid D-1$ is a little
artificial (only to emphasize the analogy of the former case with $D=1$)
because $C^{\prime\prime}_{\theta}$ is empty for $4\mid D-1$. It is also empty
for $D$ even because, by computing the determinant, $\frac{1}{2}M\in\text{\rm
PSL}_{2}({\mathbb{R}})$ implies
$\frac{1}{2}(a^{2}+b^{2})=2+\frac{1}{2}D(c^{2}+d^{2})$. If $a,\,b,\,c,\,d$ are
odd then the left side is odd and the right side is even.

In any case, (4.3) with $F$ given by (4.5) produces infinitely many solutions
of $Dx^{2}-y^{2}-z^{2}=q^{2}$ showing that $\mathcal{N}_{0}$ is infinite. As
in the one dimensional case, the theory assures [6, Lemma 6.1] that there is a
finite number of infinite families of solutions. That is to say, with a finite
number of choices of $(s,k_{0},\ell_{0})$ all the needed solutions to
construct $\mathcal{N}_{0}$ are reached with (4.3).

Let us close this section with some comments about
$q^{2}\in{\mathbb{Q}}-{\mathbb{Z}}$. Instead of deepening into the theory, we
only show a shortcut to deduce $\\#\mathcal{N}_{0}=\infty$ from our previous
knowledge about $q^{2}\in{\mathbb{Z}}$.

If $q^{2}\in{\mathbb{Q}}-{\mathbb{Z}}$, (4.1) must be replaced by

$k_{0}^{2}+\ell_{0}^{2}+q^{2}=\frac{Ds^{2}}{D_{*}s_{*}^{2}}\qquad\text{with
$D$ and $D_{*}$ squarefree and coprime}.$

Clearing the denominator and recalling (2.3), it is deduced

$D(a_{j}s/b_{j})^{2}-D_{*}s_{*}^{2}k_{j}^{2}-D_{*}s_{*}^{2}\ell_{j}^{2}=q^{2}s_{*}^{2}.$

Hence

$\mathcal{N}_{0}=\big{\\{}(y,z)\,:\,Dx^{2}-D_{*}s_{*}^{2}y^{2}-D_{*}s_{*}^{2}z^{2}=q^{2}s_{*}^{2}\text{
with }(x,y,z)\in{\mathbb{Z}}^{3}\big{\\}}.$

In principle, one could study the group leaving invariant the quadratic form
$Dx^{2}-D_{*}s_{*}^{2}y^{2}-D_{*}s_{*}^{2}z^{2}$ and proceed as before [14],
but some technical issues appear because $D_{*}s_{*}^{2}$ is not squarefree if
$s_{*}>1$. A trick to overcome this problem is to re-write the equation as

$DD_{*}x^{2}-(D_{*}s_{*}y)^{2}-(D_{*}s_{*}z)^{2}=q^{2}s_{*}^{2}.$

Then one has to look for integral solutions of
$DD_{*}x^{2}-Y^{2}-Z^{2}=q^{2}s_{*}^{2}$ with $Y$ and $Z$ multiples of
$D_{*}s_{*}$. Note that $DD_{*}$ is squarefree then (4.3) with $F$ as in (4.5)
changing $D$ by $DD_{*}$ and $(s,k_{0},\ell_{0})$ by
$(s,D_{*}s_{*}k_{0},D_{*}s_{*}\ell_{0})$ produces infinitely many solutions
$(x,Y,Z)$. The issue is to prove that infinitely many of them verify that
their two last coordinates are multiples of $D_{*}s_{*}$. As the $3\times 3$
matrices of determinant $1$ modulo $D_{*}s_{*}$ form a finite group, there
exists $N\in{\mathbb{Z}}^{+}$, its order, such that

$\big{(}F(\gamma)\big{)}^{N}=F(\gamma^{N})\equiv\text{Id}\mod{D_{*}s_{*}}.$

Hence

$F(\gamma^{N})\begin{pmatrix}s\\\ D_{*}s_{*}k_{0}\\\
D_{*}s_{*}\ell_{0}\end{pmatrix}\equiv\begin{pmatrix}s\\\ 0\\\
0\end{pmatrix}\mod{D_{*}s_{*}}.$

As $\\{\gamma^{N}\,:\,\gamma\in\Gamma\\}$ is infinite, there are infinite
solutions with $D_{*}s_{*}\mid Y,Z$ and $\mathcal{N}_{0}$ is infinite too.

Although it is possible to obtain an explicit formula for $N$, for moderate
values of $D_{*}s_{*}$ may be more useful in practice to compute the image of
successive powers of $\gamma$ until its image by $F$ fulfills the divisibility
conditions.

##  5 Some numerical simulations

The purpose of this section is to explain some numerical examples in detail
and to plot some quantum carpets555The Matlab/Octave code generating the
figures and the full size color images will be posted in the web page
https://matematicas.uam.es/~fernando.chamizo/dark/d_talbot.html along with
more related material.. The latter in the 1D case corresponds to density plots
of $|\Phi(\phi_{2},t)|$, the square root of the probability density for the
state in (2.2) with $t$ in the horizontal axis and $\phi_{2}$ in the vertical
axis. The only reason to prefer $|\Psi|$ instead of the more natural
$|\Psi|^{2}$ is to avoid the saturation giving less appealing images. These
plots involve a choice of the coefficients $\lambda_{\vec{n}}$ indicated in
each case. The angle is considered in the range $[-\pi,\pi]$, as mentioned in
the first section. For the time the natural scale is $T_{\text{rev}}$, the
revival time. Hence $[0,1]$ in the plot must be interpreted as
$[0,T_{\text{rev}}]$ in usual units. In the 2D case a density plot of
$|\Phi(\phi_{1},\phi_{2},t)|$ is not possible because it would require three
dimensions. It will replaced below by plots for close fixed values of $t$ to
exemplify the evolution of the state.

Let us start with a 1D example having a finite $\mathcal{N}_{0}$. Consider
$\ell_{0}=3$ and $q=\frac{3}{2}\sqrt{21}$. Then
$\ell_{0}^{2}+q^{2}=15^{2}/2^{2}$ which corresponds in (3.7) to $D=D_{*}=1$,
$s_{*}=3$. Using the explicit formula for $\mathcal{N}_{0}$ in terms of the
divisors of $q^{2}s_{*}^{2}=189$, it is obtained

$\mathcal{N}_{0}=\big{\\{}(0,\pm 3),(0,\pm 5),(0,\pm 15),(0,\pm 47)\big{\\}}.$

Here all the divisors of $189$ satisfy the extra condition of $189/d-d$ being
multiple of $4$. Computing the $a_{j}/b_{j}$, given by
$\sqrt{\ell_{j}^{2}+q^{2}}/\sqrt{\ell_{0}^{2}+q^{2}}$ according to (3.1), they
are $1/1$, $17/15$, $11/5$ and $19/3$. Hence $L=15$ and the revival time is
$T_{\text{rev}}=30\pi/\omega_{\vec{n}_{0}}$. It remains the same when we omit
$(0,\pm 47)$. The first two figures are the plots of the corresponding quantum
carpets in these situations with $1$ coefficients.

---
Fig. 1. Full $\mathcal{N}_{0}$ with $\lambda_{\vec{n}}=1$.
---
Fig. 2. $\mathcal{N}_{0}-\\{(0,\pm 47)\\}$ with $\lambda_{\vec{n}}=1$.

Let us now recover the example $\ell_{0}=5$ and $q=\sqrt{2}$ to show how to
obtain the result claimed in (3.6). Recall $s=D=3$. It was already shown that
there is only a family of solutions, so
$\mathcal{N}_{0}=\\{(0,x_{n})\,:\,n\in{\mathbb{Z}}\\}$ where, according to
(3.5),

$x_{n}+y_{n}\sqrt{3}=\pm(5+3\sqrt{3})(2+\sqrt{3})^{n}.$

Taking apart the $\pm$, it produces
$x_{n+1}+y_{n+1}\sqrt{3}=(2+\sqrt{3})(x_{n}+y_{n}\sqrt{3})$ that leads to the
first order recurrence

$(x_{n+1},y_{n+1})=(2x_{n}+3y_{n},x_{n}+2y_{n}).$

Considering another step and eliminating $y_{n}$ and $y_{n+1}$, it is obtained
the second order recurrence

$x_{n+2}=4x_{n+1}-x_{n}\qquad\text{with}\quad x_{0}=5,\
x_{1}=19\quad\text{for}\quad n\in{\mathbb{Z}}.$

Noting $x_{-3-n}=-x_{-n}$ for $n=0,1$, it follows for any $n$. Hence defining
$c_{n}=x_{n-1}$ for $n\geq 0$, (3.6) follows. In the same way,
$c_{n}^{\prime}=y_{n-1}$ satisfies
$c^{\prime}_{n+2}=4c^{\prime}_{n+1}-c^{\prime}_{n}$ too with starting values
$c_{0}^{\prime}=1$, $c_{1}^{\prime}=2$. An easy inductive argument using the
recurrence shows that $c_{n}^{\prime}$ and $c_{n+1}^{\prime}$ are coprime (in
particular, at least one is not divisible by $3$) then for any $\mathcal{N}$
containing $(0,c_{n})$ and $(0,c_{n+1})$ for some $n$ we have $L=3$ and the
revival time is

$T_{\rm rev}=\frac{2\pi R}{c\sqrt{3}}=\frac{6\pi}{\omega_{\vec{n}_{0}}}.$

Due to the exponential growth of $c_{n}$, even for fairly small cardinalities
of $\mathcal{N}$, the state (2.2) shows large variations for tiny changes of
$\phi_{2}$ and $t$ and the density plot is close to be a cloud of random
points. In this context it is natural to let $\lambda_{\vec{n}}$ to decay with
$|\vec{n}|$. In Fig. 3 and Fig. 4 there are examples choosing $\ell_{j}$ the
smallest values of $c_{n}$ and
$\lambda_{\vec{n}}=|\vec{n}|^{-1/4}=|\ell_{j}|^{-1/4}$.

---
Fig. 3. $\ell_{j}=\pm 1,\pm 5,\pm 19$, $\lambda_{\vec{n}}=|\ell_{j}|^{-1/4}$.
---
Fig. 4. $\ell_{j}=\pm 1,\pm 5,\pm 19,\pm 71$,
$\lambda_{\vec{n}}=|\ell_{j}|^{-1/4}$.

As a matter of fact, for each $t$ fixed, when $\\#\mathcal{N}\to\infty$ it is
obtained a lacunary series and according to known facts in fractal geometry
[13, §11.3] the graphs of the real and imaginary part give rise to a fractal,
a curve with fractional box dimension.

In this and other examples, it is apparent an underlying structure based on
interwoven straight lines with 2 different slopes. The explanation is that in
(2.2), due to the exponential growth, $\sqrt{\vec{n}^{2}+q^{2}}/n_{2}$ is very
early close to the sign of $n_{2}$ and the pairs of lines correspond to the
light cones. In fact if all the $\ell_{j}$ are chosen to be positive, commonly
the density plot looks like a dull collection of oblique parallel bands.

To mask the light cones, a possibility is to take $q^{2}$ large to avoid
$\sqrt{\ell_{j}^{2}+q^{2}}/|\ell_{j}|$ to be close to 1, at least for the
first values of $\ell_{j}$. With this idea in mind and to illustrate the
appearance of several families of solutions, consider the example

$\vec{n}_{0}=(0,3)\qquad\text{with}\qquad q=\sqrt{791}.$

From (3.2), $D=2$ and $s=20$. Consequently, the generalized Pell’s equation is
$x^{2}-2y^{2}=-791$ and by (3.5), it admits a family of solutions described by

$x+y\sqrt{2}=\pm\big{(}3+20\sqrt{2}\big{)}\big{(}3+2\sqrt{2}\big{)}^{n}\qquad\text{with}\quad
n\in{\mathbb{Z}}$

because $(3,2)$ is a minimal positive solution of $x^{2}-2y^{2}=1$. But
solutions $(x,y)$ in other families may exist. The condition
$0<a<791(1+2\sqrt{2\cdot 3})/2$ for a positive solution $(a,b)$ to be the seed
of a new family, by direct calculations, only gives $(19,24)$. So it is
obtained the second family

$x+y\sqrt{2}=\pm\big{(}19+24\sqrt{2}\big{)}\big{(}3+2\sqrt{2}\big{)}^{n}\qquad\text{with}\quad
n\in{\mathbb{Z}}$

and any solution of $x^{2}-2y^{2}=-791$ belongs to one of them.

Taking $n\geq 0$, the first positive values for $x$ in the first family are

$3,89,531,3097,18051,105209,613203,\dots$

which satisfy the recurrence $c_{k+1}=6c_{k}-c_{k-1}$. For $n<0$ the values
are

$71,429,2503,14589,85031,495597,2888551,\dots$

and, naturally, they satisfy the same recurrence with different starting
values. In the same way, for the second family it is obtained for $n\geq 0$

$19,153,899,5241,30547,178041,1037699,\dots$

and for $n<0$

$39,253,1479,8621,50247,292861,1706919,\dots$

still satisfying the same recurrence. Then

$\mathcal{N}_{0}=\big{\\{}(0,\pm c_{k}^{(j)})\,:\,1\leq j\leq 4,\
k\in{\mathbb{Z}}_{\geq 0}\big{\\}}$

where $c_{k+1}^{(j)}=6c_{k}^{(j)}-c_{k-1}^{(j)}$ with the starting values

$\begin{array}[]{r|c|c|c|c|}&j=1&j=2&j=3&j=4\\\
\hline\cr(c_{k}^{(0)},c_{k}^{(1)})&(3,89)&(71,429)&(19,153)&(39,253)\\\
\hline\cr\end{array}$

In Fig. 5 it is considered the subset

$\mathcal{N}=\big{\\{}(0,\pm 3),(0,\pm 19),(0,\pm 39)\big{\\}}$

with the coefficients $\lambda_{\vec{n}}=|\vec{n}|^{-1/4}$, as before, and in
Fig. 6 it is added $(0,\pm 71)$. In both cases the revival time is
$T_{\text{\rm rev}}=20\pi/\omega_{\vec{n}_{0}}$ because the $a_{j}/b_{j}$ are
$1/1$, $6/5$, $17/10$ and $(0,\pm 71)$ adds $27/10$ that does not change the
value of $L$ which is $10$.

---
Fig. 5. $\ell_{j}=\pm 3,\pm 19,\pm 39$, $q^{2}=791$.
---
Fig. 6. $\ell_{j}=\pm 3,\pm 19,\pm 39,\pm 71$, $q^{2}=791$.

Note that, as expected, the clear line structure in the previous example
coming from the light cones has been broken.

Recall that in the 2D setting a way to escape from the condition
$q^{2}\in{\mathbb{Q}}$ keeping $\mathcal{N}_{0}$ nontrivial was to consider
numbers with many representations as a sum of two squares. According to the
theory [16, §16.9], if $N$ is a product of $k$ distinct primes of the form
$4n+1$ then the number of representations of $N$ as a sum of two squares is
$2^{k+2}$. Choosing $k_{0}^{2}+\ell_{0}^{2}=N$, in this way
$\\#\mathcal{N}_{0}=2^{k+2}$ for any $q^{2}\not\in{\mathbb{Q}}$ while its
energy set reduces to the single element $E=c\hbar\sqrt{N+q^{2}}/R$. Many
representations are obtained applying the following 8 obvious symmetries:

$(n_{1},n_{2})\mapsto(\pm n_{1},\pm
n_{2})\qquad\text{and}\qquad(n_{1},n_{2})\mapsto(\pm n_{2},\pm n_{1}).$

For instance, for $N=1105=5\cdot 13\cdot 17$, the $2^{3+2}=32$ representations
come from applying the symmetries to the 4 pairs

$(4,33),\quad(9,32),\quad(12,31),\quad(23,24).$

The value of $q^{2}$ does not affect to $\mathcal{N}_{0}$ and letting it grow,
only the first coordinate in (2.1) is relevant, giving in the limit

$\sum\Psi^{1}_{(\pm n_{1},\pm n_{2})}+\sum\Psi^{1}_{(\pm n_{2},\pm
n_{1})}=4e^{-iEt/\hbar}\big{(}\cos(n_{1}\phi_{1})\cos(n_{2}\phi_{2})+\cos(n_{2}\phi_{1})\cos(n_{1}\phi_{2})\big{)}$

where the sum is over all sign combinations. Although this may seem very
simple, the chaotic distribution of $n_{1}$ and $n_{2}$ for $N$ large leads to
some hard questions in mathematical physics. Namely, even in this toral
situation, some general conjectures [2] about the geometry of the _nodal
lines_ (the zero level sets) of the states corresponding to the a given energy
remain unsolved.

To illustrate the richness of the situation, in Fig. 7 there is a plot of the
nodal lines of the first coordinate of the state (2.2) for $\mathcal{N}_{0}$
as in the previous example with $k_{0}^{2}+\ell_{0}^{2}=1105$ and
$\lambda_{\vec{n}}=1$ in the range $\phi_{1},\phi_{2}\in[-\pi/2,\pi/2]$.

---
Fig. 7. Nodal lines $k_{0}^{2}+\ell_{0}^{2}=1105$, $\lambda_{\vec{n}}=1$.
---
Fig. 8. $|\Phi(\phi_{1},\phi_{2})|^{0.5}$, $\lambda_{\vec{n}}=1$.

The image in Fig. 8 is the corresponding density plot of $|\Phi|^{0.5}$. As
$e^{-iEt/\hbar}$ is constant, $|\Phi|$ does not depend on $t$. The exponent
$0.5$ is only for visualization purposes. It has been introduced to mollify
the effect of the peak at $\phi_{1}=\phi_{2}=0$.

Consider now the example $(k_{0},\ell_{0})=(1,2)$ with $q^{2}=3$. The
decomposition (4.1) gives $D=s=2$ and $\mathcal{N}_{0}$ is composed by the
$(y,z)$ pairs got from the solutions of $2x^{2}-y^{2}-z^{2}=3$. If one only
want to compute a small subset of $\mathcal{N}_{0}$ then one can proceed by
direct inspection instead of using the parametrization by the Fuchsian groups.
Let us take

$\mathcal{N}=\big{\\{}(\pm 1,\pm 2),\,(\pm 2,\pm 1),\,(\pm 2,\pm 5),\,(\pm
5,\pm 2),\,(\pm 2,\pm 11),\,(\pm 11,\pm 2),\,(\pm 5,\pm 10),\,(\pm 10,\pm
5)\big{\\}}$

which comes from $\big{\\{}(1,2),\,(2,5),\,(2,11),\,(5,10)\big{\\}}$ applying
the symmetries. The resulting sequence
$a_{j}/b_{j}=\sqrt{(k_{j}^{2}+\ell_{j}^{2}+3)/8}$ is 1, 2, 4, 4. All the
denominators are 1, then the revival time is $T_{\text{\rm
rev}}=2\pi/\omega_{\vec{n}_{0}}$.

As mentioned before, the density plot of the associated state $\Phi$ is not
feasible because there are three variables $(\phi_{1},\phi_{2},t)$. To
illustrate the situation, Fig. 9 is composed by $4$ images showing the density
plots in $(\phi_{1},\phi_{2})$ for $t=0,\,T_{\text{\rm rev}}/10,\,T_{\text{\rm
rev}}/5$ and $T_{\text{\rm rev}}/10$.

---
Fig. 9. $|\Phi(\phi_{1},\phi_{2},kT_{\text{\rm rev}}/10)|$ for $k=0,1,2,3$,
$\phi_{1},\phi_{2}\in[-\pi,\pi]$, with $\lambda_{\vec{n}}=1$.

Playing with different values, it seems that the sensitivity in $t$ is not
uniform. For instance, the density plots for $t=0$ and $t=T_{\text{\rm
rev}}/20$ are quite similar while there are noticeable differences between
$t=T_{\text{\rm rev}}/2$ and $t=T_{\text{\rm rev}}/2-T_{\text{\rm rev}}/20$.
Probably part of the explanation is the large peak corresponding to
$\phi_{1}=\phi_{2}=t=0$ which masks relative differences.

##  6 Discussion of the results

In the present paper, the possible revivals for a relativistic free fermion
ruled by the Dirac equation in a one dimensional and a two dimensional torus,
were characterized. This was achieved by use of arithmetic tools such as
Pell’s equations or Fuchsian groups. The revivals shown here are exact, that
is, there is no need to take relativistic limit in order to detect them. These
results generalize those found in [25], who found revivals related to
Pythagorean triples. In the present paper, all the possible state exhibiting
revivals were characterized. The interesting feature of these revivals is that
they are known to hold for massless wave equations or for the standard
Schrödinger equation. However in the present context the dispersion relation
is modified to $E=\sqrt{m^{2}+p^{2}}$ and still the revivals are exact.
Although the torus topology may seem non realistic at first sight, it should
be remarked that the Dirac equation in several topologies may have
applications in solid state physics. The study of revivals in this context is
encouraging and perhaps may lead to interesting experimental results such as
the ones reported in [26].

## Acknowledgments

O.P. S. is grateful to the Universidad Autónoma and to the ICMAT de Madrid,
where this work was performed, by their hospitality. The present work is
supported by CONICET, Argentina and by the Grant PICT 2020-02181. This project
has received funding from the European Union’s Horizon 2020 research and
innovation program under the Marie Skłodowska-Curie grant agreement No 777822.
F. Ch. is partially supported by the PID2020-113350GB-I00 grant of the MICINN
(Spain) and by “Severo Ochoa Programme for Centres of Excellence in R&D”
(CEX2019-000904-S).

## References

  * [1] M. V. Berry.  Quantum fractals in boxes.  J. Phys. A, 29(20):6617–6629, 1996.
  * [2] M. V. Berry.  Statistics of nodal lines and points in chaotic quantum billiards: perimeter corrections, fluctuations, curvature.  J. Phys. A, 35(13):3025–3038, 2002.
  * [3] M. V. Berry and S. Klein.  Integer, fractional and fractal Talbot effects.  J. Modern Opt., 43(10):2139–2164, 1996.
  * [4] M. V. Berry, I. Marzoli, and W. P. Schleich.  Quantum carpets, carpets of light.  Physics World, 14(6):39, jun 2001.
  * [5] J. Bourgain and Z. Rudnick.  On the geometry of the nodal lines of eigenfunctions of the two-dimensional torus.  Ann. Henri Poincaré, 12(6):1027–1053, 2011.
  * [6] J. W. S. Cassels.  Rational quadratic forms, volume 13 of London Mathematical Society Monographs.  Academic Press, Inc. [Harcourt Brace Jovanovich, Publishers], London-New York, 1978.
  * [7] F. Chamizo and O. P. Santillán.  About the quantum Talbot effect on the sphere.  J. Phys. A, 56(25):Paper No. 255302, 20, 2023.
  * [8] H. Cohn.  Advanced number theory.  Dover Books on Advanced Mathematics. Dover Publications, Inc., New York, 1980.  Reprint of ıt A second course in number theory, 1962.
  * [9] K. Conrad.  Pell’s equation, II.  https://kconrad.math.uconn.edu/blurbs/, 2022.
  * [10] G. N. Copley.  Recurrence relations for solutions of Pell’s equation.  Amer. Math. Monthly, 66:288–290, 1959.
  * [11] W. Duke.  Hyperbolic distribution problems and half-integral weight Maass forms.  Invent. Math., 92(1):73–90, 1988.
  * [12] D. Eceizabarrena.  The Talbot effect as the fundamental solution to the free Schrödinger equation.  Port. Math., 78(2):233–253, 2021.
  * [13] K. Falconer.  Fractal geometry.  John Wiley & Sons, Ltd., Chichester, 1990.  Mathematical foundations and applications.
  * [14] R. Fricke and F. Klein.  Vorlesungen über die Theorie der automorphen Functionen. Erster Band; Die gruppentheoretischen Grundlagen. Leipzig. B. G. Teubner, 1897.
  * [15] J. H. Hannay and A. Lockwood.  The quantum Talbot effect on a sphere.  J. Phys. A, 41(39):395205, 9, 2008.
  * [16] G. H. Hardy and E. M. Wright.  An introduction to the theory of numbers.  The Clarendon Press, Oxford University Press, New York, fifth edition, 1979.
  * [17] M. J. Jacobson, Jr. and H. C. Williams.  Solving the Pell equation.  CMS Books in Mathematics/Ouvrages de Mathématiques de la SMC. Springer, New York, 2009.
  * [18] R.D. Klauber.  Student Friendly Quantum Field Theory: Basic Principles and Quantum Electrodynamics.  Sandtrove Press, 2013.
  * [19] W. Magnus.  Noneuclidean tesselations and their groups, volume 61 of Pure Appl. Math., Academic Press.  Academic Press, New York, NY, 1974.
  * [20] I. Marzoli, A. E. Kaplan, F. Saif, and W. P. Schleich.  Quantum carpets of a slightly relativistic particle.  Fortschr. Phys., 56(10):967–992, 2008.
  * [21] I. Marzoli, F. Saif, I. Bialynicki-Birula, O. M. Friesch, A. E. Kaplan, and W. P. Schleich.  Quantum carpets made simple.  Acta Phys. Slov., 48:323–333, 1998.
  * [22] P. J. Olver.  Dispersive quantization.  Amer. Math. Monthly, 117(7):599–610, 2010.
  * [23] K. I. Oskolkov.  The Schrödinger density and the Talbot effect.  In Approximation and probability, volume 72 of Banach Center Publ., pages 189–219. Polish Acad. Sci. Inst. Math., Warsaw, 2006.
  * [24] J. E. Shockley.  Introduction to number theory.  Holt, Rinehart and Winston, Inc., New York-Toronto-London, 1967.
  * [25] P. Strange.  Relativistic Quantum Revivals.  Phys. Rev. Lett., 104(12):120403, mar 2010.
  * [26] J. J. Torres and E. Romera.  Wave packet revivals in a graphene quantum dot in a perpendicular magnetic field.  Phys. Rev. B, 82:155419, Oct 2010.


# On the Ergodic Capacity of Reconfigurable Intelligent Surface (RIS)-Aided
MIMO Channels

Chongjun Ouyang, Sheng Wu, Chunxiao Jiang,
Yuanwei Liu, and Hongwen Yang C. Ouyang, S. Wu, and H. Yang are with the
School of Information and Communication Engineering, Beijing University of
Posts and Telecommunications, Beijing, 100876, China (e-mail: {DragonAim,
thuraya, yanghong}@bupt.edu.cn).C. Jiang is with the Tsinghua Space Center,
Tsinghua University, Beijing, 100084, China, and also with the Beijing
National Research Center for Information Science and Technology, Beijing,
100084, China (e-mail: jchx@tsinghua.edu.cn).Y. Liu is with the School of
Electronic Engineering and Computer Science, Queen Mary University of London,
London, UK (e-mail: yuanwei.liu@qmul.ac.uk).

###### Abstract

Reconfigurable intelligent surfaces (RISs) have emerged as a promising
technique to enhance the system spectral efficiency. This letter investigates
the ergodic channel capacity (ECC) of an RIS-aided multiple-input multiple-
output channel under the assumption that the transmitter-RIS, RIS-receiver,
and transmitter-receiver channels contain deterministic line-of-sight paths.
Novel expressions are derived to characterize the upper and lower bounds of
the ECC. To unveil more system insights, asymptotic analyses are performed to
the system ECC in the limit of large signal-to-noise ratio (SNR) and number of
reflecting elements (REs). Theoretical analyses suggest that the RIS’s
deployment can shape the ECC curve by influencing its high-SNR power offset
and the ECC can get improved by increasing the number of REs.

###### Index Terms:

Ergodic capacity, line-of-sight paths, multiple-input multiple-output,
reconfigurable intelligent surface.

##  I Introduction

Reconfigurable intelligent surfaces (RISs) can enhance the system spectral
efficiency (SE) with low energy consumption and deployment costs, which is
believed to be a promising technology in next-generation networks [1].
Specifically, by properly programming its adopted phase shifts, an RIS can
control the wireless physical layer channel intelligently and then give it a
performance boost [2]. One of the key scenarios where RIS shines, is for SE
improvements in multi-antenna systems [3].

Recently, there has been an increasing interest in performing theoretical
analyses on the ergodic SE (ESE) of RIS-aided multi-antenna systems [4, 5, 6,
8, 7]. For example, a closed-form lower bound was derived to characterize the
ESE of an RIS-aided multiuser multiple-input single-output (MISO) channel
under Rayleigh fading [4]. Yet, the RIS is usually deployed on the facade of a
high building in the vicinity of transceivers, which yields deterministic
line-of-sight (LoS) paths in the RIS-related channels. Thus, it makes more
sense to characterize the ESE of RIS-aided systems under fading models that
consider the influence of deterministic LoS paths, such as Rician fading.
Particularly, a tight upper bound of the ESE of a single-user RIS-aided MISO
Rician channel was formulated in [5], which was further extended to the
multiuser case [6]. Moreover, the ergodic channel capacity (ECC) of an RIS-
aided multiple-input multiple-output (MIMO) Rician channel was investigated in
[7] under the assumption of single-stream transmissions. In fact, the MIMO
channel therein can be considered equivalent to a MISO one. As a supplement of
this work, a large-system approximation of ECC of an RIS-aided MIMO Rician
channel under multi-stream transmission was derived by the replica method [8].
Unfortunately, the ECC’s calculation therein requires the solutions of a
series of fixed-point equations and it is challenging to unveil more insights
on the effect of RIS from the analytical results in [8].

To portray the effect of RIS on the ECC and thus fill the aforementioned
research gap, this letter evaluates the SE performance of an RIS-aided MIMO
Rician channel under multi-stream transmission by formulating upper and lower
bounds of the ECC. For gleaning further insights, we perform asymptotic
analyses to the ECC in the limit of large signal-to-noise ratio (SNR) and
number of reflecting elements (REs). We observe that the RIS can influence the
asymptotic ECC through its high-SNR power offset. Moreover, it is also found
that increasing the number of REs can reduce the high-SNR power offset and
thus enhance the ECC.

##  II System Model

In a MIMO channel, the transmitter sends signals to a receiver with the aid of
an RIS. Assume that the transmitter and receiver are equipped with uniform
linear arrays (ULAs) contanining $N_{\text{t}}$ and $N_{\text{r}}$ antennas,
respectively. Besides, the RIS is equipped with a uniform planar array (UPA)
containing $M$ reflecting elements with $M_{\text{h}}$ ones per row and
$M_{\text{v}}$ ones per column. The signal vector at the receiver end is given
by

$\displaystyle{\textbf{y}}=\left(\sqrt{\beta_{\text{t}}\beta_{\text{r}}}{\textbf{R}}{\bm{\Phi}}{\textbf{T}}+\sqrt{\beta_{\text{d}}}{\textbf{H}}\right){\textbf{x}}+{\textbf{n}}={\textbf{G}}{\textbf{x}}+{\textbf{n}},$
(1)

where ${\textbf{x}}\in{\mathbb{C}}^{N_{\text{t}}\times 1}$ is the transmitted
signal vector satisfying
${\mathbbmss{E}}\left\\{{\textbf{x}}^{{\dagger}}{\textbf{x}}\right\\}=P$ with
$P$ being the total transmit power; ${\textbf{T}}\in{\mathbb{C}}^{M\times
N_{\text{t}}}$, ${\textbf{R}}\in{\mathbb{C}}^{N_{\text{r}}\times M}$, and
${\textbf{H}}\in{\mathbb{C}}^{N_{\text{r}}\times N_{\text{t}}}$ denote the
channel matrices from the transmitter to the RIS, from the RIS to the
receiver, and from the transmitter to the receiver, respectively;
${\bm{\Phi}}={\text{diag}}\left\\{{\text{e}}^{\textsf{j}\theta_{1}},{\text{e}}^{\textsf{j}\theta_{2}},\cdots,{\text{e}}^{\textsf{j}\theta_{M}}\right\\}\in{\mathbb{C}}^{M\times
M}$ denotes the diagonal reflection matrix of the RIS with
${\theta_{i}}\in\left(-\pi,\pi\right]$, $\forall i=1,2,\cdots,M$;
${\textbf{n}}\sim{\mathcal{CN}}\left({\textbf{0}},\sigma^{2}{\textbf{I}}_{N_{\text{r}}}\right)$
represents the additive white Gaussian noise (AWGN) with $\sigma^{2}$ being
the noise power at each receive antenna; $\beta_{\text{t}}$,
$\beta_{\text{r}}$, $\beta_{\text{d}}$ model both path loss and shadowing in
the transmitter-RIS, RIS-receiver, and transmitter-receiver channels,
respectively;
${\textbf{G}}=\sqrt{\beta_{\text{t}}\beta_{\text{r}}}{\textbf{R}}{\bm{\Phi}}{\textbf{T}}+\sqrt{\beta_{\text{d}}}{\textbf{H}}\in{\mathbb{C}}^{N_{\text{r}}\times
N_{\text{t}}}$ denotes the effective channel matrix from the transmitter to
the receiver.

As stated before, there are deterministic LoS paths in the RIS-related
channels in general. To highlight their influence and to make the following
analyses more amenable, assume that T is a pure LoS channel [9]. Moreover, the
uncorrelated frequency-flat Rician fading is used to model the propagation
environment of R and H. We comment that the following discussions can be
directly extended to the scenario when R is a LoS channel while T is a Rician
one. By [11], we have

$\displaystyle{\textbf{T}}=\sum\nolimits_{i=1}^{P_{\textbf{T}}}\frac{t_{i}}{\sqrt{P_{\textbf{T}}}}{\textbf{a}}_{M_{\text{h}},M_{\text{v}}}\left(\varphi_{{\text{a}}_{1,i}},\theta_{{\text{a}}_{1,i}}\right){\textbf{a}}_{N_{\text{t}}}\left(\theta_{{\text{d}}_{1,i}}\right)^{\intercal},$
(2)

where $P_{\textbf{T}}$ is the number of LoS paths;
$\frac{t_{i}}{\sqrt{P_{\textbf{T}}}}\in{\mathbb{C}}$ is the complex gain of
the $i$-th path; $t_{i}$ follow the independent and identically distributed
(i.i.d.) complex Gaussian distribution ${\mathcal{CN}}\left(0,1\right)$;

$\begin{split}&{\textbf{a}}_{M_{\text{h}},M_{\text{v}}}\left(\varphi,\theta\right)^{\intercal}=\left[1,{\text{e}}^{\textsf{j}2\pi\frac{d_{\text{v}}}{\lambda}\sin{\theta}},\cdots,{\text{e}}^{\textsf{j}2\pi\left(M_{\text{v}}-1\right)\frac{d_{\text{v}}}{\lambda}\sin{\theta}}\right]\\\
&\otimes\left[1,{\text{e}}^{\textsf{j}2\pi\frac{d_{\text{h}}}{\lambda}\cos\left(\theta\right)\sin\left(\varphi\right)},\cdots,{\text{e}}^{\textsf{j}2\pi\left(M_{\text{h}}-1\right)\frac{d_{\text{h}}}{\lambda}\cos\left(\theta\right)\sin\left(\varphi\right)}\right],\end{split}$
(3)

is the array response of the RIS with $\otimes$ denoting the Kronecker product
of matrices, $\lambda$ denoting the wavelength of the carrier frequency,
$d_{\text{h}}$ and $d_{\text{v}}$ denoting the inner-element space in each row
and column of the RIS, respectively;
${\textbf{a}}_{N_{\text{t}}}\left(\theta\right)=\left[1,{\text{e}}^{\textsf{j}2\pi\frac{d_{\text{t}}}{\lambda}\sin{\theta}},\cdots,{\text{e}}^{\textsf{j}2\pi\left(N_{\text{t}}-1\right)\frac{d_{\text{t}}}{\lambda}\sin{\theta}}\right]^{\intercal}$
is the array response of the ULA at the transmitter with $d_{\text{t}}$
denoting the inner-element space; $\theta_{{\text{a}}_{i,i}}$ and
$\varphi_{{\text{a}}_{i,i}}$ are the elevation and azimuth angles of arrival
(AoAs), respectively; $\theta_{{\text{d}}_{i,i}}$ are the angles of departure
(AoDs). As for R and H, we have
${\textbf{A}}=\sqrt{L_{\textbf{A}}}\bar{\textbf{A}}+\sqrt{N_{\textbf{A}}}\tilde{\textbf{A}}$,
${\textbf{A}}\in\left\\{{\textbf{R}},{\textbf{H}}\right\\}$, where
$L_{\textbf{A}}+N_{\textbf{A}}=1$; $\tilde{\textbf{A}}$ is the non-LoS (NLoS)
component and the elements of $\tilde{\textbf{A}}$ are i.i.d. complex Gaussian
random variables, each with zero mean and unit variance. Like (2),
$\bar{\textbf{R}}$ and $\bar{\textbf{H}}$ can be written, respectively, as

$\displaystyle\bar{\textbf{R}}=\sum\nolimits_{i=1}^{P_{\textbf{R}}}\frac{r_{i}}{\sqrt{P_{\textbf{R}}}}{\textbf{a}}_{N_{\text{r}}}\left(\theta_{{\text{a}}_{2,i}}\right){\textbf{a}}_{M_{\text{h}},M_{\text{v}}}\left(\varphi_{{\text{d}}_{2,i}},\theta_{{\text{d}}_{2,i}}\right)^{\intercal},$
(4)
$\displaystyle\bar{\textbf{H}}=\sum\nolimits_{i=1}^{P_{\textbf{H}}}\frac{h_{i}}{\sqrt{P_{\textbf{H}}}}{\textbf{a}}_{N_{\text{r}}}\left(\theta_{{\text{a}}_{3,i}}\right){\textbf{a}}_{N_{\text{t}}}\left(\theta_{{\text{d}}_{3,i}}\right)^{\intercal},$
(5)

where
${\textbf{a}}_{N_{\text{r}}}\left(\theta\right)=\left[1,{\text{e}}^{\textsf{j}2\pi\frac{d_{\text{r}}}{\lambda}\sin{\theta}},\cdots,{\text{e}}^{\textsf{j}2\pi\left(N_{\text{r}}-1\right)\frac{d_{\text{r}}}{\lambda}\sin{\theta}}\right]^{\intercal}$
is the array response of the ULA at the receiver with $d_{\text{r}}$ denoting
the inner-element space; $P_{\textbf{R}}$ and $P_{\textbf{H}}$ denote numbers
of LoS paths; ${r_{i}}\sim{\mathcal{CN}}\left(0,1\right)$ and
${h_{i}}\sim{\mathcal{CN}}\left(0,1\right)$; $\theta_{{\text{a}}_{2,i}}$ and
$\theta_{{\text{a}}_{3,i}}$ denote AoAs; $\varphi_{{\text{d}}_{2,i}}$,
$\theta_{{\text{d}}_{2,i}}$, and $\theta_{{\text{d}}_{3,i}}$ denote AoDs.

##  III Ergodic Channel Capacity

Assume that the instantaneous channel state information (CSI) of G is unknown
at the transmitter and perfectly known at the receiver. We comment that this
assumption has been widely used in the MIMO literature, especially for the
frequency division duplex (FDD) systems. Then, the ECC can be written as

$\displaystyle
C={\mathbbmss{E}}\left\\{\log_{2}\det\left({{\textbf{I}}_{N_{\text{r}}}+\bar{\rho}{\textbf{G}}{\textbf{G}}^{{\dagger}}}\right)\right\\},$
(6)

where $\bar{\rho}=\frac{\rho}{N_{\text{t}}}$ with $\rho=\frac{P}{\sigma^{2}}$
being the SNR and the expectation is taken with respect to G. Note that the
phase shifts used at the RIS should be computed by either transmitter or
receiver, because the RIS has no computation power. Also, it is challenging to
estimate R or T, since the RIS cannot send or detect pilot signals. Recalling
the previous assumption, it is clear that neither the transmitter nor receiver
has access to the instantaneous CSI of R and T. As a result, if no more other
information of R and T is available at the transceivers, random phase shifts
are the only choices. In contrast to this, if the statistical CSI of H, T, and
R is available at the transmitter or receiver, then the phase shifts may be
designed more properly by utilizing this information. Yet, in both these two
cases, we can treat $\theta_{i}$ ($\forall i$) as some constants when
analyzing the ECC. We now intend to evaluate the ECC of the considered
channel.

###  III-A New Channel Statistics

######  Lemma 1.

The effective channel, G, satisfies

$\displaystyle{\text{vec}}\left({\textbf{G}}^{{\dagger}}\right)\sim{\mathcal{CN}}\left({\text{vec}}\left(\bar{\textbf{G}}^{{\dagger}}\right),{\textbf{I}}_{N_{\text{r}}}\otimes{\bm{\Psi}}\right),$
(7)

where ${\text{vec}}\left(\cdot\right)$ stacks the columns of a matrix into a
single vector,
$\bar{\textbf{G}}=\sqrt{\beta_{\text{t}}\beta_{\text{r}}L_{\textbf{R}}}\bar{\textbf{R}}{\bm{\Phi}}{\textbf{T}}+\sqrt{\beta_{\text{d}}L_{\textbf{H}}}\bar{\textbf{H}}\in{\mathbb{C}}^{N_{\text{r}}\times
N_{\text{t}}}$,
$\bm{\Psi}={\beta_{\text{t}}\beta_{\text{r}}N_{\textbf{R}}}{\textbf{T}}^{{\dagger}}{\textbf{T}}+{\beta_{\text{d}}N_{\textbf{H}}}{\textbf{I}}_{N_{\text{t}}}\in{\mathbb{C}}^{N_{\text{t}}\times
N_{\text{t}}}$, and ${\bm{\Psi}}\succ 0$.

###### Proof:

This lemma can be proven by first stacking G into a single vector and then
performing some basic mathematical manipulations. ∎

It follows from (7) that G can be written as
${\textbf{G}}=\bar{\textbf{G}}+{\textbf{X}}{\bm{\Psi}}^{1/2}$ [12], where
${\textbf{X}}\in{\mathbb{C}}^{N_{\text{r}}\times N_{\text{t}}}$ contains
independent ${\mathcal{CN}}\left(0,1\right)$ elements. This further leads to
the following remark.

######  Remark 1.

The effective channel of the considered RIS-aided MIMO system is equivalent to
a transmit correlated Rician MIMO channel with its mean matrix giving by
$\bar{\textbf{G}}$ and its transmit correlation matrix giving by $\bm{\Psi}$.

Generally, it is challenging to derive a concise as well as compact expression
to characterize the ECC of a MIMO fading channel. As a compromise, we intend
to derive some capacity bounds to evaluate the SE performance of the
considered RIS-aided system and use them to unveil more system insights.

###  III-B Capacity Upper Bound

Let us introduce some key preliminary results that will be useful in
constructing the ECC’s upper bound.

######  Lemma 2.

For ${\textbf{A}}\in{\mathbb{C}}^{n\times n}$, let
${\textbf{A}}_{\mathcal{F}}$ denote the sub-matrix formed from A by removing
columns and rows not indexed by the elements of
${\mathcal{F}}\subseteq\left\\{1,2,\cdots,n\right\\}$. Then, we have

$\displaystyle\det\left({\textbf{I}}_{n}+\lambda{\textbf{A}}\right)=\sum_{t=0}^{n}\lambda^{t}\sum_{{\mathcal{F}}\subseteq\left\\{1,2,\cdots,n\right\\},\left|{\mathcal{F}}\right|=t}\det\left({\textbf{A}}_{{\mathcal{F}}}\right),$
(8)

where ${\mathcal{F}}$ is an ordered subset of $\left\\{1,\cdots,n\right\\}$
with $\left|{\mathcal{F}}\right|=t$ and $\lambda$ is a constant.

###### Proof:

Please refer to [13] for more details. ∎

######  Lemma 3.

Given a random matrix ${\textbf{B}}\in{\mathbb{C}}^{p\times q}$ ($p\leq q$)
satisfying
${\text{vec}}\left({\textbf{B}}^{{\dagger}}\right)\sim{\mathcal{CN}}\left({\text{vec}}\left(\bar{\textbf{B}}^{{\dagger}}\right),{\bm{\Omega}}\otimes{\textbf{I}}_{q}\right)$
with ${\mathbbmss{E}}\left\\{\textbf{B}\right\\}=\bar{\textbf{B}}$ and
$\bm{\Omega}\in{\mathbb{C}}^{p\times p}$ being a positive definite matrix,
then we have

$\displaystyle{\mathbbmss{E}}\left\\{\det\left({\textbf{B}}{\textbf{B}}^{{\dagger}}\right)\right\\}=\frac{\left(q-L\right)!}{\left(q-p\right)!}\frac{{\mathsf{J}}\left(\theta_{1},\cdots,\theta_{L}\right)}{\det\left({\bm{\Omega}}^{-1}\right)},$
(9)

where
${\mathsf{J}}\left(\theta_{1},\cdots,\theta_{L}\right)=\det\left({\bm{\Delta}}\left(\theta_{1},\cdots,\theta_{L}\right)\right)/V\left(\theta_{1},\cdots,\theta_{L}\right)$
$\left\\{\theta_{i}\right\\}_{i=1}^{L}$ are the non-zero eigenvalues of
${\bm{\Omega}}^{-1}\bar{\textbf{B}}\bar{\textbf{B}}^{{\dagger}}$,
$V\left(\theta_{1},\cdots,\theta_{L}\right)=\prod_{i<j\leq
L}\left(\theta_{j}-\theta_{i}\right)$, and the $(i,j)$-th element of
${\bm{\Delta}}\left(\theta_{1},\cdots,\theta_{L}\right)\in{\mathbb{C}}^{L\times
L}$ is given by $\left(q-L+j+\theta_{i}\right)\theta_{i}^{j-1}$.

###### Proof:

Please refer to [14] for more details. ∎

######  Lemma 4.

Given ${\textbf{C}}\in{\mathbb{C}}^{p\times q}$ ($p>q$),
$\det{\left({\textbf{C}}{\textbf{C}}^{{\dagger}}\right)}=0$.

###### Proof:

By using the fact that
${\rm{rank}}\left({\textbf{C}}{\textbf{C}}^{{\dagger}}\right)\leq q<p$, this
lemma can be found. ∎

####  III-B1 $N_{\text{r}}\geq N_{\text{t}}$

When $N_{\text{r}}\geq N_{\text{t}}$, we can rewrite (6) as

$\displaystyle C$
$\displaystyle={\mathbbmss{E}}\left\\{\log_{2}\det\left({{\textbf{I}}_{N_{\text{t}}}+\bar{\rho}{\textbf{G}}^{{\dagger}}{\textbf{G}}}\right)\right\\}$
(10)
$\displaystyle\leq\log_{2}{{\mathbbmss{E}}\left\\{\det\left({{\textbf{I}}_{N_{\text{t}}}+\bar{\rho}{\textbf{G}}^{{\dagger}}{\textbf{G}}}\right)\right\\}},$
(11)

where (10) is based on the Sylvester’s determinant identity and (11) is based
on the concavity of logarithmic functions. Denoting
${\bm{\Sigma}}=\bar{\textbf{G}}^{{\dagger}}\bar{\textbf{G}}$, the following
theorem can be found.

######  Theorem 1.

When $N_{\text{r}}\geq N_{\text{t}}$, the ECC of the RIS-aided MIMO channel is
upper bounded by

$\displaystyle C$
$\displaystyle\leq\log_{2}\left(\sum\nolimits_{t=0}^{N_{\text{t}}}\bar{\rho}^{t}\sum\nolimits_{{\mathcal{Q}}\subseteq\left\\{1,\cdots,N_{\text{t}}\right\\},\left|{\mathcal{Q}}\right|=t}\det\left({\bm{\Psi}}_{{\mathcal{Q}}}\right)\right.$
$\displaystyle\times\left.{\mathsf{J}}\left(\tilde{\theta}_{1},\cdots,\tilde{\theta}_{L_{{\mathcal{Q}}}}\right){\left(N_{\text{r}}-L_{{\mathcal{Q}}}\right)!}/{\left(N_{\text{r}}-t\right)!}\right)\triangleq
C_{1},$ (12)

where ${\mathcal{Q}}$ is an ordered subset of
$\left\\{1,\cdots,N_{\text{t}}\right\\}$ with $\left|{\mathcal{Q}}\right|=t$
and $\left\\{\tilde{\theta}_{i}\right\\}_{i=1}^{L_{{\mathcal{Q}}}}$ are the
non-zero eigenvalues of
${\bm{\Psi}}_{{\mathcal{Q}}}^{-1}{\bm{\Sigma}}_{{\mathcal{Q}}}$.

###### Proof:

This theorem can be proved by firstly expanding the expectation in (11) based
on (8) and then calculating the resultant expectations based on (9). ∎

####  III-B2 $N_{\text{t}}>N_{\text{r}}$

When $N_{\text{t}}>N_{\text{r}}$, it can be deduced from Lemma 4 that the
determinants of some ${\bm{\Sigma}}={\textbf{G}}^{{\dagger}}{\textbf{G}}$
higher order sub-matrices equal zero, which further yields

$\displaystyle\det\left({{\textbf{I}}_{N_{\text{t}}}+\bar{\rho}{\bm{\Sigma}}}\right)=\sum\nolimits_{t=0}^{N_{\text{r}}}\bar{\rho}^{t}\sum\nolimits_{{\mathcal{Q}}_{i}}\det\left({\bm{\Sigma}}_{{\mathcal{Q}}_{i}}\right).$
(13)

Here, ${\mathcal{Q}}$ is an ordered subset of
$\left\\{1,\cdots,N_{\text{t}}\right\\}$ with $\left|{\mathcal{Q}}\right|=t$.
Leveraging (9), (11), and (13), the capacity upper bound for
$N_{\text{t}}>N_{\text{r}}$ can be calculated and the main result is
summarized in the following theorem.

######  Theorem 2.

When $N_{\text{t}}>N_{\text{r}}$, the ECC of the RIS-aided MIMO channel is
upper bounded by

$\begin{split}C&\leq\log_{2}\left(\sum\nolimits_{t=0}^{N_{\text{r}}}\bar{\rho}^{t}\sum\nolimits_{{\mathcal{Q}}\subseteq\left\\{1,\cdots,N_{\text{t}}\right\\},\left|{\mathcal{Q}}_{i}\right|=t}\det\left({\bm{\Psi}}_{{\mathcal{Q}}}\right)\right.\\\
&\times\left.{\mathsf{J}}\left(\tilde{\theta}_{1},\cdots,\tilde{\theta}_{L_{{\mathcal{Q}}}}\right){\left(N_{\text{r}}-L_{{\mathcal{Q}}}\right)!}/{\left(N_{\text{r}}-t\right)!}\right)\triangleq
C_{2}.\end{split}$ (14)

###  III-C Capacity Lower Bound

Turning now the capacity lower bound, we establish the following preliminary
results before the derivations,.

######  Lemma 5.

Given a random matrix ${\textbf{B}}\in{\mathbb{C}}^{p\times q}$ ($p\leq q$)
satisfying
${\text{vec}}\left({\textbf{B}}^{{\dagger}}\right)\sim{\mathcal{CN}}\left({\text{vec}}\left(\bar{\textbf{B}}^{{\dagger}}\right),{\textbf{I}}_{p}\otimes{\textbf{I}}_{q}\right)$
with ${\mathbbmss{E}}\left\\{\textbf{B}\right\\}=\bar{\textbf{B}}$. Then, we
have

$\begin{split}{\mathbbmss{E}}\left\\{\ln\det\left({\textbf{B}}{\textbf{B}}^{{\dagger}}\right)\right\\}=\sum_{k=0}^{p-1}\psi\left(q-k\right)+{\mathsf{F}}\left(\theta_{1},\cdots,\theta_{L}\right),\end{split}$
(15)

where $\psi\left(x\right)=\frac{{\rm d}}{{\rm d}x}\ln{\Gamma\left(x\right)}$
is the Digamma function [10, eq. (6.461)],
$\Gamma\left(x\right)=\int_{0}^{\infty}t^{x-1}{\text{e}}^{-t}{\rm d}t$ is the
gamma function [10, eq. (6.1.1)], $\left\\{\theta_{i}\right\\}_{i=1}^{L}$ are
the non-zero eigenvalues of $\bar{\textbf{B}}\bar{\textbf{B}}^{{\dagger}}$,

$\displaystyle{\mathsf{F}}\left(\theta_{1},\cdots,\theta_{L}\right)=\frac{\sum_{i=1}^{L}{\det{\left(\bar{\bm{\Delta}}_{i}\left(\theta_{1},\cdots,\theta_{L}\right)\right)}}}{V\left(\theta_{1},\cdots,\theta_{L}\right)},$
(16)
$\displaystyle\bar{\bm{\Delta}}_{i}\left(\theta_{1},\cdots,\theta_{L}\right)$
$\displaystyle=\left[{\textbf{g}}\left(\theta_{1}\right),\cdots,{\textbf{g}}\left(\theta_{i-1}\right),{\textbf{h}}\left(\theta_{i}\right),{\textbf{g}}\left(\theta_{i+1}\right),\cdots,{\textbf{g}}\left(\theta_{L}\right)\right],$
(17)
$\displaystyle{\textbf{g}}\left(x\right)=\left[1,x,x^{2},\cdots,x^{L-1}\right]^{\intercal},$
(18)
$\displaystyle{\textbf{h}}\left(x\right)=\left[h_{1}\left(x\right),h_{2}\left(x\right),\cdots,h_{L}\left(x\right)\right]^{\intercal},$
(19) $\displaystyle
h_{i}\left(x\right)=x^{i-1}\sum\nolimits_{k=0}^{\infty}\frac{1-\exp\left(-x\right)\sum_{n=0}^{k}\frac{x^{n}}{n!}}{q-L+i+k}.$
(20)

###### Proof:

Please refer to [14] for more details. ∎

######  Lemma 6.

Let ${\textbf{Z}}\in{\mathbb{C}}^{q\times q}$ be a positive definite matrix
whose eigenvalues are $\lambda_{1}\geq\cdots\geq\lambda_{q}>0$. Then for
arbitrary matrix ${\textbf{X}}\in{\mathbb{C}}^{p\times q}$ ($p\leq q$),

$\begin{split}\det\left({\textbf{X}}{\textbf{X}}^{{\dagger}}\right)\prod_{i=1}^{p}\lambda_{i}\geq\det\left({\textbf{X}}{\textbf{Z}}{\textbf{X}}^{{\dagger}}\right)\geq\det\left({\textbf{X}}{\textbf{X}}^{{\dagger}}\right)\prod_{i=1}^{p}\lambda_{q-i+1}.\end{split}$
(21)

###### Proof:

Please refer to [15] for more details. ∎

####  III-C1 $N_{\text{r}}\geq N_{\text{t}}$

By applying the Minkowski’s determinant inequality [16, eq. (7.8.22)] to (10)
and then leveraging the convexity of
$\log_{2}\left(1+a\exp\left(x\right)\right)$ in $x$ for $a>0$, we can obtain

$\displaystyle C$ $\displaystyle\geq
N_{\text{t}}{\mathbbmss{E}}\left\\{\log_{2}\left(1+\bar{\rho}\exp\left({N_{\text{t}}^{-1}}\ln\det\left({\textbf{G}}^{{\dagger}}{\textbf{G}}\right)\right)\right)\right\\}$
(22) $\displaystyle\geq
N_{\text{t}}\log_{2}\left(1+\bar{\rho}\exp\left({N_{\text{t}}^{-1}}{\mathbbmss{E}}\left\\{\ln\det\left({\textbf{G}}^{{\dagger}}{\textbf{G}}\right)\right\\}\right)\right).$
(23)

As stated before,
${\textbf{G}}=\bar{\textbf{G}}+{\textbf{X}}{\bm{\Psi}}^{1/2}$, which yields

$\displaystyle\ln\det\left({\textbf{G}}^{{\dagger}}{\textbf{G}}\right)=\ln\det\left(\bm{\Psi}\right)+\ln\det\left(\tilde{\textbf{G}}\tilde{\textbf{G}}^{{\dagger}}\right),$
(24)

where
${\text{vec}}\left(\tilde{\textbf{G}}^{{\dagger}}\right)\sim{\mathcal{CN}}\left({\text{vec}}\left(\bar{\textbf{G}}{\bm{\Psi}}^{-1/2}\right),{\textbf{I}}_{N_{\text{t}}}\otimes{\textbf{I}}_{N_{\text{r}}}\right)$.
By substituting (24) into (23) and calculating the resultant expectation by
(15), we arrive at the following theorem.

######  Theorem 3.

When $N_{\text{r}}\geq N_{\text{t}}$, the ECC of the RIS-aided MIMO channel is
lower bounded by

$\begin{split}&C\geq
N_{\text{t}}\log_{2}\left(1+\bar{\rho}\exp\left({N_{\text{t}}^{-1}}\left(\ln\det\left(\bm{\Psi}\right)\right.\right.\right.\\\
&\left.\left.\left.+\sum\nolimits_{k=0}^{N_{\text{t}}-1}\psi\left(N_{\text{r}}-k\right)+{\mathsf{F}}\left(\alpha_{1},\cdots,\alpha_{\tilde{L}}\right)\right)\right)\right)\triangleq
C_{3},\end{split}$ (25)

where $\left\\{\alpha_{i}\right\\}_{i=1}^{\tilde{L}}$ are the non-zero
eigenvalues of
${\bm{\Psi}}^{-\frac{1}{2}}\bar{\textbf{G}}^{{\dagger}}\bar{\textbf{G}}{\bm{\Psi}}^{-\frac{1}{2}}$.

####  III-C2 $N_{\text{t}}>N_{\text{r}}$

By following the approach in obtaining (23), we can obtain

$C\geq
N_{\text{r}}\log_{2}\left(1+\bar{\rho}\exp\left({N_{\text{r}}^{-1}}{\mathbbmss{E}}\left\\{\ln\det\left({\textbf{G}}{\textbf{G}}^{{\dagger}}\right)\right\\}\right)\right).$
(26)

It follows from ${\textbf{G}}=\bar{\textbf{G}}+{\textbf{X}}{\bm{\Psi}}^{1/2}$
that
$\ln\det\left({\textbf{G}}{\textbf{G}}^{{\dagger}}\right)=\ln\det\left(\tilde{\textbf{G}}^{{\dagger}}{\bm{\Psi}}\tilde{\textbf{G}}\right)$,
which together with (21), yields

$\ln\det\left({\textbf{G}}{\textbf{G}}^{{\dagger}}\right)\geq\ln\det\left(\tilde{\textbf{G}}^{{\dagger}}\tilde{\textbf{G}}\right)+\sum\nolimits_{i=1}^{N_{\text{r}}}\ln{\zeta_{N_{\text{t}}-i+1}},$
(27)

where $\zeta_{1}\geq\cdots\geq\zeta_{N_{\text{t}}}>0$ are the eigenvalues of
$\bm{\Psi}$. Note that
${\text{vec}}\left(\tilde{\textbf{G}}\right)\sim{\mathcal{CN}}\left({\text{vec}}\left({\bm{\Psi}}^{-1/2}\bar{\textbf{G}}^{{\dagger}}\right),{\textbf{I}}_{N_{\text{r}}}\otimes{\textbf{I}}_{N_{\text{t}}}\right)$.
By substituting (27) into (26) and calculating the resultant expectation based
on (15), we arrive at the following theorem.

######  Theorem 4.

When $N_{\text{t}}>N_{\text{r}}$, the ECC of the RIS-aided MIMO channel is
lower bounded by

$\begin{split}&C\geq
N_{\text{r}}\log_{2}\left(1+\bar{\rho}\exp\left({N_{\text{r}}^{-1}}\left(\sum\nolimits_{i=1}^{N_{\text{r}}}\ln{\zeta_{N_{\text{t}}-i+1}}\right.\right.\right.\\\
&\left.\left.\left.+\sum\nolimits_{k=0}^{N_{\text{r}}-1}\psi\left(N_{\text{t}}-k\right)+{\mathsf{F}}\left(\tilde{\alpha}_{1},\cdots,\tilde{\alpha}_{\bar{L}}\right)\right)\right)\right)\triangleq
C_{4},\end{split}$ (28)

where $\left\\{\tilde{\alpha}_{i}\right\\}_{i=1}^{\bar{L}}$ are the non-zero
eigenvalues of $\bar{\textbf{G}}{\bm{\Psi}}^{-1}\bar{\textbf{G}}^{{\dagger}}$.

######  Remark 2.

When the statistical CSI including $\bar{\textbf{G}}$ and $\bm{\Psi}$ is
available at the transmitter or receiver, Theorems 1, 2, 3, and 4 can provide
insightful guidelines for designing the reflecting matrix to improve the ECC
of the considered system.

More specifically, since the ECC defined in (6) lacks any closed-form
expressions, it makes sense to solve the problems given by
${\text{P}}_{i}:{\bm{\Phi}}^{\star}=\arg\max\nolimits_{\bm{\Phi}}\left\\{C_{i}\right\\}$
($i\in\left\\{1,2,3,4\right\\}$) for the sake of finding a suboptimal
$\bm{\Phi}$ in terms of maximizing the ECC. Generally, problems
(${\text{P}}_{i}$) can be solved by using the classical genetic algorithm (GA)
[17]. Yet, since this letter focuses more on the theoretical analyses of the
ECC instead of the design of $\bm{\Phi}$, we do not provide more discussions
about the design of $\bm{\Phi}$.

##  IV Asymptotic Analyses

We now intend to apply some asymptotic analyses to the previously derived
bounds to unveil more system insights. Specifically, two limiting cases
including $\rho\rightarrow\infty$ and $M\rightarrow\infty$ will be considered.

###  IV-A Asymptotic Analyses in the High-SNR Regime

####  IV-A1 $N_{\text{r}}\geq N_{\text{t}}$

It can be derived from (10) that

$\displaystyle\lim\nolimits_{\rho\rightarrow\infty}\left(C-N_{\text{t}}\log_{2}\bar{\rho}\right)=\frac{1}{\ln{2}}{\mathbbmss{E}}\left\\{\ln\det\left({\textbf{G}}^{{\dagger}}{\textbf{G}}\right)\right\\}.$
(29)

Interestingly, it can be found that (29) can be also derived from (23) by
setting $\rho$ as infinity, which suggests that this lower bound is exact in
the high-SNR regime. The ECC in the high-SNR regime can be generally
approximated as [18]

$C\approx{\mathcal{S}}_{\infty}\left(\log_{2}{\rho}-{\mathcal{L}}_{\infty}\right),$
(30)

where
${\mathcal{S}}_{\infty}=\lim_{\rho\rightarrow\infty}\frac{C}{\log_{2}{\rho}}$
denotes the high-SNR slope in bits/s/Hz/(3 dB), and
${\mathcal{L}}_{\infty}=\lim_{\rho\rightarrow\infty}\left(\log_{2}{\rho}-\frac{C}{{\mathcal{S}}_{\infty}}\right)$
denotes the high-SNR power offset in 3 dB units. Based on (15) and (29), we
can evaluate $S_{\infty}$ and $L_{\infty}$ in closed-form as follows.

######  Theorem 5.

When $N_{\text{r}}\geq N_{\text{t}}$, the high-SNR slope and high-SNR power
offset of the RIS-aided MIMO channel are, respectively, given by

$\displaystyle{\mathcal{S}}_{\infty}=$ $\displaystyle N_{\text{t}},$ (31)
$\displaystyle{\mathcal{L}}_{\infty}=$
$\displaystyle\log_{2}{N_{\text{t}}}-\left({N_{\text{t}}\ln{2}}\right)^{-1}\left(\ln\det\left(\bm{\Psi}\right)+\right.$
$\displaystyle\left.+\sum\nolimits_{k=0}^{N_{\text{t}}-1}\psi\left(N_{\text{r}}-k\right)+{\mathsf{F}}\left(\alpha_{1},\cdots,\alpha_{\tilde{L}}\right)\right).$
(32)

######  Remark 3.

The results shown in (31) suggest that the high-SNR slope of the ECC curve
over logarithmic SNR is $N_{\text{t}}$. In other words, when $N_{\text{r}}\geq
N_{\text{t}}$, the ECC increases by approximately $N_{\text{t}}$ bits/s/Hz for
each 3 dB increase in SNR.

It is clear that $N_{\text{t}}$ is independent of $\bm{\Phi}$ while
$\left\\{\alpha_{i}\right\\}_{i=1}^{\tilde{L}}$ are functions of $\bm{\Phi}$,
which leads to the following remark.

######  Remark 4.

Equations (31) and (32) imply that the RIS’s deployment can shape the ECC
curve over SNR by controlling its high-SNR power offset. Moreover, the RIS’s
deployment cannot influence the high-SNR slope.

####  IV-A2 $N_{\text{t}}>N_{\text{r}}$

By following the approach in obtaining (31) and (32), we find that when
$N_{\text{t}}>N_{\text{r}}$, the high-SNR slope and high-SNR power offset can
be written, respectively, as

$\displaystyle{\mathcal{S}}_{\infty}$ $\displaystyle=N_{\text{r}},$ (33)
$\displaystyle{\mathcal{L}}_{\infty}$
$\displaystyle=\log_{2}{N_{\text{t}}}-\frac{1}{N_{\text{r}}\ln{2}}{\mathbbmss{E}}\left\\{\ln{\det{\left(\tilde{\textbf{G}}^{{\dagger}}{\bm{\Psi}}\tilde{\textbf{G}}\right)}}\right\\}.$
(34)

Equation (33) implies that when $N_{\text{t}}>N_{\text{r}}$, the ECC increases
by approximately $N_{\text{r}}$ bits/s/Hz for each 3 dB increase in SNR.
Moreover, by leveraging Lemma 5 and Lemma 6 to calculate the expectation in
(34), we can obtain the following theorem.

######  Theorem 6.

When $N_{\text{t}}>N_{\text{r}}$, the high-SNR power offset of the RIS-aided
MIMO channel is upper bounded by

$\displaystyle{\mathcal{L}}_{\infty}\leq\log_{2}{N_{\text{t}}}-\frac{\sum_{i=1}^{N_{\text{r}}}\ln{\zeta_{N_{\text{t}}-i+1}}+{\mathsf{L}}\left(\tilde{\alpha}_{1},\cdots,\tilde{\alpha}_{\bar{L}}\right)}{N_{\text{r}}\ln{2}}$
(35)

and lower bounded by

$\displaystyle{\mathcal{L}}_{\infty}\geq\log_{2}{N_{\text{t}}}-\frac{\sum_{i=1}^{N_{\text{r}}}\ln{\zeta_{i}}+{\mathsf{L}}\left(\tilde{\alpha}_{1},\cdots,\tilde{\alpha}_{\bar{L}}\right)}{N_{\text{r}}\ln{2}},$
(36)

where
${\mathsf{L}}\left(\tilde{\alpha}_{1},\cdots,\tilde{\alpha}_{\bar{L}}\right)=\sum_{k=0}^{N_{\text{r}}-1}\psi\left(N_{\text{t}}-k\right)+{\mathsf{F}}\left(\tilde{\alpha}_{1},\cdots,\tilde{\alpha}_{\bar{L}}\right)$.

###  IV-B Asymptotic Analyses in the Limit of Many REs

We first rewrite (6) as

$\displaystyle
C={\mathbbmss{E}}\left\\{\log_{2}\det\left({{\textbf{I}}_{N_{\text{r}}}+\bar{\rho}M\left(\frac{{\textbf{G}}}{\sqrt{M}}\right)\left(\frac{{\textbf{G}}}{\sqrt{M}}\right)^{{\dagger}}}\right)\right\\}.$
(37)

Since
${\text{vec}}\left({\textbf{G}}^{{\dagger}}\right)\sim{\mathcal{CN}}\left({\text{vec}}\left(\bar{\textbf{G}}^{{\dagger}}\right),{\textbf{I}}_{N_{\text{r}}}\otimes{\bm{\Psi}}\right)$,
we can obtain
${\text{vec}}\left(\frac{{\textbf{G}}^{{\dagger}}}{\sqrt{M}}\right)\sim{\mathcal{CN}}\left({\text{vec}}\left(\frac{\bar{\textbf{G}}^{{\dagger}}}{\sqrt{M}}\right),{\textbf{I}}_{N_{\text{r}}}\otimes\frac{{\bm{\Psi}}}{M}\right)$.
Based on [6, Page 9],
$\left|{\textbf{a}}_{M_{\text{h}},M_{\text{v}}}\left(\varphi_{{\text{d}}_{2,i}},\theta_{{\text{d}}_{2,i}}\right)^{\intercal}{\bm{\Phi}}{\textbf{a}}_{M_{\text{h}},M_{\text{v}}}\left(\varphi_{{\text{a}}_{1,j}},\theta_{{\text{a}}_{1,j}}\right)\right|$
($\forall i,j$) is a bounded value under random phase shifts, which together
with (2), (4), and (5), yields
$\lim_{M\rightarrow\infty}\frac{\bar{\textbf{G}}^{{\dagger}}}{\sqrt{M}}={{\textbf{0}}_{N_{\text{t}}\times
N_{\text{r}}}}$. Besides,

$\lim_{M\rightarrow\infty}\frac{{\bm{\Psi}}}{M}={\beta_{\text{t}}\beta_{\text{r}}N_{\textbf{R}}}\sum_{i=1}^{P_{\textbf{T}}}\frac{\left|t_{i}\right|^{2}}{P_{\textbf{T}}}{\textbf{a}}_{N_{\text{t}}}\left(\theta_{{\text{d}}_{1,i}}\right)^{\ast}{\textbf{a}}_{N_{\text{t}}}\left(\theta_{{\text{d}}_{1,i}}\right)^{\intercal}\triangleq{\bm{\Upsilon}}.$

Thus,
$\lim_{M\rightarrow\infty}{\text{vec}}\left(\frac{{\textbf{G}}^{{\dagger}}}{\sqrt{M}}\right)\sim{\mathcal{CN}}\left({\textbf{0}},{\textbf{I}}_{N_{\text{r}}}\otimes{\bm{\Upsilon}}\right)$.
It follows that as $M\rightarrow\infty$, the ECC satisfies
$C-\dot{C}\left(\rho\right)\xrightarrow[M\rightarrow\infty]{a.s.}0$, where

$\displaystyle\dot{C}\left(\rho\right)={\mathbbmss{E}}_{\textbf{X}}\left\\{\log_{2}\det\left({\textbf{I}}_{N_{\text{r}}}+M{\rho}/{N_{\text{t}}}{\textbf{X}}{\bm{\Upsilon}}{\textbf{X}}^{{\dagger}}\right)\right\\}$
(38)

with ${\textbf{X}}\in{\mathbb{C}}^{N_{\text{r}}\times N_{\text{t}}}$
containing independent ${\mathcal{CN}}\left(0,1\right)$ elements. We now
intend to consider the case where the transmit power is scaled down by a
factor of $1/M^{\alpha}$. From (38), we have
$\lim_{M\rightarrow\infty}\dot{C}\left({\rho}{M^{-\alpha}}\right)=\infty$ for
$\alpha<1$,
$\lim_{M\rightarrow\infty}\dot{C}\left({\rho}{M^{-\alpha}}\right)=0$ for
$\alpha>1$, and
$\lim_{M\rightarrow\infty}\dot{C}\left({\rho}/{M}\right)={\mathbbmss{E}}_{\textbf{X}}\left\\{\log_{2}\det\left({\textbf{I}}_{N_{\text{r}}}+\rho/N_{\text{t}}{\textbf{X}}{\bm{\Upsilon}}{\textbf{X}}^{{\dagger}}\right)\right\\}$.

######  Remark 5.

It is clear that $\dot{C}\left({\rho}{M^{-\alpha}}\right)$ is a non-zero
constant if and only if $\alpha=1$, which suggests that in the large limit of
$M$, with no degradation in the ECC, the total SNR or transmit power can be
cut down at most ${1}/{M}$.

It is challenging to derive a close-form expression for
$\tilde{C}=\dot{C}\left(\rho/M\right)$. Fortunately, by following similar
steps outlined in Section III-B, a closed-form upper bound of $\tilde{C}$ is
available. For the sake of brevity, we omit the details here.

####  IV-B1 $N_{\text{r}}\geq N_{\text{t}}$

With the aid of (12) and (25), the capacity bounds for $N_{\text{r}}\geq
N_{\text{t}}$ can be directly obtained. Furthermore, the following corollary
can be found.

######  Corollary 1.

In the large limit of $M$, the high-SNR slope and high-SNR power offset for
the case of $N_{\text{r}}\geq N_{\text{t}}$ are given by
${\mathcal{S}}_{\infty}=N_{\text{t}}$ and
${\mathcal{L}}_{\infty}=\log_{2}\frac{N_{\text{t}}}{M}-\frac{1}{N_{\text{t}}\ln{2}}\left(\ln\det\left(\bm{\Upsilon}\right)+\sum_{k=0}^{N_{\text{t}}-1}\psi\left(N_{\text{r}}-k\right)\right)$,
respectively.

We comment that in obtaining Corollary 1, we have assumed that
$\left|\ln\det\left(\bm{\Upsilon}\right)\right|<\infty$ for simplicity. The
results in Corollary 1 suggest that increasing REs’ number can decrease the
high-SNR power offset and thus enhance the ECC.

####  IV-B2 $N_{\text{t}}>N_{\text{r}}$

By following the approach in the last section, the capacity bounds under the
scenario of $N_{\text{t}}>N_{\text{r}}$ can be also discussed. Thus, we do not
provide further detailed statements.

##  V Numerical Results

(a) $N_{\text{r}}\geq N_{\text{t}}$.

(b) $N_{\text{t}}>N_{\text{r}}$.

Figure 1: ECC Bounds of RIS-aided MIMO channels with $M=16$. The asymptotic
results in Figure 1(a) are calculated from (30), (31), and (32), whereas the
asymptotic results in Figure 1(b) are calculated from (30), (33), and (35).

(a) $\rho=10$ dB.

(b) $\rho=E/M$ and $E=10$ dB.

Figure 2: The ECC versus the number of REs with $M_{\text{h}}=4$. The
asymptotic ECC ($\tilde{C}$) in Figure 2(b) is obtained by simulation, whereas
the asymptotic upper bound is obtained by following similar steps in deriving
(12) and (14).

Numerical results are provided in this part to verify the previous analyses.
During the simulation, we set
$d_{\text{t}}=d_{\text{r}}=d_{\text{v}}=d_{\text{h}}=\frac{\lambda}{2}$,
$\beta_{\text{t}}=\beta_{\text{r}}=\beta_{\text{d}}=1$,
$P_{\textbf{R}}=P_{\textbf{T}}=1$, $P_{\textbf{H}}=2$, and
$L_{\textbf{R}}=L_{\textbf{H}}=\frac{2}{3}$. Moreover, angles of the LoS
components in the assistant channel is randomly set within
$\left[0,2\pi\right)$.

Figure 1 plots the ECC of the RIS-aided MIMO systems under random phase shifts
for $M=16$ with $M_{\text{h}}=4$. The case of $N_{\text{r}}\geq N_{\text{t}}$
is presented in Figure 1(a). As can be seen, the derived upper bound is tight
for all SNR ranges, which, thus, serves as a good performance metric.
Moreover, the lower bound is tight in the high-SNR regime and the asymptotic
results can track the simulations exactly in the high-SNR region, which is
consistent with the former analyses. Figure 1(b) plots the ECC for
$N_{\text{r}}<N_{\text{t}}$. It can be seen from this graph that the derived
upper bound is still tight for all SNR ranges. Yet, the derived capacity lower
bound is relatively looser than the upper one. Taken this two figures
together, it can be concluded that the derived upper bound can serve as a
tight bound to benchmark the performance of RIS-aided MIMO systems.

Since the derived upper bounds are tight, it makes sense to leverage them to
optimize the phase shifts. Figure 2(a) compares the ECC achieved by optimized,
random, and no phase shifts in terms of the number of REs with $\rho=10$ dB.
Particularly, the phase shifts in Figure 2(a) are optimized by applying the
genetic algorithm to the upper bounds ($C_{1}$ and $C_{2}$). As shown, the
optimized ECC is much higher than the ECC achieved by random phase shifts,
which highlights the importance of phase shifts design. Moreover, it can be
seen that the ECC can get improved by increasing the number of REs. Figure
2(b) presents the ECC with scaled-down power $\rho={E}/{M}$, where $E=10$ dB.
For simplicity, we only provide the results under random phase shifts. As
shown, when $M\rightarrow\infty$, the ECC converges to its asymptotic value,
namely $\tilde{C}$, which, thus, verifies our previous derivations.

##  VI Conclusion

Novel expressions were derived to characterize the upper and lower bounds of
the ECC of RIS-aided MIMO channels with deterministic LoS paths. To unveil
more system insights, asymptotic analyses were performed to the system ECC in
the large limit of the SNR and number of REs. It was found that the RIS can
shape the ECC curve via influencing its high-SNR power offset and the ECC
increases with the number of REs.

## References

  * [1] J. Zhang, E. Björnson, M. Matthaiou, D. W. K. Ng, H. Yang, and D. J. Love, “Prospective multiple antenna technologies for beyond 5G,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 8, pp. 1637–1660, Aug. 2020.
  * [2] Q. Wu and R. Zhang, “Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network,” _IEEE Commun. Mag._ , vol. 58, no. 1. pp 106–112, Jan. 2020.
  * [3] Y. Liu, _et al._ , “Reconfigurable intelligent surfaces: Principles and opportunities,” _IEEE Commun. Surveys Tuts._ , Early Access, 2021.
  * [4] G. Yu, _et al._ , “Design, analysis, and optimization of a large intelligent reflecting surface-aided B5G cellular internet of things,” _IEEE Internet Things J._ , vol. 7, no. 9, pp. 8902–8916, 2020.
  * [5] Y. Han, _et al._ , “Large intelligent surface-assisted wireless communication exploiting statistical CSI,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 8, pp. 8238–8242, Aug. 2019
  * [6] K. Zhi, _et al._ , “Power scaling law analysis and phase shift optimization of RIS-aided massive MIMO systems with statistical CSI,” 2020. [Online]. Available: https://arxiv.org/abs/2010.13525
  * [7] J. Wang, “Joint transmit beamforming and phase shift design for reconfigurable intelligent surface assisted MIMO systems,” _IEEE Trans. Cogn. Commun. Netw._ , Early Access, 2021.
  * [8] J. Zhang, _et al._ , “Large system achievable rate analysis of RIS-assisted MIMO wireless communication with statistical CSIT,” _IEEE Trans. Wireless Commun._ , Early Access, 2021.
  * [9] Q. Nadeem, _et al._ , “Asymptotic max–min SINR analysis of reconfigurable intelligent surface assisted MISO systems,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 12, pp. 7748–7764, Dec. 2020
  * [10] I. S. Gradshteyn and I. M. Ryzhik, _Table of Integrals, Series, and Products_ , 7th ed. Boston, MA, USA: Academic, 2007.
  * [11] H. Bölcskei, _et al._ , “Impact of the propagation environment on the performance of space-frequency coded MIMO-OFDM,” _IEEE J. Sel. Areas Commun._ , vol. 21, no. 3, pp. 427–439, Apr. 2003.
  * [12] A. K. Gupta and D. K. Nagar, _Matrix Variate Distributions_. Boca Raton, FL: Chapman and Hall/CRC, 2000.
  * [13] A. C. Aitken, _Determinants and Matrices_ , 9th ed. Edinburgh, U.K.: Oliver and Boyd, 1956.
  * [14] M. R. McKay and I. B. Collings, “General capacity bounds for spatially correlated Rician MIMO channels,” _IEEE Trans. Inf. Theory_ , vol. 51, no. 9, pp. 3121–3145, Sep. 2005.
  * [15] A. W. Marshall and I. Olkin, _Inequalities: Theory of Majorization and Its Applications_. New York: Academic, 1979.
  * [16] R. Horn and C. R. Johnson, _Matrix Analysis_ , 2nd ed.. New York, NY, USA: Cambridge Univ. Press, 2013.
  * [17] M. Mitchell, _An Introduction to Genetic Algorithms_ , Cambridge, Massachusetts, USA: MIT press, 1998.
  * [18] A. Lozano, _et al._ , “High-SNR power offset in multiantenna communication,” _IEEE Trans. Inf. Theory_ , vol. 51, pp. 4134–4151, Dec. 2005.