# Exact formula for cubic partitions

Lukas Mauth  Department of Mathematics and Computer Science
Division of Mathematics
University of Cologne
Weyertal 86-90
50931 Cologne
Germany<EMAIL_ADDRESS>
###### Abstract.

We obtain an exact formula for the cubic partition function and prove a
conjecture by Banerjee, Paule, Radu and Zeng.

###### Key words and phrases:

Circle Method, $\eta$-function, partitions

##  1\. Introduction and statement of results

Let $n$ be a non-negative integer. A non-increasing finite sequence of
positive integers that sums to $n$ is called a partition of $n.$ We denote by
$p(n)$ the number of partitions of $n,$ which can be defined by the
coefficients of the following $q$-series [1]

$\sum_{n=0}^{\infty}p(n)q^{n}=\prod_{n=1}^{\infty}\frac{1}{1-q^{n}}$ (1.1)

The function $p(n)$ is one of the most well-studied functions troughout number
theory and satisfies numerous remarkable identities. One example is the
following identity discovered by S. Ramanujan [1], [10].

$\sum_{n=0}^{\infty}p(5n+4)q^{n}=5\prod_{n=1}^{\infty}\frac{(1-q^{5n})^{5}}{(1-q^{n})^{6}}.$

An immediate consequence is one of Ramanujan’s celebrated congruences which
asserts that for any non-negative integer

$p(5n+4)\equiv 0\,\,({\rm mod}\,\,{5})$

Similar to (1.1) we define for non-negative integers the cubic partition
function $a(n)$ by

$\sum_{n=0}^{\infty}a(n)q^{n}=\prod_{n=1}^{\infty}\frac{1}{(1-q^{n})(1-q^{2n})}.$
(1.2)

H. Chan proved that the function $a(n)$ satisfies a identity similar to the
one by Ramanujan given above

$\sum_{n=0}^{\infty}a(3n+2)q^{n}=3\prod_{n=1}^{\infty}\frac{(1-q^{3})^{3}(1-q^{6})^{3}}{(1-q)^{4}(1-q^{2})^{4}},$

$a(3n+2)\equiv 0\,\,({\rm mod}\,\,{3}),$

see [6]. His result is closely related to Ramanujan’s continued cubic
fraction, see [7]. For an introduction to Ramanujan’s cubic continued
fraction, see [2],[4]. Furthermore, from the generating function (1.2) it is
obvious that $a(n)$ is the number of partition pairs $(\lambda,\mu)$ such that
$|\lambda|+|\mu|=n$ and $\mu$ consists of only even numbers. These connections
gave rise to the term ”cubic partitions”.

The first $30$ values of $a(n)$ are:
$1,1,3,4,9,12,23,31,54,73,118,159,246,329,489,651,940,1242,1751,$
$2298,3177,4142,5630,7293,9776,12584,16659,21320,27922,35532,46092.$

There is another combinatorial interpretation of the cubic partition function.
Suppose you have a real polynomial equation of degree $n.$ Then, $a(n)$ counts
the number of possibilities for the roots to be real and unequal, real and
equal (in various combinations), or simple or multiple complex conjugates. For
example, for a real cubic polynomial, one can have $3$ real distinct roots,
$3$ real roots and two are equal, $3$ real roots and $1$ real root with a
complex conjugate pair of $2$ roots. Hence, $a(3)=4$.

Another natural question about $p(n)$ is how fast it grows. In the beginnig of
the $20$th century G.H. Hardy and Ramanujan showed the following asymptotics
with their celebrated Circle method [8]

$p(n)\sim\frac{1}{4n\sqrt{3}}\exp\bigg{(}\pi\sqrt{\frac{2n}{3}}\bigg{)},\quad
n\rightarrow\infty.$ (1.3)

Hardy and Ramanujan made extensive use of the fact that the generating
function (1.1) is almost a modular form. More precisely,

$\sum_{n=0}^{\infty}p(n)q^{n}=\prod_{n=1}^{\infty}\frac{1}{1-q^{n}}=\frac{q^{\frac{1}{24}}}{\eta(\tau)}$
(1.4)

,

where $\eta(z)$ is the Dedekind eta function defined on the upper half-plane
$\mathbb{H}=\\{\tau\in\mathbb{C}\mid\text{Im}\left(\tau\right)>0\\}$ by

$\eta(\tau)=q^{\frac{1}{24}}\prod_{n=1}^{\infty}(1-q^{n}),\quad q=\exp(2\pi
i\tau).$

The function $\eta(\tau)$ is a modular form of weight $\frac{1}{2}$ and
satisfies the transformation law, for
$\gamma\in\operatorname{SL}_{2}(\mathbb{Z}),$

$\eta\bigg{(}\frac{a\tau+b}{c\tau+d}\bigg{)}=\varepsilon(a,b,c,d)\sqrt{-i(c\tau+d)}\eta(\tau),\quad\gamma=\begin{pmatrix}a&b\\\
c&d\end{pmatrix},$

$\varepsilon(a,b,c,d):=\begin{cases}e{\frac{\pi ib}{12}}&\text{if }c=0,d=1,\\\
e^{\pi i\big{(}\frac{a+d}{12}-s(d,c)\big{)}}&\text{if }c>0.\end{cases}$ (1.5)

Here $s(h,k)$ is the Dedekind sum

$s(h,k):=\sum_{n=1}^{k-1}\frac{n}{k}\bigg{(}\frac{hn}{k}-\bigg{\lfloor}\frac{hn}{k}\bigg{\rfloor}-\frac{1}{2}\bigg{)}.$

A proof can be found in [3]. A couple years later Rademacher [9] perfected the
circle method and obtained the exact formula

$p(n)=\frac{1}{\pi\sqrt{2}}\sum_{k=1}^{\infty}A_{k}(n)\sqrt{k}\cdot\frac{d}{dn}\left(\frac{\sinh\bigg{[}\frac{\pi}{k}\sqrt{\frac{2}{3}\left(n-\frac{1}{24}\right)}\bigg{]}}{\sqrt{n-\frac{1}{24}}}\right),$
(1.6)

where $A_{k}(n)$ is the Kloosterman sum given by

$A_{k}(n):=\sum_{\begin{subarray}{c}0\leq m<k\\\
\gcd(m,k)=1\end{subarray}}\exp\left(\pi
i\left(s(m,k)-\frac{2nm}{k}\right)\right).$

Note that the series (1.6) converges really fast and that the first term
recovers the result by Hardy and Ramanujan (1.3). Here we want to find an
exact formula of Rademacher type for the cubic partition function $a(n),$ by
using a result of Zuckerman [12], extending Rademacher’s work and it can be
seen as the apotheosis of the classical Circle Method. Using the Circle
Method, Zuckerman computed exact formulas for Fourier coefficients of weakly
holomorphic modular forms of arbitrary non-positive weight on subgroups of
$\operatorname{SL}_{2}(\mathbb{Z})$ with finite index in terms of the cusps of
the underlying subgroup and the negative coefficients of the form at each
cusp. We state the relevant results of Zuckerman in Section $2$ and prove the
exact formula for $a(n)$ in Section $3.$ Similar to the case above we will
obtain as an immediate consequence the asymptotics

$a(n)\sim\frac{e^{\pi\sqrt{n-\frac{1}{8}}}}{8\left(n-\frac{1}{8}\right)^{\frac{5}{4}}},n\rightarrow\infty.$

Finally, we are going to prove a conjecture by Banerjee, Paule, Radu and Zeng
[5] which predicts the following asymptotic formula for $\log a(n)$

$\log\left(a(n)\right)\sim\pi\sqrt{n}-\frac{5}{4}\log\left(n\right)-\log\left(8\right)-\left(\frac{15}{8\pi}+\frac{\pi}{16}\right)\frac{1}{\sqrt{n}},\quad
n\rightarrow\infty,\quad\frac{15}{8\pi}+\frac{\pi}{16}\approx 0.79$

## Acknowledgements

The author wishes to thank Kathrin Bringmann and Walter Bridges for suggesting
this problem, William Craig and Andreas Mono for sharing their knowledge on
modular forms and helpful suggestions, and Johann Franke for helping me verify
the results numerically. The author recieved funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (grant agreement No. 101001179).

##  2\. Zuckerman’s result

The following two sections follow Zuckerman’s work [12] closely. Most of the
results are completely taken over and translated into more modern language.
His method, due to its generality involves numerous technical parameters which
appear in the Cirlce Method and do not have obvious meaning without knowing
the context they arise in. The interested reader should thus consult [12] for
a more detailed account.

Let $\Gamma$ be a subgroup of $\operatorname{SL}_{2}(\mathbb{Z})$ of finite
index and let $F$ be a weakly holomorphic modular form of weight $k=-r,r>0.$
Thus, $F(\tau)$ satisfies a transformation equation of the form

$F\left(\frac{a\tau+b}{c\tau+d}\right)=\varepsilon(-i(c\tau+d))^{k}F(\tau),\quad\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\in\Gamma,$

where $\varepsilon=\varepsilon(a,b,c,d)$ lies on the unit circle and depends
only on the transformation. If $c\neq 0,$ then $c$ is taken to be positive and
we choose the branch of the argument such that

$-\frac{\pi}{2}<\operatorname{Arg}(-i(c\tau+d))<\frac{\pi}{2}.$

Since $\Gamma$ is of finite index in the modular group we can choose a
complete finite system of inequivalent cusps of $\Gamma,$ which we will denote
by $P_{1},\dots P_{s},$ where

$P_{g}=\frac{p_{q}}{q_{g}},\quad\gcd(p_{q},g_{q})=1,\quad g_{q}>0,\quad
g=1,\cdots,s.$

In order to treat all our cusps symetrically, we will have to assume that the
point at infinity does not belong to our chosen set of inequivalent cusps.
This can always be achieved by considering a $\Gamma$-equivlant rational point
instead. We will obtain a set of Fourier expansions $f_{1},\cdots,f_{s}$
corresponding to our set of cusps in the variable $(\tau-P_{g})^{-1}.$
Consider now any transformation $\gamma\in\operatorname{SL}_{2}(\mathbb{Z}),$
which must not necessarily belong to our subgroup $\Gamma.$ We then write

$F\left(\frac{a\tau+b}{c\tau+d}\right)=\varepsilon^{*}(-i(c\tau+d))^{k}F^{*}(\tau),\quad\gamma=\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\in\operatorname{SL}_{2}(\mathbb{Z}),$ (2.1)

where $F^{*}(\tau)$ is now a usual Fourier expansion in $\tau,$ obtained from
one of the expansions

$f_{g}(x)=\sum_{m=-\mu_{g}}^{\infty}a_{m}^{(g)}x^{m}\quad
x=\exp\left(\frac{-2\pi i}{c_{g}(\tau-P_{g})}\right),$

where $c_{g}>0$ will be specified below. Note that by choosing the identity
matrix in (2.1) we get a usual Fourier expansion $F(\tau),$ which we will
later want. Corresponding to each of the cusps $P_{g}$ we can find a
transformation in $\Gamma$

$\tau^{\prime}=\frac{a\tau+b}{c\tau+d},$

which can be written in the the form

$\frac{1}{\tau^{\prime}-P_{g}}=\frac{1}{\tau-P_{g}}+c_{g},\quad c_{g}>0.$

We can recover the transformation in the usual form as

$\tau^{\prime}=\frac{(c_{g}P_{g}+1)-c_{g}P_{g}^{2}}{c_{g}\tau+1-c_{g}P_{g}}.$

We then define $\alpha_{g}$ by the equation $\varepsilon e^{\frac{\pi
i}{2}}=e^{-2\pi i\alpha_{g}},$ where $0\leq\alpha_{g}<1$ and $\varepsilon$ is
taken with respect to $\tau^{\prime}.$

We now describe how to find $F^{*}(\tau)$ in terms of the Fourier expansions
$f_{g}.$ We therefore choose any transformation
$\gamma\in\operatorname{SL}_{2}(\mathbb{Z}).$ Then, $\frac{a}{c}$ is a
rational point if $c\neq 0$ or the point at infinity if $c=0.$ In any case
$\frac{a}{c}$ is $\Gamma$-equivalent to exactly one of our cusps $P_{g}.$ We
can therefore find a transformation

$\frac{a_{1}P_{g}+b_{1}}{c_{1}P_{g}+d_{1}}=\frac{a}{c},\quad c\geq
0,\quad\begin{pmatrix}a_{1}&b_{1}\\\ c_{1}&d_{1}\end{pmatrix}\in\Gamma,$

which takes $P_{g}$ into $\frac{a}{c}.$ We can then express $F^{*}(\tau)$ in
terms of the Fourier expansion $f_{g}$ as follows

$F^{*}(\tau)=q_{g}^{k}\exp\left(\frac{2\pi
iq_{g}^{2}}{c_{g}}\alpha_{g}\tau\right)\dot{f}_{g}\left(\exp\left(\kappa\frac{2\pi
iq_{g}}{c_{g}}(a_{1}d-c_{1}b)\right)\exp\left(\frac{2\pi
iq_{g}^{2}}{c_{g}}\tau\right)\right),$

where $\kappa=\pm 1$ is taken as the solution to the following equations

$p_{q}=\kappa(ad_{1}-cb_{1}),\quad q_{g}=\kappa(ca_{1}-ac_{1}).$

Furthermore, one obtains

$\varepsilon*=\varepsilon(a_{1},b_{1},c_{1},d_{1})\exp\left(\kappa\frac{\pi
ir}{2}\right)\exp\left(\kappa\frac{2\pi
iq_{g}}{c_{g}}\alpha_{g}(a_{1}d-c_{1}b)\right).$

This gives the Fourier expansion in spirit of (2.1) at any rational point
$\frac{a}{c},$ including the point at infinity. When executing the circle
method with these transformation equations there are numerous parameters that
appear. First of all consider as usual the Farey fractions
$\frac{h}{k},\gcd(h,k)=1,k>0,h\geq 0.$ Then, consider the point
$P=P_{g}-\frac{k}{hc_{g}}$ which is either a rational point or the point at
infinity. In any case it is $\Gamma$-equivalent to one of our cusps
$P_{\beta},$ where $\beta=\beta(h,k,g).$ There exists a transformation in
$\Gamma$ which takes this point into $P_{\beta}$

$\frac{aP+b}{cP+d}=P_{\beta},\quad\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\in\Gamma.$

Associated to this transformation is the parameter $\sigma_{h,k}^{(g)}$
implicitly defined as the solution of the equations

$a\left(\frac{c_{g}}{q_{g}}p_{g}h-k\right)+bc_{g}h=\sigma_{h,k}^{(g)}p_{\beta},$

$c\left(\frac{c_{g}}{q_{g}}p_{g}h-k\right)+dc_{g}h=\sigma_{h,k}^{(g)}q_{\beta}.$

It satisfies $-c_{g}\leq\sigma_{h,k}^{(g)}\leq c_{g},\sigma_{h,k}^{(g)}\neq
0.$ The remaining technical parameters necessary to state Zuckerman’s result
are

$\delta_{h,k}^{(g)}=\begin{cases}-1&\text{if }\sigma_{h,k}^{(g)}>0,\\\
1&\text{if }\sigma_{h,k}^{(g)}<0,\end{cases}\quad
G_{h,k}^{(g)}=-\frac{2\pi}{kc_{\beta}}\sigma_{h,k}^{(g)}q_{\beta}(cP_{g}+d),$
$\Omega_{h,k}=\varepsilon^{-1}(a,b,c,d)\exp\left((1-\delta_{h,k})\frac{\pi
ir}{2}\right)\exp\left(-\frac{2\pi
i}{k}\left(\alpha_{g}h+\frac{\sigma_{h,k}^{(g)}q_{\beta}(cP_{g}+d)\alpha_{\beta}}{c\beta}\right)\right).$

Finally, Zuckerman’s result reads now as follows.

######  Theorem 2.1 (Zuckerman’s exact formula).

Assume the notation and hypotheses above and denote by $I_{r}$ the Bessel
function of order $r.$ If $n+\alpha_{g}>0,$ then we have

$a_{n}^{(g)}=\frac{2\pi
c_{g}^{\frac{r-1}{2}}}{(\alpha_{g}+n)^{\frac{r+1}{2}}}\sum_{k=1}^{\infty}\sum_{\begin{subarray}{c}0\leq
h<k\\\ (h,k)=1\end{subarray}}\Omega_{h,k}^{(g)}\exp\bigg{(}-2\pi
in\frac{h}{k}\bigg{)}\\\
\times\sum_{\nu=1}^{\mu_{\beta}}a_{-\nu}^{\beta}\exp\left(-\nu
G_{h,k}^{(g)}i\right)\frac{\left|\sigma_{h,k}^{(g)}\right|q_{\beta}(\nu-\alpha_{\beta})^{\frac{r+1}{2}}}{kc_{\beta}^{\frac{r+1}{2}}}\dot{I}_{r+1}\left(\frac{4\pi|\sigma_{h,k}^{(g)}|q_{\beta}\sqrt{(\nu-\alpha_{\beta})(\alpha_{g}+n)}}{k\sqrt{c_{g}c_{\beta}}}\right).$

##  3\. Exact formula for the cubic partition function

We recall that the cubic partition function was defined by the identity

$\sum_{n=0}^{\infty}a(n)q^{n}=\prod_{n=1}^{\infty}\frac{1}{(1-q^{n})(1-q^{2n})},$

which is up to a power of $q$ equal to $(\eta(\tau)\eta(2\tau))^{-1}.$ Using
the relation

$2\frac{a\tau+b}{c\tau+d}=\frac{a(2\tau)+2b}{\frac{c}{2}(2\tau)+d}$

we can view $\eta(2\tau)$ in light of (1.5) as a modular form of weight
$\frac{1}{2}$ on the subgroup $\Gamma_{0}(2)$ of
$\operatorname{SL}_{2}(\mathbb{Z})$ with multiplier

$\varepsilon=\varepsilon\left(a,2b,\frac{c}{2},d\right),$

Hence, we find that

$F(\tau)=\frac{1}{\eta(\tau)\eta(2\tau)}$

is a weakly modular form of weight $-1$ with multiplier and thus we can apply
Zuckerman’s result to find the Fourier coefficients of $F(\tau)$ and therefore
a formula of $a(n).$ The group $\Gamma_{0}(2)$ has two inequivalent cusps $0$
and $\frac{1}{2}.$ We consider the transformations

$\tau^{\prime}=\frac{\tau}{2\tau+1},\quad\tau^{\prime}=\frac{3\tau-1}{4\tau-1}$

which we write as

$\frac{1}{\tau^{\prime}}=\frac{1}{\tau}+2,\quad\frac{1}{\tau^{\prime}-\frac{1}{2}}=\frac{1}{\tau-\frac{1}{2}}+4.$

Thus, we find that $c_{1}=2$ and $c_{2}=4.$ Moreover we can now determine
$\alpha_{1}$ and $\alpha_{2}.$ Since $F(\tau)$ has the multiplier

$\varepsilon=\varepsilon(a,b,c,d)^{-1}\varepsilon\left(a,2b,\frac{c}{2},d\right)^{-1},$

a short calculation shows that $\alpha_{1}=\alpha_{2}=\frac{7}{8}.$ Since the
point at infinity is $\Gamma_{0}(2)$-equivalent to $\frac{1}{2},$ we obtain
that $F^{*}(\tau))$ is given in terms of the Fourier expansion $f_{2}.$ To
find the exact expression for $F^{*}(\tau))$ we follow Zuckerman’s method by
choosing a transformation that maps $\frac{1}{2}$ to the point at infinity.
Obvioulsy

$\begin{pmatrix}-1&0\\\ 2&-1\end{pmatrix}\frac{1}{2}=i\infty.$

Hence, $\kappa=-1$ and

$F^{*}(\tau)=\frac{1}{2}\exp\left(\frac{7\pi i}{8}\tau\right)f_{1}(-\exp(2\pi
i\tau))=\frac{1}{2}\exp\left(\frac{14\pi
i}{8}\right)\sum_{m=-\mu_{1}}^{\infty}(-1)^{m}a_{m}^{(1)}.$

Furthermore, we consider the point $\frac{a}{c}=\frac{0}{1},$ which
corresponds to the transformation $\tau\mapsto-\frac{1}{\tau}.$ In this case
the obvious equivalent point is $0$ itself and the identity is the
corresponding transformation. Thus, in this case

$F^{*}(\tau)=\exp\left(\frac{7\pi i}{8}\tau\right)f_{2}(\exp(\pi i\tau)).$

We now determine $a_{-1}^{(g)}$ and $\mu_{-2}^{(g)}.$ In the first case we
have on the one hand

$F\left(\begin{pmatrix}1&0\\\
0&1\end{pmatrix}\tau\right)=F(\tau)=\varepsilon^{*}(-1,0,2,-1)\frac{1}{2}\exp\left(\frac{14\pi
i}{8}\tau\right)\sum_{m=-\mu_{1}}^{\infty}(-1)^{m}a_{m}^{(1)}$

with $\varepsilon^{*}=e^{\frac{7\pi i}{8}}.$ On the other hand we have by
definition that

$F(\tau)=e^{-\frac{\pi i\tau}{4}}\left(a(0)+a(1)e^{2\pi
i\tau}+\dots\right)=e^{\frac{2\pi i7}{8}\tau}\left(a(0)e^{-2\pi
i\tau}+a(1)+\dots\right)$ (3.1)

Comparing coefficients we see that $\mu_{-1}^{(2)}=1$ and
$a_{-1}^{(2)}=-2e^{\frac{\pi i}{8}}.$ Furthermore, (3.1) implies that the
Fourier coefficients are exactly $a(n)$ shifted down by one. Thus, we can
recover $a(n)$ through the relation $n=m+1.$ A similar calculation shows that
$\mu_{-1}^{(1)}=1$ and $a_{-1}^{(2)}=\frac{e^{-\frac{\pi i}{2}}}{\sqrt{2}}.$

We now continue to determine the remaining parameters only for $g=2$ as the
Fourier expansion at $\frac{1}{2}$ gives us the exact formula. For
$\beta(h,k,2)$ it is required to find transformations that send
$P_{g}-\frac{k}{c_{g}h}$ into $P_{\beta}.$ Let
$h^{\prime},h^{\prime\prime},h^{\prime\prime\prime}$ be any solutions of

$\left(\frac{k}{2}-h\right)h^{\prime}\equiv 1\,\,({\rm mod}\,\,{2h}),\quad
h^{\prime}>0,\quad\text{ for }k\equiv 0\,\,({\rm mod}\,\,{4}),$
$kh^{\prime\prime}\equiv-1\,\,({\rm mod}\,\,{h}),\quad
h^{\prime\prime}>0,\quad\text{ for }k\equiv 0\,\,({\rm mod}\,\,{2}),$
$(k-2h)h^{\prime\prime\prime}\equiv 1\,\,({\rm mod}\,\,{4h}),\quad
h^{\prime\prime\prime}>0,\quad\text{ for }k\equiv 1\,\,({\rm mod}\,\,{2}).$

Then, the following table gives the remaining parameters for needed to apply
Theorem 2.1

$k$ |  | $a$ |  | $b$ |  | $c$ |  | $d$ | $\beta$ | $\sigma_{h,k}^{(2)}$ | $\delta_{h,k}^{(2)}$ |  | $G_{h,k}^{(2)}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
$\equiv 1\,\,({\rm mod}\,\,{2})$ |  | $h^{\prime\prime\prime}$ |  | $\frac{(k-2h)h^{\prime\prime\prime}-1}{4h}$ |  | $2h^{\prime\prime\prime}+4h$ |  | $\frac{(k-2h)h^{\prime\prime\prime}-1}{2h}-2h+k$ | $2$ | $-1$ | $1$ |  | $\left(\frac{h^{\prime\prime\prime}k-1}{2kh}+1\right)\pi$
$\equiv 0\,\,({\rm mod}\,\,{4})$ |  | $h^{\prime}$ |  | $\frac{\left(\frac{k}{2}-h\right)h^{\prime}-1}{2h}$ |  | $2h^{\prime}+2h$ |  | $\frac{\left(\frac{k}{2}-h\right)h^{\prime}-1}{h}+\frac{k}{2}-h$ | $2$ | $-2$ | $1$ |  | $\left(\frac{h^{\prime}k-2}{kh}+1\right)\pi$
$\equiv 2\,\,({\rm mod}\,\,{4})$ |  | $h$ |  | $\frac{k-2h}{4}$ |  | $4h^{\prime\prime}$ |  | $\frac{kh^{\prime\prime}+1}{h}-2h^{\prime\prime}$ | $1$ | $4$ | $-1$ |  | $-\left(\frac{4h^{\prime\prime}k+4}{kh}\right)\pi$

Since the series in Theorem 2.1 converges absolutely, the terms can be
rearranged and putting everything together, keeping in mind the shift $n=m+1,$
it follows that

$\displaystyle a(n)=$ $\displaystyle\frac{e^{-\frac{\pi
i}{2}}}{2\sqrt{2}}\sum_{\begin{subarray}{c}k\equiv 2\,\,({\rm mod}\,\,{4})\\\
k>0\end{subarray}}A_{k}(n+1)\frac{1}{n-\frac{1}{8}}I_{2}\left(\frac{2\pi\sqrt{n-\frac{1}{8}}}{k}\right)$
$\displaystyle-\frac{\pi e^{\frac{\pi
i}{8}}}{2}\sum_{\begin{subarray}{c}k\equiv 0\,\,({\rm mod}\,\,{4})\\\
k>0\end{subarray}}\frac{1}{k}A_{k}(n+1)\frac{1}{n-\frac{1}{8}}I_{2}\left(\frac{\sqrt{2}\pi\sqrt{n-\frac{1}{8}}}{k}\right)$
$\displaystyle-\frac{\pi e^{\frac{\pi
i}{8}}}{4}\sum_{\begin{subarray}{c}k\equiv 1\,\,({\rm mod}\,\,{2})\\\
k>0\end{subarray}}\frac{1}{k}A_{k}(n+1)\frac{1}{n-\frac{1}{8}}I_{2}\left(\frac{\sqrt{2}\pi\sqrt{n-\frac{1}{8}}}{2k}\right),$
(3.2)

where

$A_{k}(n)=\sum_{\begin{subarray}{c}0\leq h<k\\\
(h,k)=1\end{subarray}}\varepsilon\left(h,\frac{k-2h}{4},4h^{\prime\prime},\frac{kh^{\prime\prime}+1}{h}-2h^{\prime\prime}\right)^{-1}\exp(\pi
i)\\\ \times\exp\left(-\frac{2\pi
i7}{8k}h+\frac{1}{8}\left(\frac{2h^{\prime\prime}k+2}{kh}\right)\pi i-2\pi
in\frac{h}{k}\right),\quad k\equiv 2\,\,({\rm mod}\,\,{4}),$ (3.3)

$A_{k}(n)=\sum_{\begin{subarray}{c}0\leq h<k\\\
(h,k)=1\end{subarray}}\varepsilon\left(h^{\prime},\frac{\left(\frac{k}{2}-h\right)h^{\prime}-1}{2h},2h^{\prime}+2h^{\prime},\frac{\left(\frac{k}{2}-h\right)h^{\prime}-1}{h}+\frac{k}{2}-h\right)^{-1}\exp(\pi
i)\\\ \times\exp\left(-\frac{2\pi
i7}{8k}h-\frac{1}{8}\left(\frac{h^{\prime}k-2}{kh}+1\right)\pi i-2\pi
in\frac{h}{k}\right),\quad k\equiv 0\,\,({\rm mod}\,\,{4}),$ (3.4)

$A_{k}(n)=\sum_{\begin{subarray}{c}0\leq h<k\\\
(h,k)=1\end{subarray}}\varepsilon\left(h^{\prime\prime\prime},\frac{(k-2h)h^{\prime\prime\prime}-1}{4h},2h^{\prime\prime\prime}+4h,\frac{(k-2h)h^{\prime\prime\prime}-1}{2h}-2h+k\right)^{-1}\exp(\pi
i)\\\ \times\exp\left(-\frac{2\pi
i7}{8k}h-\frac{1}{8}\left(\frac{h^{\prime\prime\prime}k-1}{2kh}+1\right)\pi
i-2\pi in\frac{h}{k}\right),\quad k\equiv 2\,\,({\rm mod}\,\,{4}),$ (3.5)

Changing the order of summation one last time we arrive at the following

######  Theorem 3.1.

For $n\geq 0$ we have the exact formula

$\displaystyle a(n)=\frac{\pi e^{-\frac{\pi
i}{2}}}{4\sqrt{2}\left(n-\frac{1}{8}\right)}\sum_{\begin{subarray}{c}l\equiv
1\,\,({\rm mod}\,\,{2})\\\
l>0\end{subarray}}\frac{1}{l}A_{l^{\prime}}(n+1)I_{2}\left(\frac{\pi\sqrt{n-\frac{1}{8}}}{l}\right)$
$\displaystyle-\frac{\pi e^{\frac{\pi
i}{8}}}{2\left(n-\frac{1}{8}\right)}\sum_{\begin{subarray}{c}l\equiv
0\,\,({\rm mod}\,\,{2})\\\
l>0\end{subarray}}\frac{1}{l}A_{l^{\prime}}(n+1)I_{2}\left(\frac{\sqrt{2}\pi\sqrt{n-\frac{1}{8}}}{l}\right),$
(3.6)

where $l^{\prime}=l$ if $l\equiv 1\,\,({\rm mod}\,\,{2}),$ $l^{\prime}=l$ if
$l\equiv 0\,\,({\rm mod}\,\,{4}),$ $l^{\prime}=\frac{l}{2}$ if $l\equiv
2\,\,({\rm mod}\,\,{4}).$

It is easy to see that the series converges rapidly and that the term $l=1$
dominates the rest of the sum. Thus, we may conclude that

$a(n)\sim\frac{\pi}{4\sqrt{2}\left(n-\frac{1}{8}\right)}I_{2}\left(\pi\sqrt{n-\frac{1}{8}}\right).$

For the Bessel function of order $\nu\in\mathbb{N}$ the following asymptotics
are known,

$I_{\nu}(z)\sim\frac{e^{z}}{\sqrt{2\pi
z}}\sum_{k=0}^{\infty}(-1)^{k}\frac{a_{k}(\nu)}{z^{k}},$ (3.7)

see [11], where $a_{k}(\nu)$ is defined as

$a_{k}(\nu)=(-1)^{k}\frac{\left(\frac{1}{2}-\nu\right)\left(\frac{1}{2}+\nu\right)}{2^{k}k!},\quad
k=0,1,2,\dots$ (3.8) $(\lambda)_{n}=\frac{\Gamma(\lambda+n)}{\Gamma(\lambda)}$
(3.9)

$(\lambda)_{0}=1$ (3.10)

Thus, we find that the first order approximation is

$I_{2}(z)=\frac{e^{z}}{\sqrt{2\pi z}}\left(1+O\left(\frac{1}{z}\right)\right)$

We now obtain the following asymptotics.

######  Corollary 3.2.

We have

$a(n)\sim\frac{e^{\pi\sqrt{n-\frac{1}{8}}}}{8\left(n-\frac{1}{8}\right)^{\frac{5}{4}}}$
(3.11)

As a final consequence we prove a conjecture by Banerjee, Paule, Radu and
Zeng.

######  Corollary 3.3.

We have

$\log a(n)\sim\pi\sqrt{n}-\frac{5}{4}\log n-\log
8-\left(\frac{15}{8\pi}+\frac{\pi}{16}\right)\frac{1}{\sqrt{n}},\quad
n\rightarrow\infty.$

###### Proof.

We need to shift our approximation from $n-\frac{1}{8}$ to $n.$ Thus, write

$e^{\pi\sqrt{n-\frac{1}{8}}}=e^{\pi\sqrt{n}}\left(1-\frac{\pi}{16\sqrt{n}}+O\left(\frac{1}{n}\right)\right),$

$\frac{1}{8\left(n-\frac{1}{8}\right)^{\frac{5}{4}}}=\frac{1}{8n^{\frac{5}{4}}}\left(1+O\left(\frac{1}{n}\right)\right).$

Together with the second order approximation for the Bessel function

$I_{2}(z)=\frac{e^{z}}{\sqrt{2\pi
z}}\left(1-\frac{15}{8z}+O\left(\frac{1}{z^{2}}\right)\right)$

we obtain

$a(n)=\frac{e^{\pi\sqrt{n}}}{8n^{\frac{5}{4}}}\left(1-\frac{15}{8n}+O\left(\frac{1}{n^{2}}\right)\right)\left(1-\frac{\pi}{16\sqrt{n}}+O\left(\frac{1}{n}\right)\right)\left(1-\frac{\pi}{16\sqrt{n}}+O\left(\frac{1}{n}\right)\right)\left(1+O\left(\frac{1}{n}\right)\right)$

Taking logarithms, we find as conjectured

$\log a(n)\sim\pi\sqrt{n}-\frac{5}{4}\log n-\log
8-\left(\frac{15}{8\pi}+\frac{\pi}{16}\right)\frac{1}{\sqrt{n}},\quad
n\rightarrow\infty.$

∎

######  Remark.

We can refine our proof and obtain stronger asymptotics of the form

$\log a(n)\sim\pi\sqrt{n}-\frac{5}{4}\log n-\log
8-\left(\frac{15}{8\pi}+\frac{\pi}{16}\right)\frac{1}{\sqrt{n}}+\sum_{j=2}^{\infty}\frac{c_{j}}{n^{\frac{j}{2}}}$

for some real $c_{i}$ by taking higher order expansions of the involved
functions.

## References

  * [1] G. Andrews, _The Theory of Partitions_ , Encycl. Math. and Its Appl., vol. 2, G.-C. Rota ed., Addison-Wesley, Reading, 1976 (Reissued: Cambridge University Press, Cambridge 1998).
  * [2] G. Andrews and B. Berndt, _Ramanujan’s Lost Notebook_ , Part I, SpringerVerlag, New York, 2005.
  * [3] T. Apostol, _Modular Functions and Dirichlet Series in Number Theory_ , Springer-Verlag, 1976
  * [4] B. Berndt, _Ramanujan’s Notebooks, Part III_ , Springer-Verlag, New York, 1991\.
  * [5] K. Banerjee, P. Paule, C. Radu and W. Zeng, _New inequalities for $p(n)$ and $\log p(n)$_, Ramanujan J. (2022).
  * [6] H. Chan, _Ramanujan’s cubic continued fraction and an analog of his ”most beautiful identity”_ , International Journal of Number Theory Vol. 06, No. 03, pp. 673-680 (2010).
  * [7] W. Chen and B. Lin, _Congruences for the Number of Cubic Partitions Derived from Modular Forms_ , arxiv:0910.1263
  * [8] G. Hardy and S. Ramanujan, _Asymptotic formulae in combinatory analysis_ , Proc. London Math. Soc. Ser. 2 17 (1918), 75–115.
  * [9] H. Rademacher, _A convergent series for the partition function $p(n)$_, PNAS February 1, 1937 23 (2), 78–84.
  * [10] S. Ramanujan, _Collected Papers_ , Cambridge University Press, Cambridge, 1927; reprinted by Chelsea, New York, 1962; reprinted by American Mathematical Society, Providence, 2000.
  * [11] N. Temme, _Asymptotic Methods for Integrals_ , World Scientific (2014).
  * [12] H. Zuckerman, On the coefficients of certain modular forms belonging to subgroups of the modular group, Trans. Amer. Math. Soc. 45 (1939), 298a321.


# Robust Beamforming for IRS Aided MIMO Full Duplex Systems

Chandan Kumar Sheemar, Jorge Querol, Sourabh Solanki, Sumit Kumar, and Symeon
Chatzinotas
SnT, University of Luxembourg, emails<EMAIL_ADDRESS>
###### Abstract

In this paper, a novel robust beamforming for an intelligent reflecting
surface (IRS) assisted FD system is presented. Since perfect channel state
information (CSI) is often challenging to acquire in practice, we consider the
case of imperfect CSI and adopt a statistically robust beamforming approach to
maximize the ergodic weighted sum rate (WSR). We also analyze the achievable
WSR of an IRS-assisted FD with imperfect CSI, for which the lower and the
upper bounds are derived. The ergodic WSR maximization problem is tackled
based on the expected Weighted Minimum Mean Squared Error (WMMSE), which is
guaranteed to converge to a local optimum. The effectiveness of the proposed
design is investigated with extensive simulation results. It is shown that our
robust design achieves significant performance gain compared to the naive
beamforming approaches and considerably outperforms the robust Half-Duplex
(HD) system.

###### Index Terms:

full duplex, intelligent reflecting surfaces, robust beamforming, ergodic
weighted sum rate.

##  I Introduction

Full Duplex (FD) is a promising wireless transmission technology offering
simultaneous transmission and reception in the same frequency band, which
theoretically doubles the spectral efficiency [1, 2]. Beyond spectral
efficiency, FD can be beneficial to improve security, enable advanced joint
communication and sensing, and reduce end-to-end delays. However, FD systems
suffer from Self-Interference (SI), which could be $90-110~{}$dB higher than
the received signal of interest, and it is a major challenge for achieving an
ideal FD operation. Advanced SI Cancellation (SIC) techniques are pivotal to
mitigate the SI power close to the noise floor and enable correct reception of
the received signal. In parallel to FD, another emerging technology is the
Intelligent Reflecting Surfaces (IRSs), which can provide a smart and
reconfigurable wireless environment [3, 4, 5].

The IRS-aided FD (IRS-FD) systems hold great promise in supporting the
forthcoming traffic demands while being highly spectrally and energy efficient
[6, 7, 8]. In [3], the authors investigated the performance of the IRS-FD
system by analyzing the outage and error probabilities. In [9], the authors
present a novel beamforming design for an FD relay assisted with one IRS to
maximize the minimum achievable rate with the max-min optimization. In [10],
the authors studied the IRS-FD systems to improve the worst-case achievable
security rate under imperfect channel state information (CSI). Furthermore, a
deep neural network is also developed to reduce the computational burden of
the proposed beamforming solution. Finally, in [11], the authors presented a
joint beamforming design for weighted sum rate (WSR) maximization in a single-
cell multiple-input single-output (MISO) IRS-FD system. Note that the existing
literature only provides the secrecy rate expression for IRS-FD systems under
imperfect Channel State Information (CSI). However, there is currently no
analytical expression available for the ergodic weighted sum rate (WSR) or a
novel robust beamforming design considering imperfect CSI in this context.

Motivated by the aforementioned considerations, we address the challenge of
maximizing the ergodic WSR for IRS-FD systems with multi-antenna users under
imperfect CSI. Our objective is to enhance both the uplink (UL) and downlink
(DL) sum rates by jointly optimizing the digital beamformers and the IRS phase
response while accounting for imperfect CSI. It is important to note that
incorporating CSI errors in the beamforming designs, with multi-antenna users,
makes the problem of MIMO beamforming for IRS-FD systems highly complex. The
CSI errors are modelled by using the Gaussian-Kronecker model and then an
analytical approximation for the ergodic achievable WSR using this model is
derived. Subsequently, we formulate the problem of maximizing the ergodic WSR,
considering average sum-power and unit-modulus constraints for the IRS phase
response under imperfect CSI. To solve this problem, we leverage the
connection between the ergodic WSR and the expected weighted minimum mean
squared error (EWMMSE), where the expectation is taken with respect to the CSI
errors. Simulation results validate the accuracy of the analytically derived
ergodic WSR expression. It is further shown that our proposed approach
significantly outperforms the naive schemes in the presence of imperfect
CSI.111 _Notations:_ Boldface lower and upper case characters denote vectors
and matrices, respectively. $\mathbb{E}\\{\cdot\\}$, $\mbox{Tr}\\{\cdot\\}$,
$\mathbf{I}$, $\otimes$, $\widetilde{\mathbf{X}}$ denote expectation, trace,
identity matrix, Kronecker product, and transmit covariance matrix,
respectively. Finally, $\mbox{diag}(\mathbf{x})$ and denote a diagonal matrix
with vector $\mathbf{x}$ on its main diagonal and $\mathbf{I}$ denotes
identity matrix.

##  II System Model

Let $j$ and $k$ denote the multi-antenna DL and UL user served by the MIMO FD
base station (BS), respectively, and let $M_{k}$ and $N_{j}$ denote their
number of transmit and receive antennas, respectively. The FD BS is assumed to
have $M_{0}$ transmit and $N_{0}$ receive antennas. We consider a multi-stream
approach, and the number of data streams for the UL user $k$ and DL user $j$
is denoted as $u_{k}$ and $v_{j}$, respectively. Let
$\mathbf{U}_{k}\in\mathbb{C}^{M_{k}\times u_{k}}$ and
$\mathbf{V}_{j}\in\mathbb{C}^{M_{t}\times v_{j}}$ denote the precoders for
white unitary-variance data streams $\mathbf{s}_{k}\in\mathbb{C}^{u_{k}\times
1}$ and $\mathbf{s}_{j}\in\mathbb{C}^{v_{j}\times 1}$, respectively. We assume
that the considered FD system is aided with one IRS of size $R\times C$. Let
$\bm{\theta}=[e^{i\theta_{1}},....,e^{i\theta_{RC}}]$ denote the vector
containing the phase-shift response of its $RC$ elements, and let
$\mathbf{\Theta}=\mbox{diag}(\bm{\theta})$ denote a diagonal matrix containing
$\bm{\theta}$ on its main diagonal. Let $\mathbf{n}_{0}$ and $\mathbf{n}_{j}$
denote the noise vectors at the FD BS and DL user $j$, respectively, which are
modelled as $\mathbf{n}_{0}=\mathcal{CN}(\bm{0},\sigma_{0}^{2}\mathbf{I})$
$\mathbf{n}_{j}=\mathcal{CN}(\bm{0},\sigma_{j}^{2}\mathbf{I})$, where
$\sigma_{0}^{2}$ and $\sigma_{j}^{2}$ denote the noise variances at the FD
node and the DL user $j$, respectively. The channel responses from the UL user
$k$ to the BS and from the BS to the DL user $j$ are denoted with
$\mathbf{H}_{k}\in\mathbb{C}^{N_{0}\times M_{k}}$ and
$\mathbf{H}_{j}\in\mathbb{C}^{N_{j}\times M_{0}}$, respectively. Let
$\mathbf{H}_{0}\in\mathbb{C}^{N_{0}\times M_{0}}$ and
$\mathbf{H}_{j,k}\in\mathbb{C}^{N_{j}\times M_{k}}$ denote the SI channel
response for the FD BS and cross-interference channel response between the UL
user $k$ and the DL user $j$, respectively. The channel responses from the
transmit antenna array of the FD BS to the IRS and from the IRS to the receive
antenna array of the FD BS are denoted with
$\mathbf{H}_{\theta,0}\in\mathbb{C}^{RC\times M_{0}}$ and
$\mathbf{H}_{0,\theta}\in\mathbb{C}^{N_{0}\times RC}$, respectively. Let
$\mathbf{H}_{j,\theta}\in\mathbb{C}^{N_{j}\times RC}$ and
$\mathbf{H}_{\theta,k}\in\mathbb{C}^{RC\times M_{k}}$ denote the channel
responses from the IRS to the DL user $j$ and from the UL user $k$ to the IRS,
respectively.

###  II-A Imperfect CSI Modelling

Let
$\Delta\mathbf{H}_{k},\Delta\mathbf{H}_{j},\Delta\mathbf{H}_{0},\Delta\mathbf{H}_{j,k},\Delta\mathbf{H}_{\theta,0},\Delta\mathbf{H}_{0,\theta},\Delta\mathbf{H}_{j,\theta},\Delta\mathbf{H}_{\theta,k}$
denote the estimation errors for the channel responses
$\mathbf{H}_{k},\mathbf{H}_{j},\mathbf{H}_{0},\mathbf{H}_{j,k},\mathbf{H}_{\theta,0},\mathbf{H}_{0,\theta},\mathbf{H}_{j,\theta},\mathbf{H}_{\theta,k}$,
respectively. The available CSI matrices can be written as a sum of the
channel estimates and CSI errors as

$\displaystyle\mathbf{H}_{k}=\hat{\mathbf{H}}_{k}+\Delta\mathbf{H}_{k},\quad\quad\quad\quad\quad\quad\quad\;\mathbf{H}_{j}=\hat{\mathbf{H}}_{j}+\Delta\mathbf{H}_{j},$
(1)
$\displaystyle\mathbf{H}_{0}=\hat{\mathbf{H}}_{0}+\Delta\mathbf{H}_{0},\quad\quad\quad\quad\quad\quad\quad\;\;\mathbf{H}_{j,k}=\hat{\mathbf{H}}_{j,k}+\Delta\mathbf{H}_{j,k},$
$\displaystyle\mathbf{H}_{\theta,0}=\hat{\mathbf{H}}_{\theta,0}+\Delta\mathbf{H}_{\theta,0},\quad\quad\quad\quad\quad\;\mathbf{H}_{0,\theta}=\hat{\mathbf{H}}_{0,\theta}+\Delta\mathbf{H}_{0,\theta}$
$\displaystyle\mathbf{H}_{j,\theta}=\hat{\mathbf{H}}_{j,\theta}+\Delta\mathbf{H}_{j,\theta},\quad\quad\quad\quad\quad\;\mathbf{H}_{\theta,k}=\hat{\mathbf{H}}_{\theta,k}+\Delta\mathbf{H}_{\theta,k},$

where the channel matrices of the form $\hat{\mathbf{X}}$ denote the channel
estimates. To model the estimation errors, we adopt the Gaussian Kronecker
model [12], which dictates that

$\displaystyle\Delta\mathbf{H}_{k}=\mathcal{CN}(\bm{0},\mathbf{J}_{k}\otimes\mathbf{K}_{k}),\quad\;\;\Delta\mathbf{H}_{\theta,0}=\mathcal{CN}(\bm{0},\mathbf{J}_{\theta,0}\otimes\mathbf{K}_{\theta,0}),$
(2)
$\displaystyle\Delta\mathbf{H}_{0}=\mathcal{CN}(\bm{0},\mathbf{J}_{0}\otimes\mathbf{K}_{0}),\quad\quad\Delta\mathbf{H}_{j,k}=\mathcal{CN}(\bm{0},\mathbf{J}_{j,k}\otimes\mathbf{K}_{j,k}),$
$\displaystyle\Delta\mathbf{H}_{j}=\mathcal{CN}(\bm{0},\mathbf{J}_{j}\otimes\mathbf{K}_{j}),\quad\;\quad\quad\Delta\mathbf{H}_{0,\theta}=\mathcal{CN}(\bm{0},\mathbf{J}_{0}\otimes\mathbf{K}_{\theta}),$
$\displaystyle\Delta\mathbf{H}_{j,\theta}=\mathcal{CN}(\bm{0},\mathbf{J}_{j,\theta}\otimes\mathbf{K}_{j,\theta}),\;\;\Delta\mathbf{H}_{\theta,k}\mathcal{CN}(\bm{0},\mathbf{J}_{\theta,k}\otimes\mathbf{K}_{\theta,k}),$

where the matrices $\mathbf{J}$ and $\mathbf{K}$ are the covariance matrices
seen from the transmitter and receiver for each link, respectively. The
estimation errors are uncorrelated with the estimated channel matrices, and
hence we can write

$\displaystyle\mathbf{H}_{k}=\mathcal{CN}(\hat{\mathbf{H}}_{k},\mathbf{J}_{k}\otimes\mathbf{K}_{k}),\;\;\;\mathbf{H}_{\theta,0}=\mathcal{CN}(\hat{\mathbf{H}}_{\theta,0},\mathbf{J}_{\theta,0}\otimes\mathbf{K}_{\theta,0}),$
(3)
$\displaystyle\mathbf{H}_{0}=\mathcal{CN}(\hat{\mathbf{H}}_{0},\mathbf{J}_{0}\otimes\mathbf{K}_{0}),\;\;\;\mathbf{H}_{j,k}=\mathcal{CN}(\hat{\mathbf{H}}_{j,k},\mathbf{J}_{j,k}\otimes\mathbf{K}_{j,k}),$
$\displaystyle\mathbf{H}_{j}=\mathcal{CN}(\hat{\mathbf{H}}_{j},\mathbf{J}_{j}\otimes\mathbf{K}_{j}),\;\;\;\;\mathbf{H}_{\theta,k}=\mathcal{CN}(\hat{\mathbf{H}}_{\theta,k},\mathbf{J}_{\theta,k}\otimes\mathbf{K}_{\theta,k})$
$\displaystyle\mathbf{H}_{j,\theta}=\mathcal{CN}(\hat{\mathbf{H}}_{j,\theta},\mathbf{J}_{j,\theta}\otimes\mathbf{K}_{j,\theta})\mathbf{H}_{0,\theta}=\mathcal{CN}(\hat{\mathbf{H}}_{0,\theta},\mathbf{J}_{0,\theta}\otimes\mathbf{K}_{0,\theta})$

Let
$\overline{\mathbf{H}}_{0,k},\overline{\mathbf{H}}_{0,0},\overline{\mathbf{H}}_{j},\overline{\mathbf{H}}_{j,k}$
denote the effective channel responses affected by the estimation errors
defined as

$\overline{\mathbf{H}}_{k}=(\hat{\mathbf{H}}_{k}+\Delta\mathbf{H}_{k})+(\hat{\mathbf{H}}_{0,\theta}+\Delta\mathbf{H}_{0,\theta})\mathbf{\Theta}(\hat{\mathbf{H}}_{\theta,k}+\Delta\mathbf{H}_{\theta,k})$
(4a)
$\overline{\mathbf{H}}_{0}=(\hat{\mathbf{H}}_{0}+\Delta\mathbf{H}_{0})+(\hat{\mathbf{H}}_{0,\theta}+\Delta\mathbf{H}_{0,\theta})\mathbf{\Theta}(\hat{\mathbf{H}}_{\theta,0}+\Delta\mathbf{H}_{\theta,0})$
(4b)
$\overline{\mathbf{H}}_{j}=(\hat{\mathbf{H}}_{j}+\Delta\mathbf{H}_{j})+(\hat{\mathbf{H}}_{j,\theta}+\Delta\mathbf{H}_{j,\theta})\mathbf{\Theta}(\hat{\mathbf{H}}_{\theta,0}+\Delta\mathbf{H}_{\theta,0})$
(4c)
$\overline{\mathbf{H}}_{j,k}=(\hat{\mathbf{H}}_{j,k}+\Delta\mathbf{H}_{j,k})+(\hat{\mathbf{H}}_{j,\theta}+\Delta\mathbf{H}_{j,\theta})\mathbf{\Theta}(\hat{\mathbf{H}}_{\theta,k}+\Delta\mathbf{H}_{\theta,k})$
(4d)

Let $\mathbf{y}_{k}$ and $\mathbf{y}_{j}$ denote the signals received by the
FD BS from UL user $k$ and by the DL user $j$, respectively, given as

$\mathbf{y}_{k}=\overline{\mathbf{H}}_{k}\mathbf{U}_{k}\mathbf{s}_{k}+\overline{\mathbf{H}}_{0}\mathbf{V}_{j}\mathbf{s}_{j}+\mathbf{n}_{0},$
(5a)
$\mathbf{y}_{j}=\overline{\mathbf{H}}_{j}\mathbf{V}_{j}\mathbf{s}_{j}+\overline{\mathbf{H}}_{j,k}\mathbf{U}_{k}\mathbf{s}_{k}+\mathbf{n}_{j}.$
(5b)

###  II-B Problem Formulation

We aim to maximize the ergodic WSR of the IRS-FD system under imperfect CSI
given the CSI error statistics and the channel estimates. Let
$\mathcal{R}=\mathcal{R}_{k}+\mathcal{R}_{j}$ denote the WSR of the system in
the case of perfect CSI, with $\mathcal{R}_{k}$ and $\mathcal{R}_{j}$ denoting
the weighted rate of the users $k$ and $j$, respectively. The ergodic WSR of
the system, i.e., the average WSR with respect to the CSI errors, can be
written as $\mathbb{E}_{\mathbf{H}|\hat{\mathbf{H}}}[\mathcal{R}]$. However,
such a problem is infeasible to be solved and can be tackled by applying
Jensen’s inequality, which allows moving the expectation operator insider as
$\mathbb{E}_{\mathbf{H}|\hat{\mathbf{H}}}[\mathcal{R}]\geq\mathcal{R}(\mathbb{E}_{\mathbf{H}|\hat{\mathbf{H}}})$
[13]. The optimization problem for ergodic WSR
$\mathcal{R}(\mathbb{E}_{\mathbf{H}|\hat{\mathbf{H}}})$, under the sum-power
and the unit-modulus constraint can be written as

$\underset{\begin{subarray}{c}\mathbf{V}_{j},\mathbf{U}_{k}\end{subarray},\mathbf{\Theta}}{\max}\quad\mathcal{R}_{k}(\mathbb{E}_{\mathbf{H}|\hat{\mathbf{H}}})+\mathcal{R}_{j}(\mathbb{E}_{\mathbf{H}|\hat{\mathbf{H}}})\vspace{-4mm}$
(6a)
$\vspace{-1mm}\text{s.t.}\quad\mbox{Tr}(\mathbf{U}_{k}\mathbf{U}_{k}^{H})\preceq\alpha_{k},\;\&\;\mbox{Tr}(\mathbf{V}_{j}\mathbf{V}_{j}^{H})\leq\alpha_{0},$
(6b) $\vspace{-1mm}\quad\quad\quad|\bm{\theta}(i)|=1,\quad\forall i,$ (6c)

where $\alpha_{0}$ and $\alpha_{k}$ denote the sum-power constraint at the BS
and UL user, respectively.

###  II-C Ergodic WSR Analysis with Imperfect CSI

In the following, we derive the expression for the ergodic WSR
$\mathcal{R}(\mathbb{E}_{\mathbf{H}|\hat{\mathbf{H}}})$ in the presence of
estimated channel responses and CSI errors for the IRS-FD system with multi-
antenna UL and DL users. Let
$\widetilde{\mathbf{U}_{k}}=\mathbf{U}_{k}\mathbf{U}_{k}^{H}$ and
$\widetilde{\mathbf{V}_{j}}=\mathbf{V}_{j}\mathbf{V}_{j}^{H}$ denote the
transmit covariance matrices for UL user $k$ and DL user $j$, respectively. In
the case of imperfect CSI, the received signal plus interference plus noise
covariance matrices $\mathbf{R}_{k}$ and $\mathbf{R}_{j}$, including the CSI
errors and the IRS phase response $\mathbf{\Theta}$, are given as

$\displaystyle\mathbf{R}_{k}=$
$\displaystyle\overline{\mathbf{H}_{k}}\widetilde{\mathbf{U}_{k}}\overline{\mathbf{H}_{k}}^{H}+\overline{\mathbf{H}_{0}}\widetilde{\mathbf{V}_{j}}\overline{\mathbf{H}_{0}}^{H}+\sigma_{0}^{2}\mathbf{I},$
(7a) $\displaystyle\mathbf{R}_{j}=$
$\displaystyle\overline{\mathbf{H}_{j}}\widetilde{\mathbf{V}_{j}}\overline{\mathbf{H}_{j}}^{H}+\overline{\mathbf{H}}_{j,k}\widetilde{\mathbf{U}_{k}}\overline{\mathbf{H}}_{j,k}^{H}+\sigma_{j}^{2}\mathbf{I}.$
(7b)

The interference plus noise covariance matrices can be obtained as
$\mathbf{R}_{\overline{k}}=\mathbf{R}_{k}-\mathbf{S}_{k},\mathbf{R}_{\overline{j}}=\mathbf{R}_{j}-\mathbf{S}_{j}$,
with $\mathbf{S}_{k}$ and $\mathbf{S}_{j}$ denoting the useful received signal
covariance part.

Given the statistical distribution of the CSI errors (2) and the channel
estimates, the lower bound for the ergodic WSR
$\mathcal{R}(\mathbb{E}_{\mathbf{H}|\hat{\mathbf{H}}})$ of an IRS-FD system
with multi-antenna users with imperfect CSI can be approximated as

$\displaystyle\mathcal{R}(\mathbb{E}_{\mathbf{H}|\hat{\mathbf{H}}})$
$\displaystyle=w_{k}\mbox{ln}[\mbox{det}(\mathbf{I}+\mathbf{U}_{k}^{H}(\hat{\mathbf{H}}_{k}+\hat{\mathbf{H}}_{0,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,k})^{H}\mathbf{\Sigma}_{\overline{k}}^{-1}$
(8)
$\displaystyle(\hat{\mathbf{H}}_{k}+\hat{\mathbf{H}}_{0,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,k})\mathbf{U}_{k})]+w_{j}\mbox{ln}[\mbox{det}(\mathbf{I}+\mathbf{V}_{j}^{H}(\hat{\mathbf{H}}_{j}+$
$\displaystyle\hat{\mathbf{H}}_{j,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,0})^{H}\mathbf{\Sigma}_{\overline{j}}^{-1}(\hat{\mathbf{H}}_{j}+\hat{\mathbf{H}}_{j,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,0})\mathbf{V}_{j})],$

$\displaystyle\mathbf{\Sigma}_{\overline{k}}=$
$\displaystyle\mbox{Tr}(\widetilde{\mathbf{U}_{k}}\mathbf{J}_{k}^{T})\mathbf{K}_{k}+\hat{\mathbf{H}}_{0,\theta}\mathbf{\Theta}\mbox{Tr}(\widetilde{\mathbf{U}_{k}}\mathbf{J}_{\theta,k}^{T})\mathbf{K}_{\theta,k}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{0,\theta}^{H}+\mbox{Tr}(\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,k}\widetilde{\mathbf{U}_{k}}\hat{\mathbf{H}}_{\theta,k}\mathbf{\Theta}^{H}\mathbf{J}_{0,\theta}^{T})\mathbf{K}_{0,\theta}+\mbox{Tr}(\widetilde{\mathbf{U}_{k}}\mathbf{J}_{\theta,k}^{T})\mbox{Tr}(\mathbf{\Theta}\mathbf{K}_{\theta,k}\mathbf{\Theta}^{H}\mathbf{J}_{0,\theta}^{T})\mathbf{K}_{0,\theta}$
(9a)
$\displaystyle+\hat{\mathbf{H}}_{0}\widetilde{\mathbf{V}_{j}}\hat{\mathbf{H}}_{0}^{H}+\hat{\mathbf{H}}_{0}\widetilde{\mathbf{V}_{j}}\mathbf{H}_{\theta,0}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{0,\theta}+\mbox{Tr}(\widetilde{\mathbf{V}_{j}}\mathbf{J}_{0}^{T})\mathbf{K}_{0}^{T}+\hat{\mathbf{H}}_{0,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,0}\widetilde{\mathbf{V}_{j}}\hat{\mathbf{H}}_{0}^{H}+\hat{\mathbf{H}}_{0,\theta}\mathbf{\Theta}\mathbf{H}_{\theta,0}\widetilde{\mathbf{V}_{j}}\mathbf{H}_{\theta,0}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{0,\theta}^{H}$
$\displaystyle+\hat{\mathbf{H}}_{0,\theta}\mathbf{\Theta}\mbox{Tr}(\widetilde{\mathbf{V}_{j}}\mathbf{J}_{\theta,0}^{T})\mathbf{K}_{\theta,0}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{0,\theta}^{H}+\mbox{Tr}(\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,0}\widetilde{\mathbf{V}_{j}}\hat{\mathbf{H}}_{\theta,0}^{H}\mathbf{\Theta}^{H}\mathbf{J}_{0,\theta}^{T})\mathbf{K}_{0,\theta}+\mbox{Tr}(\widetilde{\mathbf{V}_{j}}\mathbf{J}_{\theta,0}^{T})\mbox{Tr}(\mathbf{\Theta}\mathbf{K}_{\theta,0}\mathbf{\Theta}^{H}\mathbf{J}_{0,\theta}^{T})\mathbf{K}_{0,\theta}+\sigma_{0}^{2}\mathbf{I},$
$\displaystyle\mathbf{\Sigma}_{\overline{j}}=$
$\displaystyle\mbox{Tr}(\widetilde{\mathbf{V}_{j}}\mathbf{J}_{j}^{T})\mathbf{K}_{j}+\hat{\mathbf{H}}_{j,\theta}\mathbf{\Theta}\mbox{Tr}(\widetilde{\mathbf{V}_{j}}\mathbf{J}_{\theta,0}^{T})\mathbf{K}_{\theta,0}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{j,\theta}^{H}+\mbox{Tr}(\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,0}\widetilde{\mathbf{V}_{j}}\hat{\mathbf{H}}_{\theta,0}^{H}\mathbf{\Theta}^{H}\mathbf{J}_{j,\theta}^{T})\mathbf{K}_{j,\theta}+\mbox{Tr}(\widetilde{\mathbf{V}_{j}}\mathbf{J}_{\theta,0}^{T})\mbox{Tr}(\mathbf{\Theta}\mathbf{K}_{\theta,0}\mathbf{\Theta}^{H}\mathbf{J}_{j,\theta}^{T})\mathbf{K}_{j,\theta}$
(9b)
$\displaystyle+\hat{\mathbf{H}}_{j,k}\widetilde{\mathbf{U}_{k}}\hat{\mathbf{H}}_{j,k}^{H}+\hat{\mathbf{H}}_{j,k}\widetilde{\mathbf{U}_{k}}\hat{\mathbf{H}}_{\theta,k}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{j,\theta}^{H}+\mbox{Tr}(\widetilde{\mathbf{U}_{k}}\mathbf{J}_{j,k}^{T})\mathbf{K}_{j,k}+\hat{\mathbf{H}}_{j,\theta}\mathbf{\Theta}\mathbf{H}_{\theta,k}\widetilde{\mathbf{U}_{k}}\hat{\mathbf{H}}_{j,k}+\hat{\mathbf{H}}_{j,\theta}\mathbf{\Theta}\mathbf{H}_{\theta,k}\widetilde{\mathbf{U}_{k}}\hat{\mathbf{H}}_{\theta,k}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{j,\theta}^{H}$
$\displaystyle+\hat{\mathbf{H}}_{j,\theta}\mathbf{\Theta}\mbox{Tr}(\widetilde{\mathbf{U}_{k}}\mathbf{J}_{\theta,k}^{T})\mathbf{K}_{\theta,k}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{j,\theta}^{H}+\mbox{Tr}(\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,k}\widetilde{\mathbf{U}_{k}}\hat{\mathbf{H}}_{\theta,k}^{H}\mathbf{\Theta}^{H}\mathbf{J}_{j,\theta}^{T})\mathbf{K}_{j,\theta}+\mbox{Tr}(\widetilde{\mathbf{U}_{k}}\mathbf{J}_{\theta,k}^{T})\mbox{Tr}(\mathbf{\Theta}\mathbf{K}_{\theta,k}\mathbf{\Theta}^{H}\mathbf{J}_{j,\theta}^{T})\mathbf{K}_{j,\theta}+\sigma_{j}^{2}\mathbf{I}.$

where $\mathbf{\Sigma}_{\overline{k}}$ and $\mathbf{\Sigma}_{\overline{j}}$
are given as in (9). This result is derived based on the identities: For any
$\mathbf{H}\sim\mathcal{CN}(\hat{\mathbf{H}},\mathbf{J}\otimes\mathbf{K})$,
there is
$\mathbb{E}[\mathbf{H}\mathbf{X}\mathbf{H}^{H}]=\hat{\mathbf{H}}\mathbf{X}\hat{\mathbf{H}}^{H}+\mbox{Tr}(\mathbf{X}\mathbf{J}^{T})\mathbf{K}$
and
$\mathbb{E}[\mathbf{H}^{H}\mathbf{X}\mathbf{H}]=\hat{\mathbf{H}}^{H}\mathbf{X}\hat{\mathbf{H}}+\mbox{Tr}(\mathbf{K}\mathbf{X})\mathbf{J}^{T}$.
Note that an upper bound for ergodic WSR can be easily derived from (9) in the
presence of ideal CSI, i.e., by setting $\mathbf{J}$ and $\mathbf{I}$ equal to
the identity.

##  III Robust Beamforming Via Expected WMMSE

The problem of ergodic WSR maximization (6) is non-convex due to interference
and in the case of perfect CSI it can be solved with an equivalent problem
formulation of EWMMSE, which exploits the relationship between the WSR and the
expected error covariance matrices [13].

###  III-A Digital Combining

We assume that the FD BS for UL user $k$ and the DL user $j$ apply the
combiners $\mathbf{F}_{k}$ and $\mathbf{F}_{j}$ to estimate their data streams
as $\hat{\mathbf{s}}_{k}=\mathbf{F}_{k}\mathbf{y}_{k}$ and
$\hat{\mathbf{s}}_{j}=\mathbf{F}_{j}\mathbf{y}_{j}$. Let
$\mathbf{E}_{\tilde{k}}$ and $\mathbf{E}_{\tilde{j}}$ denote the MSE error
matrices for instantaneous CSI for UL user $k$ and DL user $j$, respectively,
which can be written as
$\mathbf{E}_{\tilde{k}}=\mathbb{E}[\mathbf{e}_{k}\mathbf{e}_{k}^{H}],\quad\mathbf{E}_{\tilde{j}}=\mathbb{E}[\mathbf{e}_{j}\mathbf{e}_{j}^{H}]$
where $\mathbf{e}_{k}=\mathbf{F}_{k}\mathbf{y}_{k}-\hat{\mathbf{s}}_{k}$ and
$\mathbf{e}_{j}=\mathbf{F}_{j}\mathbf{y}_{j}-\hat{\mathbf{s}}_{j}$ denote the
error vectors. Let $\mathbf{Q}_{k},\mathbf{T}_{j},\mathbf{T}_{0}$ and
$\mathbf{Q}_{j,k}$ denote the matrices defined as

$\mathbf{Q}_{k}=\mathbb{E}_{\mathbf{H}|\hat{\mathbf{H}}}(\overline{\mathbf{H}}_{k}\widetilde{\mathbf{U}_{k}}\overline{\mathbf{H}}_{k}^{H}),\,\mathbf{T}_{j}=\mathbb{E}_{\mathbf{H}|\hat{\mathbf{H}}}(\overline{\mathbf{H}}_{j}\widetilde{\mathbf{V}_{j}}\overline{\mathbf{H}}_{j}^{H}),$
(10a)
$\mathbf{T}_{0}=\mathbb{E}_{\mathbf{H}|\hat{\mathbf{H}}}(\overline{\mathbf{H}}_{0}\widetilde{\mathbf{V}_{j}}\overline{\mathbf{H}}_{0}^{H}),\,\mathbf{Q}_{j,k}=\mathbb{E}_{\mathbf{H}|\hat{\mathbf{H}}}(\overline{\mathbf{H}}_{j,k}\widetilde{\mathbf{U}_{k}}\overline{\mathbf{H}}_{j,k}^{H}).$
(10b)

Since we adopt the EWMMSE approach, consider the expected MSE of the UL and DL
user as

$\displaystyle\mathbf{E}_{k}$
$\displaystyle=\mathbb{E}_{\mathbf{H}|\hat{\mathbf{H}}}[\mathbf{E}_{\tilde{k}}]$
(11a)
$\displaystyle=\mathbf{F}_{k}\mathbf{Q}_{k}\mathbf{F}_{k}^{H}-\mathbf{F}_{k}(\hat{\mathbf{H}}_{k}+\hat{\mathbf{H}}_{0,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,k})\mathbf{U}_{k}+\mathbf{F}_{k}\mathbf{T}_{0}\mathbf{F}_{k}^{H}$
$\displaystyle\quad+\sigma_{0}^{2}\mathbf{F}_{k}\mathbf{F}_{k}^{H}-\mathbf{U}_{k}^{H}\hat{\mathbf{H}}_{k}^{H}\mathbf{F}_{k}^{H}-\mathbf{U}_{k}^{H}\hat{\mathbf{H}}_{\theta,k}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{0,\theta}^{H}\mathbf{F}_{k}^{H}-\mathbf{I}.$
$\displaystyle\mathbf{E}_{j}$
$\displaystyle=\mathbb{E}_{\mathbf{H}|\hat{\mathbf{H}}}[\mathbf{E}_{\tilde{j}}]$
(11b)
$\displaystyle=\mathbf{F}_{j}\mathbf{T}_{j}\mathbf{F}_{j}^{H}-\mathbf{F}_{j}(\hat{\mathbf{H}}_{j}+\hat{\mathbf{H}}_{j,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,0})\mathbf{V}_{j}+\mathbf{F}_{j}\mathbf{Q}_{j,k}\mathbf{F}_{j}^{H}$
$\displaystyle\quad+\sigma_{j}^{2}\mathbf{F}_{j}\mathbf{F}_{j}^{H}-\mathbf{V}_{j}^{H}\hat{\mathbf{H}}_{j}^{H}\mathbf{F}_{j}^{H}-\mathbf{V}_{j}^{H}\hat{\mathbf{H}}_{\theta,0}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{j,\theta}^{H}\mathbf{F}_{j}^{H}-\mathbf{I}.$

By minimizing the trace of the (11), the optimal EWMMSE combiners can be
derived as

$\mathbf{F}_{k}=\mathbf{U}_{k}^{H}(\hat{\mathbf{H}}_{k}^{H}+\hat{\mathbf{H}}_{\theta,k}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{0,\theta}^{H})(\mathbf{Q}_{k}+\mathbf{T}_{0}+\sigma_{0}^{2}\mathbf{I})^{-1},$
(12a)
$\mathbf{F}_{j}=\mathbf{V}_{j}^{H}(\hat{\mathbf{H}}_{j}^{H}+\hat{\mathbf{H}}_{\theta,0}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{j,\theta}^{H})(\mathbf{T}_{j}+\mathbf{Q}_{j,k}+\sigma_{j}^{2}\mathbf{I})^{-1}.$
(12b)

###  III-B Active Digital Beamforming Under Imperfect CSI

Given the optimal combiners, the EWMMSE problem with respect to the digital
beamformers under the average total sum-power constraint and given the IRS
phase response $\mathbf{\Theta}$ fixed, can be formally stated as

$\underset{\begin{subarray}{c}\mathbf{V}_{j},\mathbf{U}_{k}\end{subarray}}{\min}\quad\mbox{Tr}(\mathbf{W}_{k}\mathbf{E}_{k})+\mbox{Tr}(\mathbf{W}_{j}\mathbf{E}_{j}),\vspace{-2mm}$
(13a)
$\vspace{-2mm}\text{s.t.}\;\mbox{Tr}(\mathbf{U}_{k}\mathbf{U}_{k}^{H})\preceq\alpha_{k},\;\&\;\;\mbox{Tr}(\mathbf{V}_{j}\mathbf{V}_{j}^{H})\leq\alpha_{0},$
(13b)

where $\mathbf{W}_{i}$ is a constant weight matrix associated with node $i$.
The problem (13) and the WSR maximization problem (6) are equivalent if their
gradient results to be the same, which can be assured if the weight matrices
are chosen as

$\mathbf{W}_{k}=\frac{w_{k}}{\mbox{ln}\;2}(\mathbf{E}_{k})^{-1},\quad\mathbf{W}_{j}=\frac{w_{j}}{\mbox{ln}\;2}(\mathbf{E}_{j})^{-1}.$
(14)

Such a result can be easily shown by following a similar proof provided in
[14, Appendix A], but carried out for the case of average MSE error under the
imperfect CSI case. To optimize the digital beamformers
$\mathbf{V}_{j},\mathbf{U}_{k}$ we take the partial derivative of the
Lagrangian function of (13) with respect to their conjugate, which leads to
the following optimal beamformers

$\displaystyle\mathbf{X}_{k}=$
$\displaystyle(\hat{\mathbf{H}}_{k}^{H}+\hat{\mathbf{H}}_{\theta,k}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{0,\theta}^{H})\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\hat{\mathbf{H}}_{k}+\mbox{Tr}(\mathbf{K}_{k}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k})\mathbf{J}_{k}^{T}+\hat{\mathbf{H}}_{\theta,k}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{0,\theta}^{H}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\hat{\mathbf{H}}_{0,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,k}$
(15a)
$\displaystyle+\mathbf{H}_{k}^{H}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\hat{\mathbf{H}}_{0,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,k}+\mbox{Tr}(\mathbf{K}_{\theta,k}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{0,\theta}^{H}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\hat{\mathbf{H}}_{0,\theta}\mathbf{\Theta})\mathbf{J}_{\theta,k}^{T}+\hat{\mathbf{H}}_{\theta,k}^{H}\mathbf{\Theta}^{H}\mbox{Tr}(\mathbf{K}_{0,\theta}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k})\mathbf{J}_{0,\theta}^{T}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,k}$
$\displaystyle+\hat{\mathbf{H}}_{j,k}^{H}\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j}\hat{\mathbf{H}}_{j,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,k}+(\hat{\mathbf{H}}_{j,k}^{H}+\hat{\mathbf{H}}_{\theta,k}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{j,\theta}^{H})\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j}\hat{\mathbf{H}}_{j,k}+\mbox{Tr}(\mathbf{K}_{j,k}\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j})\mathbf{J}_{j,k}^{T}$
$\displaystyle+\mbox{Tr}(\mathbf{K}_{0,\theta}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k})\mbox{Tr}(\mathbf{K}_{\theta,k}\mathbf{\Theta}^{H}\mathbf{J}_{0,\theta}^{T}\mathbf{\Theta})\mathbf{J}_{\theta,k}^{T}+\hat{\mathbf{H}}_{\theta,k}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{j,\theta}^{H}\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j}\hat{\mathbf{H}}_{j,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,k}+\mbox{Tr}(\mathbf{K}_{\theta,k}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{j,\theta}^{H}\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j}\hat{\mathbf{H}}_{j,\theta}\mathbf{\Theta})\mathbf{J}_{\theta,k}^{T}$
$\displaystyle+\hat{\mathbf{H}}_{\theta,k}^{H}\mathbf{\Theta}^{H}\mbox{Tr}(\mathbf{K}_{j,\theta}\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j})\mathbf{J}_{j,\theta}^{T}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,k}+\mbox{Tr}(\mathbf{K}_{j,\theta}\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j})\mbox{Tr}(\mathbf{K}_{\theta,k}\mathbf{\Theta}^{H}\mathbf{J}_{j,\theta}^{T}\mathbf{\Theta})\mathbf{J}_{\theta,k}^{T},$
$\displaystyle\mathbf{X}_{j}=$
$\displaystyle(\hat{\mathbf{H}}_{j}^{H}+\hat{\mathbf{H}}_{\theta,0}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{j,\theta}^{H})\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j}\hat{\mathbf{H}}_{j}+\mbox{Tr}(\mathbf{K}_{j}\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j})\mathbf{J}_{j}^{T}+\hat{\mathbf{H}}_{\theta,0}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{j,\theta}^{H}\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j}\hat{\mathbf{H}}_{j,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,0}+\hat{\mathbf{H}}_{j}^{H}\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j}\hat{\mathbf{H}}_{j,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,0}$
(15b)
$\displaystyle+\mbox{Tr}(\mathbf{K}_{\theta,0}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{j,\theta}^{H}\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j}\hat{\mathbf{H}}_{j,\theta}\mathbf{\Theta})\mathbf{J}_{\theta,0}^{T}+\hat{\mathbf{H}}_{\theta,0}^{H}\mathbf{\Theta}^{H}\mbox{Tr}(\mathbf{K}_{j,\theta}\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j})\mathbf{J}_{j,\theta}^{T}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,0}+\mathbf{H}_{0}^{H}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\hat{\mathbf{H}}_{0,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,0}$
$\displaystyle+(\hat{\mathbf{H}}_{0}^{H}+\hat{\mathbf{H}}_{\theta,0}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{0,\theta}^{H})\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\hat{\mathbf{H}}_{0}+\mbox{Tr}(\mathbf{K}_{0}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k})\mathbf{J}_{0}^{T}+\mbox{Tr}(\mathbf{K}_{j,\theta}\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j})\mbox{Tr}(\mathbf{K}_{\theta,0}\mathbf{\Theta}^{H}\mathbf{J}_{j,\theta}^{T}\mathbf{\Theta})\mathbf{J}_{\theta,0}^{T}$
$\displaystyle+\hat{\mathbf{H}}_{\theta,0}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{0,\theta}^{H}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\hat{\mathbf{H}}_{0,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,0}+\mbox{Tr}(\mathbf{K}_{\theta,0}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{0,\theta}^{H}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\hat{\mathbf{H}}_{0,\theta}\mathbf{\Theta})\mathbf{J}_{\theta,0}^{T}+\hat{\mathbf{H}}_{\theta,0}^{H}\mathbf{\Theta}^{H}\mbox{Tr}(\mathbf{K}_{0,\theta}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k})\mathbf{J}_{0,\theta}^{T}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,0}$
$\displaystyle+\mbox{Tr}(\mathbf{K}_{0,\theta}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k})\mbox{Tr}(\mathbf{K}_{\theta,0}\mathbf{\Theta}^{H}\mathbf{J}_{0,\theta}^{T}\mathbf{\Theta})\mathbf{J}_{\theta,0}^{T}.$

$\mathbf{U}_{k}=(\mathbf{X}_{k}+\lambda_{k}\mathbf{I})^{-1}(\hat{\mathbf{H}}_{k}^{H}+\hat{\mathbf{H}}_{\theta,k}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{0,\theta}^{H})\mathbf{F}_{k}^{H}\mathbf{W}_{k},$
(16a)
$\mathbf{V}_{j}=(\mathbf{X}_{j}+\lambda_{0}\mathbf{I})^{-1}(\hat{\mathbf{H}}_{j}^{H}+\hat{\mathbf{H}}_{\theta,0}^{H}\mathbf{\Theta}^{H}\hat{\mathbf{H}}_{j,\theta}^{H})\mathbf{F}_{j}^{H}\mathbf{W}_{j},$
(16b)

where $\mathbf{X}_{k}$ and $\mathbf{X}_{j}$ are defined in (15a) and (15b),
respectively, and the scalars $\lambda_{k}$ and $\lambda_{0}$ denote the
Lagrange multiplier for the uplink user $k$ and the FD BS. The multipliers can
be searched while performing power allocation for the users given the average
total sum-power constraints. Namely, consider the singular value decomposition
(SVD) of the matrices as
$\mathbf{X}_{k}=\mathbf{A}_{k}\mathbf{\Lambda}_{k}\mathbf{B}_{k}$ and
$\mathbf{X}_{j}=\mathbf{A}_{j}\mathbf{\Lambda}_{j}\mathbf{B}_{j}$, where
$\mathbf{A}_{i}$ and $\mathbf{B}_{i}$ denote the left and right unitary
matrices obtained with SVD and $\mathbf{\Lambda}_{i}$ denote the singular
values. The average power constraints (13b), after some simplifications, can
be written as

$\mbox{Tr}(\mathbf{V}_{j}\mathbf{V}_{j}^{H})=\frac{\sum_{i=i}^{M_{0}}\mathbf{S}_{j}(i,i)}{(\lambda_{0}+\mathbf{\Lambda}_{j}(i,i))^{2}},$
(17a)
$\mbox{Tr}(\mathbf{U}_{k}\mathbf{U}_{k}^{H})=\frac{\sum_{i=i}^{M_{k}}\mathbf{S}_{k}(i,i)}{(\lambda_{j}+\mathbf{\Lambda}_{k}(i,i))^{2}},$
(17b)

where the matrices $\mathbf{S}_{k}$ and $\mathbf{S}_{j}$ are defined as

$\displaystyle\mathbf{S}_{j}=\mathbf{B}_{j}(\hat{\mathbf{H}}_{j}+\hat{\mathbf{H}}_{j,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,0})^{H}$
$\displaystyle\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{W}_{j}\mathbf{F}_{j}$
(18a)
$\displaystyle(\hat{\mathbf{H}}_{j}+\hat{\mathbf{H}}_{j,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,0})\mathbf{A}_{j},$
$\displaystyle\mathbf{S}_{k}=\mathbf{B}_{k}(\hat{\mathbf{H}}_{k}+\hat{\mathbf{H}}_{0,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,k})^{H}$
$\displaystyle\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{W}_{k}\mathbf{F}_{k}$
(18b)
$\displaystyle(\hat{\mathbf{H}}_{k}+\hat{\mathbf{H}}_{0,\theta}\mathbf{\Theta}\hat{\mathbf{H}}_{\theta,k})\mathbf{A}_{k}.$

The optimal Lagrange multipliers satisfying the average total power
constraints can be searched with a linear search, and in this work, we adopt
the Bisection method. If the obtained values of the multipliers are negative,
then we replace them with zero. Note that the average power constraint is met
based on the CSI estimates.

###  III-C Passive Beamforming Under Imperfect CSI

In this section, we consider optimizing the phase response of the IRS to
jointly assist the UL and DL channels for rate enhancement. Let
$\mathbf{S},\mathbf{T}$ and $\mathbf{Z}$ denote the matrices independent of
the IRS phase response, given in Appendix A. The EWMMSE optimization problem
for the IRS phase response $\mathbf{\Theta}$ under imperfect CSI, given the
matrices $\mathbf{S},\mathbf{T}$ and $\mathbf{Z}$, can be formally stated as

$\underset{\begin{subarray}{c}\mathbf{\Theta}\end{subarray}}{\min}\quad\mbox{Tr}(\mathbf{\Theta}^{H}\mathbf{Z}\mathbf{\Theta}\mathbf{T})+\mbox{Tr}(\mathbf{\Theta}^{H}\mathbf{S}^{H})+\mbox{Tr}(\mathbf{\Theta}\mathbf{S})+c,$
(19a) $\text{s.t.}|\bm{\theta}(i)|=1,\quad\forall i,$ (19b)

where the scalar $c$ denotes the constant terms, independent of
$\mathbf{\Theta}$. In (19a), the problem is stated with respect to the matrix
$\mathbf{\Theta}$. However, we wish to maximize only the diagonal response of
such matrix to maximize the ergodic WSR or minimize the expected MSE, because
the off-diagonal elements result to be zero. Therefore, we first consider
restating the problem (19a) with respect to $\bm{\theta}$, being a vector made
of the diagonal elements of $\mathbf{\Theta}$. For doing so, we write the
trace terms as

$\mbox{Tr}(\mathbf{\Theta}^{H}\mathbf{Z}\mathbf{\Theta}\mathbf{T})=\bm{\theta}^{H}\mathbf{\Sigma}\bm{\theta},\quad\mbox{where}\;\mathbf{\Sigma}=\mathbf{Z}\odot\mathbf{T}^{T},$
(20a)
$\mbox{Tr}(\mathbf{\Theta}^{H}\mathbf{S}^{H})={\mathbf{s}}^{H}\bm{\theta}^{*},\quad\mbox{Tr}(\mathbf{\Theta}\mathbf{S})={\mathbf{s}}^{T}\bm{\theta},$
(20b)

where $\mathbf{s}$ denotes the vector containing only the diagonal elements of
the matrix $\mathbf{S}$. Problem (19a) can be restated with respect to the
vector $\bm{\theta}$ as

$\underset{\begin{subarray}{c}\bm{\theta}\end{subarray}}{\min}\quad\bm{\theta}^{H}\mathbf{\Sigma}\bm{\theta}+{\mathbf{s}_{r}}^{H}\bm{\theta}^{*}+{\mathbf{s}}^{T}\bm{\theta}_{l},\vspace{-2mm}$
(21a) $\text{s.t.}\quad|\bm{\theta}(i)|=1,\quad\forall i.\vspace{-2mm}$ (21b)

Problem (21) is still very challenging as it is a non-convex problem due to
the unit-modulus constraint. To solve it we adopt the majorization-
maximization method [15]. Such method aims to solve a more difficult problem
by constructing a series of more tractable problems stated with the upper
bound. Let $\mathcal{R}(\bm{\theta}^{(n)})$ denote the function evaluating
(21) at the $n$-th iteration for computed $\bm{\theta}$. Let
$\mathcal{R}_{u}(\bm{\theta}|\bm{\theta}^{(n)})$ denote an upper bound
constructed at $n$-th iteration for $\mathcal{R}$. According to [15], for the
problem of the form (21), an upper bound can be constructed as

$\mathcal{R}_{u}(\bm{\theta}|\bm{\theta}^{(n)})=2\mbox{Re}\\{{\mathbf{s}}^{H}\mathbf{q}^{(n)}\\}+c_{u},\vspace{-1mm}$
(22)

where $c_{u}$ denote constant terms in the upper bound independent of
$\bm{\theta}$ and $\mathbf{q}^{(n)}$ is given by

$\mathbf{q}^{(n)}=(\lambda^{max}\mathbf{I}-\mathbf{\Sigma})\bm{\theta}^{(n)}-{\mathbf{s}}^{*},$
(23)

with $\lambda^{max}$ denoting the maximum eigenvalues of $\mathbf{\Sigma}$.
Based on the result above, the hard non-convex optimization problem (21a)
simplifies to a series of the following problem

$\underset{\begin{subarray}{c}\bm{\theta}\end{subarray}}{\min}\quad
2\mbox{Re}\\{{\mathbf{s}}^{H}\mathbf{q}^{(n)}\\},\vspace{-3mm}$ (24a)
$\vspace{-3mm}\text{s.t.}\quad|\bm{\theta}(i)|=1,\quad\forall i,$ (24b)

which need to be solved iteratively until convergence for each update of
$\bm{\theta}$. By solving (24a), we get the optimal solution of $\bm{\theta}$
at $n+1$-th iteration, given the $\bm{\theta}$ at the $n$-th iteration as

$\bm{\theta}^{(n+1)}=e^{i\angle\mathbf{q}^{(n)}}.$ (25)

Formally, the optimization algorithm to optimize the IRS phase response is
given in Algorithm 1 and the joint optimization procedure is provided in
Algorithm 2.

Algorithm 1 Optimization of IRS Phase Response

Initialize: $\bm{\theta}^{0}$, iteration index $n=1$, accuracy $\epsilon$.
Evaluate: $\mathcal{R}(\bm{\theta}^{0})$.
Repeat until convergence

Calculate $\mathbf{q}_{i}^{(n)}$ with (23).

Update $\bm{\theta}_{i}^{(n+1)}$ with (25).

if
$|\mathcal{R}^{(n+1)}-\mathcal{R}(\bm{\theta}_{i}^{(n)})|/\mathcal{R}(\bm{\theta}_{i}^{(n+1)})\leq\epsilon$

Stop and return $\bm{\phi}_{i}^{(n+1)}$.

Algorithm 2 Robust Beamforming for MIMO IRS-FD system

Initialize the iteration index $n$, accuracy $\epsilon$, beamformers and
combiners.
Repeat until convergence

for $i$, where $i=k$ or $i=j$

Update $\mathbf{F}_{i}$ with (12).

Update $\mathbf{W}_{i}$ with (14).

if i=k $\rightarrow$ update $\mathbf{U}_{k}$ with (16a).

else $\rightarrow$ update $\mathbf{V}_{j}$ with (16b).

Update $\mathbf{\Theta}$ with Algorithm 2.

if convergence condition is satisfied

Stop and return the optimized variables.

Note that each update of the beamformers leads to a monotonic decrease in the
EWMMSE which assures the convergence of the proposed methods. A more formal
proof can be presented by following a similar strategy as [16]. However, due
to space limitations, we omit it here, which will be available in the extended
version of the paper.

##  IV Numerical Results

In this section, we present extensive simulation results to evaluate the
performance of the proposed robust beamforming design. We consider the FD BS
and the IRS to be centered in the positions $(0\mbox{m},0\mbox{m},0\mbox{m})$
and $(20\mbox{m},10\mbox{m},0\mbox{m})$, respectively, in the three-
dimensional coordinates. The UL and the DL user are assumed to be randomly
located in circles of radius $r=8$ m centered in the positions
$(20\mbox{m},0\mbox{m},30\mbox{m})$ and $(30\mbox{m},0\mbox{m},20\mbox{m})$,
respectively. We assume that the FD BS and the users are equipped with the
uniform linear array (ULA), placed at the distance of half-wavelength. The IRS
of $10\times 10=100$ elements is assumed to assist the MIMO FD communication
unless otherwise stated. The MIMO FD BS is assumed to have $M_{0}=15$ transmit
and $N_{0}=8$ receive antennas. The UL user $k$ and the DL user $j$ is assumed
to have $M_{k}=5$ transmit and $5$ receive antennas, respectively, and served
with $d_{k}=d_{j}=2$ data streams. The large scale path loss is modelled as
[15] and the SI channel is modelled as a Rician fading [17] with Rician factor
$1$. The direct links between the users and the FD BS are also modelled with a
Rician fading channel model with a Rician factor of $1$. The channels
involving the IRS are modelled according to Rayleigh fading [18]. We define
the transmit signal-to-noise (SNR) of our system as
$SNR=\frac{\alpha_{0}}{\sigma_{j}^{2}}=\frac{\alpha_{k}}{\sigma_{0}^{2}}$,
where $\alpha_{0}$ and $\alpha_{k}$ are the average total transmit power at
the FD BS and the multi-antenna UL user $k$. We assume that the CSI errors to
be i.i.d zero-mean circularly symmetric complex Gaussian distribution with the
same variance $\sigma_{csi}^{2}$, and therefore, the error covariance matrices
set chosen as $\mathbf{J}_{i}=\mathbf{I}$ and
$\mathbf{K}_{j}=\sigma_{csi}^{2}\mathbf{I}$, $\forall i,j$. Since the variance
of the CSI errors strictly depend on the average transmit power (used also to
estimate the channels) and the noise variance, we assume that
$\sigma_{csi}^{2}$ decays as $\mathcal{O}(SNR^{-\alpha})$, for some constant
$\alpha$, satisfying the CSI error decay rate inversely proportional to the
SNR. As the SNR represents the transmit SNR, note that $SNR\rightarrow\infty$
is equivalent to $\alpha_{k}=\alpha_{0}\rightarrow\infty$, which enhances the
CSI quality and reduces the CSI error variance as $\sigma_{csi}^{2}\rightarrow
0$. The CSI errors variance is set as $\sigma_{csi}^{2}=\rho/SNR^{\alpha}$,
where $\rho$ denotes a scale factor and $\alpha\in[0,1]$ determines the
quality of the CSI. Maintaining the quality of the CSI with $\alpha=1$ could
be exhausting in terms of resources required to acquire the CSI and therefore
we set $\alpha=0.6$.

Figure 1: Average WSR as a function of $\rho$ with SNR$=30~{}$dB.

We label our proposed design as a _FD-IRS-RB_. For comparison, we define the
following benchmark schemes: 1) _FD-IRS-Non-RB_ : FD system assisted with IRS
and non-robust beamforming, i.e., CSI is treated as perfect; 2) _FD-No-IRS-RB_
: FD system with no IRS and robust beamforming; 3) _FD-No-IRS-Non-RB_ : FD
system with no IRS and non-robust beamforming; 4)_HD-IRS-RB_ : HD system
assisted with IRS and robust beamforming; 5) _HD-IRS-Non-RB_ : HD system
assisted with IRS and non-robust beamforming; 6) _HD-No-IRS-RB_ : HD system
with no IRS and robust beamforming; 7) _HD-No-IRS-Non-RB_ : HD system with no
IRS and non-robust beamforming. We also compare our approximation proposed for
the analytical ergodic WSR derived above, which is labelled as _Analytical FD-
IRS-RB_ and compare it with the ergodic rate achieved with the EWMMSE robust
beamforming approach.

Figure 2: Average WSR as a function of SNR with $\rho=0.4$.

Fig. 1 shows the performance of the proposed robust joint beamforming design
as a function of the scale factor $\rho$ dictating the CSI error variance, by
means of Monte Carlo simulations at SNR$=30~{}$dB as a function of $\rho$. It
is shown that the proposed robust design achieves significant performance gain
compared to the naive FD-IRS-Non-RB scheme, which does not account for the CSI
errors. Moreover, we can also see that the achieved ergodic WSR accurately
matches the analytical result derived above. It is to be noted that as the CSI
error variance gets extremely large, the proposed statistically robust
beamforming design preserves significant robustness against the uncertainties
in the CSI quality.

Fig. 2 shows the ergodic WSR of the proposed design as a function of the
transmit SNR with $\rho=0.4$. We can see that our robust design achieves
significant gain in the presence of CSI errors compared to the other schemes.
Theoretically, the FD systems offer a twofold gain in the WSR compared to the
HD systems. However, Fig. 2 for the scheme FD-IRS-Non-RB, at low SNR, i.e.,
high CSI error variance, shows that the IRS-FD system may achieve no gain
compared to the IRS-aided HD system at any SNR level if non-robust beamforming
approach is adopted. This is due to the fact that the residual SI can become
extremely large in the presence of large CSI errors and can significantly
degrade the performance. However, we can observe that as the CSI error
variance decreases, the performance of the non-robust beamforming scheme tends
toward the performance of the robust beamforming schemes, which tends towards
the two-fold improvement in the WSR due to quasi-ideal simultaneous
transmission and reception.

##  V Conclusions

In this paper, a novel and statistically robust beamforming design for an IRS-
FD system under imperfect CSI is presented. A closed-form approximation for
the ergodic WSR is first derived given the statistical distribution of the
errors. Then the ergodic WSR maximization problem is handled based on the
EWMMSE method leading to two layers of sub-problems which are solved
iteratively based on alternating optimization. Simulation results show that
imperfect CSI can potentially limit the performance of the IRS-FD systems, and
adopting robust beamforming leads to significant improvement in the
performance gains.

##  Appendix A Matrices to optimize Θ

$\displaystyle\mathbf{S}=$
$\displaystyle\hat{\mathbf{H}}_{0,\theta}^{H}\mathbf{F}_{k}^{H}(\mathbf{W}_{k}\mathbf{F}_{k}\hat{\mathbf{H}}_{k}\widetilde{\mathbf{U}_{k}}\hat{\mathbf{H}}_{\theta,k}^{H}+\mathbf{W}_{k}\mathbf{F}_{k}\hat{\mathbf{H}}_{0}\widetilde{\mathbf{V}_{j}}\hat{\mathbf{H}}_{\theta,0}^{H}$
(26)
$\displaystyle+\mathbf{W}_{k}\mathbf{U}_{k}^{H}\hat{\mathbf{H}}_{\theta,k}^{H})+\hat{\mathbf{H}}_{j,\theta}^{H}\mathbf{F}_{j}^{H}(\mathbf{W}_{j}\mathbf{F}_{j}\mathbf{H}_{j}\widetilde{\mathbf{V}_{j}}\mathbf{H}_{\theta,0}^{H}$
$\displaystyle+\mathbf{W}_{j}\mathbf{F}_{j}\mathbf{H}_{j,k}\widetilde{\mathbf{U}_{k}}\mathbf{H}_{\theta,k}^{H}-\mathbf{W}_{j}\mathbf{V}_{j}^{H}\hat{\mathbf{H}}_{\theta,0}^{H}),$

$\displaystyle\footnotesize\mathbf{Z}=$
$\displaystyle\hat{\mathbf{H}}_{0,\theta}^{H}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\mathbf{H}_{0,\theta}+\hat{\mathbf{H}}_{0,\theta}^{H}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\mathbf{H}_{0,\theta}$
(27)
$\displaystyle+\mbox{Tr}(\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\mathbf{J}_{0,\theta}^{T})\mathbf{K}_{0,\theta}+\mbox{Tr}(\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\mathbf{J}_{0,\theta}^{T})\mathbf{K}_{0,\theta}$
$\displaystyle+\mbox{Tr}(\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\mathbf{J}_{0,\theta}^{T})\mathbf{K}_{0,\theta}+\hat{\mathbf{H}}_{0,\theta}^{H}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\hat{\mathbf{H}}_{0,\theta}$
$\displaystyle+\hat{\mathbf{H}}_{0,\theta}^{H}\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\hat{\mathbf{H}}_{0,\theta}+\mbox{Tr}(\mathbf{F}_{k}^{H}\mathbf{W}_{k}\mathbf{F}_{k}\mathbf{J}_{0,\theta}^{T})\mathbf{K}_{0,\theta}$
$\displaystyle+\hat{\mathbf{H}}_{j,\theta}^{H}\mathbf{F}_{j}\mathbf{W}_{j}\mathbf{F}_{j}\hat{\mathbf{H}}_{j,\theta}+\hat{\mathbf{H}}_{j,\theta}^{H}\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j}\hat{\mathbf{H}}_{j,\theta}$
$\displaystyle+\mbox{Tr}(\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j}\mathbf{J}_{j,\theta}^{T})\mathbf{K}_{j,\theta}+\mbox{Tr}(\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j}\mathbf{J}_{j,\theta}^{T})\mathbf{K}_{j,\theta}$
$\displaystyle+\hat{\mathbf{H}}_{j,\theta}^{H}\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j}\hat{\mathbf{H}}_{j,\theta}+\hat{\mathbf{H}}_{j,\theta}^{H}\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j}\hat{\mathbf{H}}_{j,\theta}$
$\displaystyle+\mbox{Tr}(\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j}\mathbf{J}_{j,\theta}^{T})\mathbf{K}_{j,\theta}+\mbox{Tr}(\mathbf{F}_{j}^{H}\mathbf{W}_{j}\mathbf{F}_{j}\mathbf{J}_{j,\theta}^{T})\mathbf{K}_{j,\theta},$

$\displaystyle\footnotesize\mathbf{T}=$
$\displaystyle\hat{\mathbf{H}}_{\theta,k}\widetilde{\mathbf{U}_{k}}\hat{\mathbf{H}}_{\theta,k}^{H}+\mbox{Tr}(\widetilde{\mathbf{U}_{k}}\mathbf{J}_{\theta,k}^{T})\mathbf{K}_{\theta,k}+\hat{\mathbf{H}}_{\theta,k}\widetilde{\mathbf{U}_{k}}\hat{\mathbf{H}}_{\theta,k}^{H}$
(28)
$\displaystyle+\mbox{Tr}(\widetilde{\mathbf{U}_{k}}\mathbf{J}_{\theta,k}^{T})\mathbf{K}_{\theta,k}+\hat{\mathbf{H}}_{\theta,0}\widetilde{\mathbf{V}_{j}}\hat{\mathbf{H}}_{\theta,0}^{H}+\mbox{Tr}(\widetilde{\mathbf{V}_{j}}\mathbf{J}_{\theta,0}^{T})\mathbf{K}_{\theta,0}$
$\displaystyle+\hat{\mathbf{H}}_{\theta,0}\widetilde{\mathbf{V}_{j}}\hat{\mathbf{H}}_{\theta,0}^{H}+\mbox{Tr}(\widetilde{\mathbf{V}_{j}}\mathbf{J}_{\theta,0}^{T})\mathbf{K}_{\theta,0}+\hat{\mathbf{H}}_{\theta,0}\widetilde{\mathbf{V}_{j}}\hat{\mathbf{H}}_{\theta,0}^{H}$
$\displaystyle+\mbox{Tr}(\widetilde{\mathbf{V}_{j}}\mathbf{J}_{\theta,0}^{T})\mathbf{K}_{\theta,0}+\hat{\mathbf{H}}_{\theta,k}\widetilde{\mathbf{U}_{k}}\hat{\mathbf{H}}_{\theta,k}^{H}+\mbox{Tr}(\widetilde{\mathbf{U}_{k}}\mathbf{J}_{\theta,k}^{T})\mathbf{K}_{\theta,k}$
$\displaystyle+\hat{\mathbf{H}}_{\theta,k}\widetilde{\mathbf{U}_{k}}\hat{\mathbf{H}}_{\theta,k}^{H}+\mbox{Tr}(\widetilde{\mathbf{U}_{k}}\mathbf{J}_{\theta,k}^{T})\mathbf{K}_{\theta,k}.$

## Acknowledgement

This work is supported by the Luxembourg National Fund (FNR)-RISOTTI–the
Reconfigurable Intelligent Surfaces for Smart Cities under Project
FNR/C20/IS/14773976/RISOTTI.

## References

  * [1] C. K. Sheemar, C. K. Thomas, and D. Slock, “Practical hybrid beamforming for millimeter wave massive MIMO full duplex with limited dynamic range,” _IEEE OJ-COMS_ , vol. 3, pp. 127–143, 2022.
  * [2] C. K. Sheemar and D. Slock, “Beamforming for bidirectional MIMO full duplex under the joint sum power and per antenna power constraints,” in _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021, pp. 4800–4804.
  * [3] P. K. Sharma and P. Garg, “Intelligent reflecting surfaces to achieve the full-duplex wireless communication,” _IEEE Communications Letters_ , vol. 25, no. 2, pp. 622–626, 2020.
  * [4] C. K. Sheemar and D. Slock, “Hybrid beamforming and combining for millimeter wave full duplex massive MIMO interference channel,” in _2021 IEEE Global Communications Conference (GLOBECOM)_. IEEE, 2021, pp. 1–6.
  * [5] C. K. Sheemar, “Hybrid beamforming techniques for massive mimo full duplex radio systems,” Ph.D. dissertation, Ph. D. dissertation, EURECOM, 2022.
  * [6] C. K. Sheemar, G. C. Alexandropoulos, D. Slock, J. Querol, and S. Chatzinotas, “Full-duplex-enabled joint communications and sensing with reconfigurable intelligent surfaces,” _arXiv preprint arXiv:2306.10865_ , 2023.
  * [7] C. K. Sheemar, S. Tomasin, D. Slock, and S. Chatzinotas, “Near-field intelligent reflecting surfaces for millimeter wave MIMO full duplex,” _arXiv preprint arXiv:2211.10700_ , 2022.
  * [8] C. K. Sheemar and D. Slock, “Receiver design and agc optimization with self interference induced saturation,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 5595–5599.
  * [9] Z. Abdullah, G. Chen, S. Lambotharan, and J. A. Chambers, “Optimization of intelligent reflecting surface assisted full-duplex relay networks,” _IEEE Wireless Communications Letters_ , vol. 10, no. 2, pp. 363–367, 2020\.
  * [10] Y. Ge and J. Fan, “Robust secure beamforming for intelligent reflecting surface assisted full-duplex MISO systems,” _IEEE Transactions on Information Forensics and Security_ , vol. 17, pp. 253–264, 2021.
  * [11] M. A. Saeidi, M. J. Emadi, H. Masoumi, M. R. Mili, D. W. K. Ng, and I. Krikidis, “Weighted sum-rate maximization for multi-IRS-assisted full-duplex systems with hardware impairments,” _IEEE Trans. Cogn. Commun. Netw._ , vol. 7, no. 2, pp. 466–481, 2021.
  * [12] Y. Rong, “Robust design for linear non-regenerative MIMO relays with imperfect channel state information,” _IEEE transactions on signal processing_ , vol. 59, no. 5, pp. 2455–2460, 2011.
  * [13] F. Negro, I. Ghauri, and D. T. Slock, “Sum rate maximization in the noisy MIMO interfering broadcast channel with partial CSIT via the expected weighted MSE,” in _IEEE ISWCS_ , 2012, pp. 576–580.
  * [14] A. C. Cirik, R. Wang, Y. Hua, and M. Latva-aho, “Weighted sum-rate maximization for full-duplex MIMO interference channels,” _IEEE Transactions on communications_ , vol. 63, no. 3, pp. 801–815, 2015.
  * [15] C. Pan, H. Ren, K. Wang, W. Xu, M. Elkashlan, A. Nallanathan, and L. Hanzo, “Multicell MIMO communications relying on intelligent reflecting surfaces,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 8, pp. 5218–5233, 2020.
  * [16] S. S. Christensen, R. Agarwal, E. De Carvalho, and J. M. Cioffi, “Weighted sum-rate maximization using weighted MMSE for MIMO-BC beamforming design,” _IEEE Transactions on Wireless Communications_ , vol. 7, no. 12, pp. 4792–4799, 2008.
  * [17] C. K. Sheemar and D. Slock, “Hybrid beamforming for bidirectional massive MIMO full duplex under practical considerations,” in _2021 IEEE 93rd Vehicular Technology Conference (VTC2021-Spring)_. IEEE, 2021, pp. 1–5.
  * [18] H. Guo, Y.-C. Liang, J. Chen, and E. G. Larsson, “Weighted sum-rate optimization for intelligent reflecting surface enhanced wireless networks,” _arXiv preprint arXiv:1905.07920_ , 2019.