# Time-Dependent AGN Disc Winds I - X-ray Irradiation

Sergei Dyda,1 Shane W. Davis,1 and Daniel Proga2
1 Department of Astronomy, University of Virginia, 530 McCormick Rd.,
Charlottesville, VA 22904, USA
2 Department of Physics & Astronomy, University of Nevada, Las Vegas, 4505 S.
Maryland Pkwy, Las Vegas, NV, 89154-4002, USA<EMAIL_ADDRESS>
###### Abstract

We study AGN line driven disc winds using time-dependent radiation
hydrodynamics. The key criterion for determining wind launching is the
coupling strength of the UV radiation field via the spectral lines of the gas.
The strength of these lines in turn relies crucially on the gas ionization
state, determined by the local X-ray intensity.

We consider a suite of models where the central ionizing radiation is affected
by scattering, absorption and re-emission by the intervening gas. In a pure
attenuation model, the disc launches an episodic wind, as previous studies
have shown. Including scattering or re-emission tends to weaken the wind,
lowering the mass flux and outflow velocity and if sufficiently dominant,
suppressing the outflow entirely. However, the exponential nature of radiative
attenuation means only a modest, factor of a few, increase in the absorption
cross section can overcome the wind suppression due to scattering and re-
emission. We find mass outflow rates of $\sim 20$% or more of the assumed
inflow rate through the disk, indicating that radiation driven winds may
significantly alter the structure of the accretion flow. The winds also supply
a large, time-varying column of material above the nominal constant disk scale
height, which will determine the geometry of reprocessed emission from the
central source. Our results suggest the need for accurate photoionization
modeling, radiation transport as well as accretion disc physics, to study
their effects on the AGN disc winds.

###### keywords:

galaxies: active - methods: numerical - hydrodynamics - radiation: dynamics

††pagerange: Time-Dependent AGN Disc Winds I - X-ray Irradiation–References

##  1 Introduction

Active Galactic Nuclei (AGN) are the central regions of galaxies that exhibit
high luminosity and activity, which is attributed to the presence of a
supermassive black hole at their centers. One of the most important features
of AGNs is the presence of powerful outflows of gas, called AGN winds, that
can reach speeds of up to tens of thousands of kilometers per second. These
winds can have a range of effects on their host galaxies. They can transport
large amounts of energy and momentum away from the central black hole, which
can regulate the growth of the black hole and influence the star formation
rate in the galaxy. The winds can also heat and ionize the gas in the
interstellar medium, which can affect the chemical composition of the gas and
the formation of stars. In some cases, the winds can even be powerful enough
to expel gas from the galaxy altogether, leading to a decrease in the overall
mass of the galaxy.

A key unanswered question is what physical processes are responsible for
launching and accelerating these strong outflows? Observations of UV
absorption lines in AGN spectra has strongly suggested the possibility that
winds are accelerated via radiation pressure on lines, so called line driving.
For line driving to successfully accelerate the wind two conditions must be
met. Firstly the gas in the AGN environment must be sufficiently ionized so
that the gas can interact with and gain momentum from the UV radiation field.
Secondly, the AGN must provide a high enough UV flux to transfer sufficient
momentum to the gas and overcome gravity.

Line driving was originally proposed as a driving mechanism for OB stars (
Lucy & Solomon (1970); Castor et al. (1975), hereafter CAK) and successfully
predicted global wind properties such as mass loss rates and outflow
velocities. These early successes lead to line driving being applied to study
outflows from other compact object system. In the case of cataclysmic
variables (CVs), disc wind models (Proga et al. (1998); Proga et al. (1999))
have been able to reproduce observed absorption profiles (Proga, 2003) and the
formation of small scale clumps (Dyda & Proga (2018a); Dyda & Proga (2018b)).
In the case of X-ray binaries, the ionization state has been observed to be
too high for line driving to be the primary acceleration mechanism, thereby
requiring an alternative such as thermal or magnetic driving (Proga & Kallman
(2002); Miller et al. (2016); Higginbottom et al. (2017); Tomaru et al.
(2020)). The presence of some optical (e.g. Muñoz-Darias et al. (2019); Mata
Sánchez et al. (2018)) and ultraviolet Castro Segura et al. (2022) wind
signatures in XRBs indicates the presence of some lower ionization material,
indicating that line driving may play some role in ejecting material.

Early photoionization studies found the acceleration due to the line force is
remarkably constant for stellar temperatures $10000\rm{K}\lesssim T\lesssim
50000$ K Abbott (1982). The physical environment around AGN is however
substantially different than in the OB star context. The gas can be more
strongly ionized, and much hotter leading to a severe reduction in the force
due to lines. Studying line driven outflows from AGN therefore requires
accurate modeling of the ionization state of the gas, the radiation transport
and the gas hydrodynamics.

The earliest attempts at photoionization modeling Stevens & Kallman (1990)
parametrized the force multiplier in terms of gas temperature and ionization.
Further efforts were made to alter the orignial CAK formulation due to changes
in gas temperature and ionization, and line saturation Gayley (1995). These
force multiplier prescriptions could then be used in analytic (Arav & Li
(1994); Arav et al. (1994)) and hydrodynamical modeling (Murray et al. (1995);
Pereyra et al. (1997); Proga et al. (2000) hereafter PSK00; Proga & Kallman
(2004) hereafter PK04; Nomura et al. (2016); Nomura & Ohsuga (2017)). More
recent attempts are in the same spirit, Dannen et al. (2019), but rather than
parametrizing results, force multiplier data tables are used by hydro
simulations directly (Dannen et al., 2023). Further they benefited from
substantially more complete data cataloguing line transitions and incident
SEDs. The conclusion from these hydrodynamical simulations is that AGN disc
winds can be launched and accelerated via line driving, provided the gas is
not over-ionized. Either the gas had to be clumpy Blondin et al. (1990), or
some intervening gas had to shield the wind from the central source (de Kool &
Begelman (1995); Murray et al. (1995)).

In order to further refine these photoionization studies, efforts were made to
understand the radiation transfer of the X-rays through the wind. Schurch et
al. (2009) used 1D radiation transfer through the PK04 solution to confirm
that the viewing angle has a significant effect on spectra due to absorption
by the wind. Subsequent development of multi-dimensional Monte-Carlo radiation
transfer codes (Sim et al. (2005); Sim et al. (2008) Sim et al. (2010)
hereafter Sim10) allowed for scattering and reprocessing effects to be taken
into account. Sim10 showed that the PK04 could reproduce complex spectral
features such as broad emission lines, narrow absorption lines and a Compton
hump. However, they found that the Monte Carlo derived ionization structure
was inconsistent with the simulation’s ionization structure. Their work was
extended by Higginbottom et al. (2014), who considered longer wavelengths and
lower ionization stages that are critical to allow effective line-driving and
the formation of broad absorption lines (BALs). They similarly concluded that
once the effects of scattering and reprocessing are included, the ionization
state corresponding to PK04 was too high to sustain a line driven wind.
Nevertheless, disc wind solutions persist in being able to reproduce observed
spectral features when using modern radiation transfer methods (see for
example Matthews et al. (2023)).

Efforts are therefore currently underway to _simultaneously_ simulate the
radiation transport, gas dynamics and photoionization modeling. Higginbottom
et al. (submitted to MNRAS) have used a Monte Carlo radiation-hydro code
(Higginbottom et al., 2020) to compute line driven disc winds from accreting
white dwarfs and found that scattering processes over-ionize the flow and
inhibit wind formation. Perliminary work indicates a similar problem will be
present in AGN simulations (N. Scepi et al., private communication). Here we
use the radiation hydrodynmaics (rad-hydro) code Athena++ to study AGN disc
winds. We use a two-band model for the radiation, where X-rays evolve
according to the time-dependent radiation transport equation whereas the UV
field is assumed to be fixed via a geometric optics approximation. The X-rays
serve to ionize the gas, which determines the effective number of optically
thick lines which couple to the UV field and extract momentum from it.

In this work, our focus is on the treatment of the X-ray band and its coupling
to the gas. We study several coupling regimes between the gas and X-ray field,
quantified via the gas opacity. We initially study a pure absorption regime as
was first considered in PSK00 and PK04 and recover their non-stationary, disc
wind solutions. We then consider scattering effects as studied in Sim10 and
show that winds are suppressed due to over-ionization. We further show that
re-radiation effects can likewise over-ionize a wind and suppress an outflow.
Finally we show that with both scattering and re-radiation operating the wind
can still be launched if the absorption opacity are a few times stronger than
the electron scattering opacity.

In the interest of focusing on the question of over/under ionization of the
gas due to X-rays we have chosen a correspondingly simplified geometric optics
treatment of the UV band. In a future study (Paper II) we will leverage the
fully multi-frequency capabilities of Athena++ (Jiang, 2022) to evolve both
the X-ray and UV fields using the time-dependent radiation transport equation.

The layout of this paper is as follows. In §2 we describe our simulation set-
up, in particular the treatment of the X-ray band which is the novel aspect of
this work. In §3 we first describe the results of our pure attenuation model
which most closely resembles the treatment of previous works. We then describe
the effects of including scattering and re-radiation. Finally we conclude in
§4 where we discuss our results, in particular in the context of other
treatments of the X-ray radiation using Monte Carlo methods and layout the
rational for our future study using time-dependent methods for the UV band.

##  2 Theory

###  2.1 Basic Equations

The basic equations for single fluid radiation hydrodynamics are

$\frac{\partial\rho}{\partial t}+\nabla\cdot\left(\rho\mathbf{v}\right)=0,$
(1a) $\frac{\partial(\rho\mathbf{v})}{\partial
t}+\nabla\cdot\left(\rho\mathbf{vv}+\sf{P}\right)=\mathbf{G}+\rho\mathbf{g}_{\rm{grav}},$
(1b) $\frac{\partial E}{\partial
t}+\nabla\cdot\left((E+P)\mathbf{v}\right)=cG^{0}+\rho\mathbf{v}\cdot\mathbf{g}_{\rm{grav}}+\rho\mathcal{L},$
(1c)

where $\rho$ is the fluid density, $\mathbf{v}$ the velocity, $\sf{P}$ a
diagonal tensor with components $P$ the gas pressure. The total gas energy is
$E=\frac{1}{2}\rho|\mathbf{v}|^{2}+\mathcal{E}$ where
$\mathcal{E}=P/(\gamma-1)$ is the internal energy and $\gamma$ the gas
constant. The radiation momentum and energy source terms are $\mathbf{G}$ and
$G^{0}$ respectively and receive contributions from the X-ray (see §2.2) and
UV bands (see §2.3). $\mathcal{L}$ is an optically thin radiative heating term
(see §2.4) and $\mathbf{g}_{\rm{grav}}$ is the gravitational acceleration (see
§2.5)

The temperature is $T=(\gamma-1)\mathcal{E}\mu m_{\rm{p}}/\rho k_{\rm{b}}$
where $\mu$ is the mean molecular weight and other symbols have their standard
meaning.

The radiation source terms $\mathbf{G}$ and $cG^{0}$ are assumed to receive
contributions from a UV and an X-ray band

$\mathbf{G}=\mathbf{G}_{\rm{UV}}+\mathbf{G}_{\rm{X}},$ (2a)
$G^{0}=G^{0}_{\rm{UV}}+G^{0}_{\rm{X}}.$ (2b)

One can alternatively think of these components as representing emission from
a compact central source (the X-ray component) and emission from a disk that
is extended in radius (the UV component). In principle, there could be X-rays
from the disk and UV from the central source, but the approximation used here
is that the ionizing radiation (X-rays) are centrally concentrated and the
radiation that potentially drives outflow comes from an extended disk alone.
We will therefore model the ionizing central source radiation (X-ray
component) via the radiation hydrodynamics module in Athena++, assuming this
component is only incident into the domain on the inner radial boundary. We
will model the disk emission (UV component) as an optically thin radiation
force computed using the modified CAK formalism.

This work is meant to be the first in a series of papers to explore the impact
of ionizing radiation on the launching of line driven winds in AGN. In order
to connect with earlier work, we have endeavored to match the assumptions and
methodological approach of PK04 as nearly as possible. This includes the
treatment of the modified CAK line driving prescription, the heating rate,
boundary conditions, and gravity, which are described below. As PK04 note,
their prescription was a conservative one in that they made a number of
unfavorable assumptions (e.g. ignoring the X-ray contribution to the radiation
force) that were not conducive to driving an outflow.

The key difference between this work and PK04 is that we leverage the
radiation hydrodynamics module of Athena++ to evolve the X-ray contribution to
the radiation field by solving the time dependent radiation transfer equation
directly. In contrast, PK04 solved the radially attenuated X-ray radiation
field with an equation of the form

$F_{x}=\exp\left(-\tau\right)\frac{L_{x}}{4\pi r^{2}},$ (3)

where $\tau=\int\rho\kappa dr$ is the optical depth integrated along radial
directed rays. We adopt a prescription that allows us to closely emulate this
pure attenuation treatment of X-ray flux while also allowing us to include
scattered and reradiated radiation, as described below.

###  2.2 X-ray Radiation

The X-ray radiation field is treated by directly solving the gray (frequency
averaged) time dependent radiation transport equation using the implicit
implementation in Athena++ (Jiang, 2021). We refer the reader to Jiang (2021)
for the specific formulation but the equation solved is equivalent to

$\frac{\partial I}{\partial t}+c\mathbf{n}\cdot\nabla I=cS_{I},$ (4)

where $I$ is the frequency integrated specific intensity and the source term

$S_{I}=\Gamma^{-3}\rho\left[\left(\kappa_{P}\frac{caT^{4}}{4\pi}-\kappa_{E}J_{c}\right)-\left(\kappa_{s}+\kappa_{F}\right)\left(I_{c}-J_{c}\right)\right],$
(5)

where $\kappa_{s}$ is the scattering opacity, $\kappa_{F}$ is the absorption
contribution to the flux mean opacity, $\kappa_{P}$ the Planck mean,
$\kappa_{E}$ the energy mean opacity and
$\Gamma=\gamma(1-\mathbf{n}\cdot\mathbf{v}/c)$ is a frame transformation
factor with $\gamma=1/\sqrt{1-v^{2}/c^{2}}$ the Lorentz factor. We use a
subcript "c" to denote variables evaluated in the comoving frame, including

$J_{c}=\frac{1}{4\pi}\int I_{c}\;d\Omega_{c},$ (6)

the angle averaged comoving frame mean intensity. The X-ray momentum and
energy source terms are then

$\mathbf{G}_{\rm{X}}=\frac{1}{c}\int\mathbf{n}S_{I}d\Omega,$ (7)
$cG^{0}_{\rm{X}}=c\int S_{I}d\Omega.$ (8)

The form of equation (5) is designed so that in the limit $v\rightarrow 0$

$\mathbf{G}_{\rm{X}}\rightarrow\frac{\rho(\kappa_{F}+\kappa_{s})}{c}\mathbf{F}_{r},$
(9) $cG^{0}_{\rm{X}}\rightarrow-c\rho(\kappa_{P}aT^{4}-\kappa_{E}E_{r}).$ (10)

For this general scheme outlined above, a common approximation for the
opacities is to set $\kappa_{E}$ equal to $\kappa_{P}$ so that radiative
equilibrium corresponds to $E_{r}=aT^{4}$. Similarly, $\kappa_{F}$ is often
set so that the sum $\kappa_{F}+\kappa_{s}$ corresponds to Rosseland mean
opacity, with $\kappa_{F}$ being the contribution from absorption opacity and
$\kappa_{s}$ is set to the electron scattering opacity $\kappa_{\rm es}$,
which is nearly constant in the Thomson limit. This scheme is usually a
sensible choice when the radiation variables represent the total frequency
integrated radiation field, but here we consider a case where the radiation
hydrodynamics solver is only being used to model the X-ray (ionizing)
component of the radiation field. Since X-rays alone do not determine the net
heating/cooling or radiation force, we need to consider an alternative
prescription that accounts for this complexity.

In our current setup, we choose to turn off the direct coupling of the X-ray
radiation field with the gas by setting
$\mathbf{G}_{\rm{X}}=cG^{0}_{\rm{X}}=0$. The impact of the radiation on
heating/cooling the gas is instead accounted for via the optically thin
prescription described in § 2.4. In principle, we could still retain the
momentum coupling given the X-ray contribution to the radiation force, but we
drop this term to be consistent with previous work. Hence, only the radiation
force from the UV disc component contributes directly to the radiation force.
Instead, the radiation transfer of the X-rays is included solely to determine
the associated ionization parameter

$\xi=\frac{4\pi cE_{r}}{n},$ (11)

where $n=\rho/\mu m_{p}$ is the gas number density with $m_{p}$ the proton
mass and $\mu=0.6$ the mean molecular mass. Since both the heating
prescription and the UV line driving force are sensitive to the value of
$\xi$, the X-ray component still indirectly impacts the radiation force and
radiative heating/cooling.

These considerations also motivate an alternative prescription for the
opacities that accounts for the fact that some fraction of the absorbed X-ray
radiation field may be thermalized and reemitted outside the X-ray (ionizing)
band. Determining the fraction reradiated in the X-rays would require a
detailed non-LTE photoionization calculation that is beyond the scope of this
work (and current numerical capabilities, at least for fully coupled radiation
hydrodynamics). We intend to explore more sophisticated approximations in
future work but for this study we consider a restricted class of models where
we parameterize reradiation as a constant fraction $f_{\sigma}$ of the
absorbed X-ray radiation that is remitted in the X-ray band. We assume this
reradiated radiation is isotropically emitted in the comoving frame.

We implement the prescription described above by setting
$\kappa_{s}=\kappa_{\rm es}$, $\kappa_{F}=\kappa_{a}$, and
$\kappa_{E}=\kappa_{a}(1-f_{\sigma})$. Here, $\kappa_{a}$ is an attenuation
opacity that we can vary to see how the radiation field responds under
different assumptions. We also set $\kappa_{P}=0$ to decouple the reemitted
radiation from gas temperature, which is instead set by our heating
prescription. Our source term then reduces to

$S_{I}=\Gamma^{-3}\rho\left[\kappa_{\rm
es}\left(J_{c}-I_{c}\right)+\kappa_{a}(f_{\sigma}J_{c}-I_{c})\right].$ (12)

With these assumptions $\kappa_{a}$ is a form of pure absorption opacity that
when $f_{\sigma}=0$, simply attenuates the radiation field. When
$f_{\sigma}>0$, a fraction $f_{\sigma}$ of this attenuated radiation is
reradiated isotropically in the comoving frame. When $f_{\sigma}=1$ our
treatment of absorption and reemission is equivalent to our treatment of
scattering.

Finally, we assume a constant radiation flux at the inner radial boundary
directed along outgoing radial rays. To mimic the effects of a very compact
source at the origin, only the most nearly radial outgoing rays are given a
non-zero intensity on the inner boundary. This is accomplished by using a
spherical polar discretization of the unit sphere (see e.g. White et al.,
2023) with a polar axis corresponds to a unit vector in the local $r$
direction and selecting the ring of ray with smallest polar angle. The
intensity of the selected rays is chosen such that the total flux in X-rays is
a fraction $f_{X}=0.10\Gamma_{D}$ of the disc Eddington fraction.

###  2.3 UV Radiation

The UV flux is assumed to be time-independent and the wind is assumed to be
optically thin to the UV continuum. The momentum transfer can be broken up
into contributions from electron scattering and radiation pressure on spectral
lines,

$\mathbf{G}_{\rm{UV}}=\frac{\rho\kappa_{\rm
es}}{c}\int\mathbf{n}I_{D}f_{UV}\left[1+M(t)\right]d\Omega,$ (13)

where the surface integral is over the disc, assumed to be radiating like a
Shakura-Sunyaev disc (Shakura & Sunyaev, 1973) with a constant accretion rate
$\dot{M}$. Even though mass is injected into the domain, there is no feedback
on the assummed mass flow rate, which remains constant. Hence, the radiation
at the disk surface is modeled as

$I_{D}=\frac{3}{\pi}\frac{GM_{BH}}{r_{*}^{2}}\frac{c}{\kappa_{\rm
es}}\Gamma_{D}\left[\left(\frac{r_{*}}{r}\right)^{3}\left(1-\left[\frac{r_{*}}{r}\right]^{1/2}\right)+\frac{x}{3\pi}\left\\{\sin^{-1}\left(\frac{r_{*}}{r}\right)-\frac{r_{*}}{r}\left(1-\left[\frac{r_{*}}{r}\right]^{2}\right)^{1/2}\right\\}\right].$
(14)

where $r_{*}=6GM/c^{2}$, $x$ is the re-radiation factor ($x=1$ in this work),
and $\Gamma_{D}=\dot{M}/\dot{M}_{\rm Edd}$, where

$\dot{M}_{\rm Edd}=\frac{4\pi GM_{BH}}{\eta\kappa_{\rm es}c}.$ (15)

We assume an efficiency $\eta=1/12$ in this work. In this model, the spectrum
is blackbody, with a temperature

$T_{D}=\left[\frac{\pi I_{D}}{\sigma}\right]^{1/4}.$ (16)

The fraction of the intensity in the UV is then found by integrating the
Planck function

$f_{UV}=\frac{\pi}{\sigma
T_{D}^{4}}\int_{\lambda_{1}}^{\lambda_{2}}\frac{2hc^{2}}{\lambda^{5}}\frac{1}{\exp\left\\{hc/\lambda
k_{b}T_{D}\right\\}-1},$ (17)

where $\lambda_{1}=200$ Å and $\lambda_{2}=3200$ Å and $\sigma$ is the Stefan-
Boltzmann constant.

The strength of the spectral lines, relative to electron scattering is
quantified via the force multiplier Owocki et al. (1988)

$M(t)=kt^{\alpha}\left[\frac{\left(1+\tau_{\rm{max}}\right)^{1-\alpha}-1}{\tau_{\rm{max}}^{1-\alpha}}\right],$
(18)

with $\alpha=0.6$ the ratio of optically thick to optically thin lines,

$k=0.03+0.385\ \exp\left\\{-1.4\ \xi^{0.6}\right\\},$ (19)

the number of lines and $\tau_{\rm{max}}=t\eta_{\rm{max}}$ with the parameter

$\log_{10}\eta_{\rm{max}}=\begin{cases}&6.9\exp\left\\{0.16\
\xi^{0.4}\right\\}\\\ &9.1\exp\left\\{-7.96\times 10^{-3}\ \xi\right\\}\\\
\end{cases}$ (20)

determining the maximum number of lines available and the optical depth
parameter

$t=\frac{\kappa_{\rm es}\rho v_{\rm{th}}}{\left|dv_{l}/dl\right|},$ (21)

with the thermal velocity of the gas $v_{\rm{th}}=4.2\times 10^{5}$ cm/s. The
disc as a source of radiation is taken to extend outside the computational
domain to $R_{d}=1500r_{g}$.

Finally, the work done by the UV radiation force is then

$cG^{0}_{\rm{UV}}=v\cdot\mathbf{G}_{\rm{UV}}.$ (22)

###  2.4 Radiative Heating

We assume a locally optically thin radiative heating and cooling due to
Compton, X-ray photoionization, Bremsstrahlung and lines. For a 10keV
Bremsstrahlung, Blondin (1994) found a heating function of the form

$\rho\Lambda=n^{2}\left(G_{\rm{Comp}}+G_{\rm{Xray}}-L_{Brem}-L_{\rm{lines}}\right),$
(23)

where the contributions in cgs units, $\left[\rm{erg\ cm^{-3}\
s^{-1}}\right]$, due to Compton, X-ray photoionization, Bremsstrahlung and
lines are respectively

$G_{\rm{Comp}}=8.9\times 10^{-36}\xi\left(T_{X}-T\right),$ (24a)
$G_{\rm{Xray}}=1.5\times 10^{-21}\xi^{1/4}T^{-1/2}\left(1-T/T_{X}\right),$
(24b) $L_{\rm{B}}=3.3\times 10^{-27}T^{1/2},$ (24c)
$\begin{split}L_{\rm{lines}}=\delta\Big{[}1.7\times
10^{-18}\exp\left\\{-\frac{T_{L}}{T}\right\\}\xi^{-1}T^{-1/2}+10^{-24}\Big{]},\end{split}$
(24d)

with the line temperature $T_{L}=1.3\times 10^{5}$ and $\delta=1$ the
parameter that controls the efficiency of line cooling.

###  2.5 Gravitational Potential

We treat the disc as a mass reservoir of fixed density and fixed radiation
intensity. Thus, we are implicitly assuming that our disc boundary is at a
height $h$ above the disc midplane. For radiation dominated discs (Shakura &
Sunyaev, 1973) the photosphere should be roughly a few times the Eddington
fraction in units of $r_{\rm{ISCO}}$. We take the shift to be

$h=3\Gamma_{D}r_{*}f_{h},$ (25)

above the disc midplane. In this work we choose $f_{h}=0.85$, where $f_{h}=1$
would correspond to a displacement by one Shakura & Sunyaev scale height. In
these shifted coordinates, the gravitational acceleration acquires additional
non-radial components

$g_{r}=\frac{GM_{\rm{BH}}}{r^{{}^{\prime}2}}\sin\alpha,$ (26a)
$g_{\theta}=\frac{GM_{\rm{BH}}}{r^{{}^{\prime}2}}\cos\alpha,$ (26b)

where

$\sin\alpha=\frac{h\sin\theta}{r^{\prime}},$ (27a)
$\cos\alpha=\frac{r+h\cos\theta}{r^{\prime}},$ (27b)

and

$r^{{}^{\prime}2}=r^{2}+2rh\cos\theta+h^{2}.$ (28)

###  2.6 Simulation Parameters

Follwoing PK04 we consider a central black hole mass $M_{BH}=10^{8}M_{\odot}$,
corresponding to an innermost stable circular orbit (ISCO)
$r_{*}=6r_{g}=8.8\times 10^{13}\rm{cm}$. We report some results in terms of
gravitational radii $r_{g}=GM/c^{2}=1.4\times 10^{13}\rm{cm}$ and the orbital
period at the inner boundary $t_{0}=2\pi((60r_{g})^{3}/GM)^{1/2}=1.44\times
10^{6}s$

We impose inflow (outflow) boundary conditions at the inner (outer) radial
boundaries and axis boundary conditions along the $\theta$ = 0 axis. We assume
a reflection symmetry about the $\theta=\pi/2$ midplane. We use a vacuum
boundary condition for the radiation along the disc midplane and outer radial
boundaries and keep the X-ray flux fixed at the inner radial boundary. After
every full time step we reset $\rho_{d}=10^{-8}\rm{gcm^{-3}}$, $v_{r}=0$ and
$v_{\phi}=v_{K}=\sqrt{GM/r}$. We also impose that the vertical velocity
component $v_{\theta}$ is unchanged due to resetting density.

We choose a domain size $n_{r}$ × $n_{\theta}$ = 96 × 140. The radial domain
extends over the range $10\ r_{*}<r<500\ r_{*}$ with logarithmic spacing
$dr_{i+1}/dr_{i}=1.05$. The polar angle range is $0<\theta<\pi/2$ and has
logarithmic spacing $d\theta_{j+1}/d\theta{j}=0.938$ which ensures that we
have sufficient resolution near the disc midplane to resolve the acceleration
of the flow.

Initially, the cells along the disc are set to have $\rho=\rho_{d}$,
$v_{r}=v_{\theta}=0$, $v_{\phi}=v_{K}$. In the rest of the domain
$\rho=10^{-20}\ \rm{gcm^{-3}}$ and $v_{r}=v_{\theta}=v_{\phi}=0$. Everywhere
the temperature is constant along vertical cylinders corresponding to the
Shakura-Sunyaev disc temperature at the base given by Shakura & Sunyaev
(1973).

We take the disc Eddington number $\Gamma_{D}=0.5$ and re-radiation factor
$x=1$ (see eq 14). This choice of $\Gamma_{D}$ corresponds to an accretion
rate $\dot{M}=1.3\;\rm M_{\odot}/yr$ for $M_{BH}=10^{8}M_{\odot}$. We assume
the luminosity of the central source radiation in X-rays is 10% of the
corresponding disk luminosity. We assume the central source does not emit in
the UV. The UV radiation from the disc is assumed to extend all the way from
the ISCO (effectively outside the simulation domain) with $r_{*}\leq R_{d}\leq
1500r_{*}$. A sphere of radius $r_{*}$, effectively the black hole and nearby
corona, is assumed to be optically thick for purposes of shielding the wind
from the backside of the disc. . We impose a density floor
$\rho_{\rm{floor}}=10^{-22}\rm{gcm^{-3}}$ which adds matter to stay above this
floor, while conserving momentum, if the density ever drops below it. In
addition, we have a temperature floor at the same cylindrical radius.

##  3 Results

Model |  | Opacity [$\kappa_{\rm{es}}$] | Wind Properties | Description
---|---|---|---|---
$\kappa_{a}$ | $\kappa_{s}$ | $f_{\sigma}$ | $\log\dot{m}\ [M_{\odot}/\rm{yr}]$ |  $v_{\rm{out}}$ [km/s] | $\omega$ | $\Delta t_{\rm{high}}\ [t_{0}]$ |
A1 | Pure attenuation | 1 | 0 | 0 | -0.5 | 11 000 | $76^{\circ}$ | – | Episodic Wind
S1 | Pure scattering | 0 | 1 | 0 | – | – | – | – | No wind
A1S1 | Attenuation & Scattering | 1 | 1 | 0 | -1.3 | 6 000 | $69^{\circ}$ | 85 | Episodic Wind
A2S1 | Attenuation & Scattering | 2 | 1 | 0 | 0.53 | 10 000 | $60^{\circ}$ | – | Stationary wind
A10S1 | Attenuation & Scattering | 10 | 1 | 0 | 1.1 | 20 000 | $70^{\circ}$ | – | Stationary wind
A1S1f50 | Scattering & 50% re-emission | 1 | 1 | 0.5 | – | – | – | – | No Wind
A2S1f50 | Scattering & 50% re-emission | 2 | 1 | 0.5 | -1.4 | 6 000 | 64 | 102 | Episodic wind
A1f5 | 5% re-emission | 1 | 0 | 0.05 | -1.7 | 9 000 | $70^{\circ}$ | 92 | Episodic wind
A1f10 | 10% re-emission | 1 | 0 | 0.10 | -1.6 | 5000 | $67^{\circ}$ | 107 | Episodic wind
Table 1: Summary of models listing the gas/radiation coupling parameters that
were varied (absorption opacity $\kappa_{a}$, scattering opacity $\kappa_{s}$
and re-emission fraction $f_{\sigma}$) and the resulting wind properties
including mass outflow rate $\dot{m}$, outflow velocity $v_{\rm{out}}$,
inclination angle $\omega$ and time between high state outflows $\Delta
t_{\rm{high}}$.

###  3.1 Absorption Effects

Figure 1: Wind solution at $t=196t_{0}$ showing density (and gas velocity
vectors and $M(t)=2$ contour in red), temperature, poloidal velocity (and unit
gas velocity), force multiplier (and radiation force vectors), ionization
parameter (and unit UV flux vectors) and X-ray radiation energy density ( and
X-ray flux). We only plot vectors in parts of the domain above a
representative density floor of $\rho=10^{-19}\ \rm{gcm^{-3}}$ for clarity. We
see a strong outflow at $\omega\sim 75^{\circ}$, with acceleration primarilly
occuring at the base of the wind where the ionization is low, $\xi\lesssim 2$
and the force multiplier correspondingly large, $M(t)\sim 10$

Figure 2: Model A1. Top panel \- Time dependent mass flux through the outer
boundary for $0^{\circ}\leq\theta\leq 85^{\circ}$. We highlight a typical low
(high) state in blue (red) where the wind is weak (strong). Center panels \-
Dynamical wind variables for the low (left) and high (right) states. The
shading shows the one standard deviation variations of the respective
quantities with respect to time. The low state is dominated by a fast, low
density component with moderate inclinations. The high state has a slower,
less dense component at high inclinations and a fast, dense component at low
inclinations that is responsible for the enhanced mass outflow rate. Bottom
panels \- Time-averaged density and velocity fields for the low (left) and
high (right) states. The high state shows evidence of a failed wind at
moderate inclination angles.

For our fiducial model, we utilize a pure attenuation ($\kappa_{a}=\kappa_{\rm
es}$ and $f_{\sigma}=0$) treatment to emulate equation (3). We first test our
results by using an angular grid with a single radial ray. This single ray
approximation is the closest we come to equation (3) with the primary
difference being that $\tau$ is accumulated over the light crossing timescale
to the current radius rather than constructed from the current instantaneous
value of the column. The wind properties obtained from this single ray
approximation were consistent with PK04 results.

Since we are ultimately interested in the effects of scattering and
reemission, we perform our simulations with $n=256$ rays uniformly covering
the unit sphere. We checked for convergence of our model using $n_{r}=128,512$
rays and concluded this was a sufficiently resolved angular resolution.

We find radial, time variable disc wind solutions, broadly in agreement with
the results in PK04. In Fig. 1 we show (going counter-clockwise) density (and
gas velocity vectors), temperature, poloidal velocity (and unit gas velocity),
force multiplier (and radiation force vectors), ionization parameter (and unit
UV flux vectors) and X-ray radiation energy density ( and X-ray flux) at
$t=196t_{0}$. On long time-scales, the outflow has two components, a faster,
less dense flow at $\theta\sim 50^{\circ}$ and a denser, slower component at
$\theta\sim 75^{\circ}$.

The fast component has a density $\rho\sim 10^{-18}\rm{g\ cm^{-3}}$ and
velocities as high as 11 000 km/s. The slow component is an order of magnitude
denser with $\rho\sim 10^{-17}\rm{g\ cm^{-3}}$ and velocities of $v\sim 5\
000$ km/s.

The wind is optically thick with $N_{H}\gtrsim 10^{24}$ for $\theta\lesssim
60^{\circ}$. We can understand the wind launching via the X-ray radiation
energy density. In the wind, $E_{r}$ is low, ensuring the ionization is also
sufficiently low and the force multiplier large enough for the line driving to
overcome gravity. From the force multiplier plot, we see the radiation force
is sufficiently strong ($M(t)>1$, indicated by the green contours) in the very
base of the disc, close to the central region. Within this contour the
ionization $\xi\lesssim 10^{4}$ is sufficiently low so as to not suppress the
force multiplier. In the limit of an optically thin wind,
$\tau_{\rm{max}}\rightarrow 0$, the maximum force multiplier $M_{\rm{max}}=2$
for $\log\xi=2.23$. For our Eddington parameter $\Gamma=0.5$, this is roughly
the minimum force multiplier to launch a wind, so our best case scenario is
wind launching in the orange colored contours of the ionization plot.

The flow is time variable and we see periods with stronger and weaker
outflows. In the top panel of Fig 2 we show the total mass flux out the outer
boundary in $0^{\circ}\leq\theta\leq 85^{\circ}$. We excised the gas close to
the disc because the disc density is quite high relative to the wind and
dominated the signal. We have identified two epochs representative phases of
the wind, a “weak" state (shaded in blue) and a “high" state (shaded in red).
For each epoch, we show the time-averaged wind fluxes exiting the outer
boundary as a function of $\theta$ (middle panels). We plot the column density
$N_{H}=\int ndr$, the momentum flux $\rho v_{r}$ and the mass flux density
$dm_{\theta}=4\pi\rho v_{r}r^{2}d\theta$. The shading shows the one standard
deviation variations of the respective quantities with respect to time. The
lower panels show the time-averaged density and velocity field.

During the low state, mass leaves the domain only at $\theta\sim 60^{\circ}$.
During the high state, the mass flux is dominated by the more radial component
where the velocities reach $\gtrsim 11\ 000\ \rm{km/s}$ but there is also a
slower, less dense component at higher latitudes. The density profiles show
that even in the high state the flow is unsteady. We see evidence of a failed
wind with matter accreting both at latitudes just above the slow part of the
wind.

Close to the disc, where the gravitational force scales as
$|\mathbf{g}_{z}|\approx GMz/r^{3}$ matter can be supported against gravity
via the radiation force due to electron scattering, which is approximately
constant near the photosphere. As the matter is launched, the relative
contribution of gravity and electron scattering changes and further
acceleration requires line driving to take over. If line driving is too weak,
due to overionization for example, the wind will fail and the launched gas
will fall back towards the central object. The wind launching is a highly
dynamic process, with gas being launched, failing, and then being swept back
up in the outflow.

In Table 1 we show a summary of our other wind models. We indicate the
gas/radiation coupling parameters that were varied (absorption opacity
$\kappa_{a}$, scattering opacity $\kappa_{s}$ and re-emission fraction
$f_{\sigma}$) and the resulting wind properties including mass outflow rate
$\dot{m}$, outflow velocity $v_{\rm{out}}$, inclination angle $\omega$ and
time between high state outflows $\Delta t_{\rm{high}}$. The goal is to
explore the effects of scattering and re-radiation on the fiducial attenuation
model. As an initial condition, we use the wind from the fiducial run at
$t=2\times 10^{8}\rm{s}\approx 140.5t_{0}$, where a wind has been well
established. We list the main radiation variables varied, namely the
absorption and scattering opacitites and the re-radiation fraction. For each
set of parameters, the solution is described as either a stationary wind, an
episodic wind or having no wind. In the case of an outflow, we list the time
averaged mass outflow rate, the peak outflow velocity and the azimuthal
position of this velocity peak. In the case of episodic winds, the wind
parameters are expressed during the “high" state during which an outflow is
occurring. For episodic winds we also list the average period between outflow
maxima.

###  3.2 Scattering Effects

Figure 3: Same as Fig 2 but for model A1S1.

We next consider cases in the pure scattering limit, where
$\kappa_{a}=\kappa_{P}=0$, $f_{\sigma}=0$ and $\kappa_{s}=\kappa_{es}$. The
radiation field changes on timescales of $\sim 10\ t_{0}$ as radiation is
scattered from the low density funnel into the higher density wind. In the
wind launching region, the X-ray energy density increases by $\sim 4$ orders
of magnitude, leading to a commensurate increase in the ionization parameter.
This effectively turns off the line driving and the wind gradually fails as
new gas does not launch and the already launched gas exits the simulation
domain.

To facilitate wind launching, we try the alternate scenario where matter is
allowed to enter the initially empty domain from the disc midplane, as we used
to launch our fiducial run. Matter enters the simulation domain, forming a
thin disc, but the X-ray energy density is even higher than when using our
restart prescription. Our results suggest that launching AGN disc winds is not
possible in a pure scattering regime.

We therefore turn our attention to a regime with both scattering and
absorption, $\kappa_{a}=\kappa_{s}=\kappa_{es}$. Unlike our fiducial model, we
find a quasi episodic wind. In the top panel of Fig 3 we plot the mass flux
exiting the simulation domain as a function of time. We identify several
outbursts and periods of quiescence in the outflow, and highlight one such low
(high) state shaded in blue (red). For these respective states we plot the
time-averaged values of the dynamical variables at the outer boundary, and the
density and velocity fields. As with the fiducial model, the low state is
characterized by a slower, less dense and higher inclination flow than the
high state.

The quasi periodicity is due to the ionization structure of the launching
region. After outburst, the wind is too ionized and gas accretes towards the
central region. This causes the density near the inner radius to peak, and the
ionization to decrease. This allows for a larger force multiplier and a new
outflow is launched. The time scale of this process is thus set by a
combination of the inflow time and outflow time of the accreting/outflowing
gas.

Increasing the absorption opacity has the effect of making the wind stationary
and strengthening the outflow. Increasing $\kappa_{a}=2\kappa_{es}$ causes the
wind to become stationary, with a fourfold increase in the mass flux compared
to the quasi episodic winds high state and a doubling of the peak velocity.
Further increasing $\kappa_{a}=10\kappa_{es}$, the values explored in PSK00,
yield winds with 15 times the mass flux and peak velocities of 20 000 km/s.

###  3.3 Re-radiation Effects

We consider models where a fraction of absorbed radiation is re-emitted in the
X-rays with $\kappa_{E}=\kappa_{a}\left(1-f_{\sigma}\right)$ and
$0<f_{\sigma}\leq 1$.

Firstly consider models with only an absorption opacity. With only 5% re-
radiation, $f_{\sigma}=0.05$, we find a very weak, episodic wind.
Qualitatively it resembles the $A1S1$ model, with slightly weaker $\sim 50\%$
lower mass flux but 50% greater peak outflow velocities. Increasing the re-
radiation to 10%, the qualitative behaviour is unchanged and we again find a
weak, episodic outflow with peak velocities of $\sim 5\ 000\rm{km/s}$.
Unsurprisingly, allowing gas to re-radiate has a similar physical effect to
the wind as scattering. The X-rays tend to isotropize, and penetrate into the
wind acceleration region. This tends to weaken the wind as gas becomes over-
ionized. Further increasing the re-radiation fraction eventually leads to the
wind being fully suppressed, which we see for $f_{\sigma}=0.5$.

We now consider models with both scattering, and absorption opacity. We set
the re-radiation to 50%. When $\kappa_{a}=\kappa_{s}=\kappa_{es}$ we find that
no wind is launched. This is expected as the model with these parameters and
no scattering failed to launch a wind, so scattering could serve only to
further inhibit wind acceleration.

We therefore increase the absorption opacity to $\kappa_{a}=2\kappa_{es}$ and
find that we again enter a regime of episodic outflows. The wind is both
quantitatively and qualitatively similar to the model without re-emission,
A2S1.

##  4 Discussion

Figure 4: Time averaged radiation field for different disc wind models. The
grey contour indicates where $M(t)=2$, where we expect the radiation force due
to lines to overcome gravity.The radiation flux is only plotted for regions
with $\rho>10^{-17}\ \rm{gcm^{-3}}$ Models with an extended region of low
X-ray radiation density (A1, A2S1) are able to sustain stronger. Winds with
only a small such region (A1f5, A1S1, A2S1f50) can only launch episodic winds
when failed winds shield the innermost part of the disc. Model S1 is always
over-ionized and cannot launch winds.

We have used time-dependent radiation transfer to reproduce well established
results for AGN disc winds that accounted only for radiative attenuation. We
then use our more sophisticated radiation treatment to study the effects of
scattering and re-radiation. We have used an idealized treatment of radiative
absorption, scattering and reprocessing to develop an intuition for their
relative effects on radiation driven disc winds.

Our study demonstrates that radiative absorption of X-ray photons must be
sufficiently strong to prevent over-ionization and allow wind launching. If
scattering and re-radiation effects are neglected, absorption opacities equal
to electron scattering is sufficient. In a pure scattering regime no wind is
launched. If scattering and absorption opacities are equal, we find a much
weaker, episodic outflow. The time scale of the variability is set by the
dynamical time for the wind to fail, fall back towards the black hole and
shield the outer parts of the disc to allow for a new outflow to be launched.
If we allow re-radiation within the X-ray band, a modest 5% re-radiation
fraction makes the wind become episodic as in the scattering case. If re-
radiation fraction is 50% we require an enhancement of a few times electron
scattering to the absorption opacity to maintain a wind.

These results show that absorption and re-processing of X-rays is crucial to
launch an AGN line driven disc wind. Any process which enhances the ionization
parameter by increasing the X-ray energy density such as scattering or re-
emission will tend to weaken and ultimately suppress the wind. In Fig. 4 we
plot the time-averaged X-ray energy density and radiation flux in units of the
central object $F_{0}=L_{BH}/4\pi r_{0}^{2}$. The grey contour indicates the
$M(t)=2$ surface, interior to which we expect the force due to line driving
$f_{LD}\sim\Gamma M(t)>f_{\rm{grav}}$. We see that the models with stationary
winds, A1 and A2S1 have a geometrically large region where the radiation force
dominates, whereas models with periodic outflows (A1f5, A1S1, A2S1f50) have
only a geometrically small region. The model with no wind, S1 has no such
region. This region of enhanced launching has low radiation energy density,
consistent with a low ionization parameter.

Fortunately, attenuation effects scale like $\exp\left(-\rho\kappa dr\right)$,
so a modest, factor of a few, increase in the opacity can exponentially
suppress X-ray photons. Thus, while we find that for $\kappa_{a}=\kappa_{es}$
a mere 5% re-radiation effect causes the wind to become weak and episodic and
50% re-radiation suppresses the wind, simply increasing the absorption opacity
to $\kappa_{a}=2\kappa_{es}$ allows for an episodic wind with 50% re-
radiation. Photoionization studies find that a boost of a factor of a few
above electron scattering is reasonable in our temperature and ionization
ranges, suggesting that more self consistent treatments of the opacity will
yield robust wind solutions.

We have assumed here the wind is optically thin in the UV and correspondingly
the UV radiation field is time-independent. This assumption may be violated by
both attenuation and re-radiation effects. We see that parts of the flow are
optically thick, and therefore obscuration of the outer parts of the wind may
be important. PSK00 had shown that accounting for UV attenuation in the
approximation that the radiation is attenuated by radial columns of gas did
not qualitatively change their solutions. However, optical depth effects of
the UV continuum should not be too challenging to account for using our new
methods. Also we expect that X-rays, attenuated by the gas, may re-radiate in
the UV, and enhance the driving flux. Further, as some of our solutions are
episodic, correctly accounting for light travel time effects may be important.

Reprocessing effects may also be important. Assuming local thermodynamic
equilibrium, gas at temperature $T$ will emit radiation in a Planck
distribution. We may therefore account for the number of UV and X-ray photons
emitted by the gas. The UV photons will change the radiation force by changing
the local radiation flux. Likewise, the X-ray photons will alter the
ionization state. As this work shows, the success of line driving for
launching outflows is crucially related to the gas ionization so determining
the ionization state and force multiplier are key.

Monte Carlo simulations in the context of CVs have shown that ionization
effects can severely inhibit line driving, suppressing or even preventing a
wind (Sim10, Higginbottom et al. (2014)). Their work solved for the radiation
field, assuming that the hydrodynamics was fixed. Here, we find that, indeed
the scattering can weaken a line driven wind, but the wind is still launched
and accelerated due to its ability to adjust in a time-dependent fashion to a
given radiation environment. A key lesson from their work is that suppression
of the force multiplier is not due only to X-rays but also the FUV
(Higginbottom et al. submitted to MNRAS). It highlights the non-linearity of
these processes and the limitation of using the mean photon energy i.e
ionization parameter approach. The shape of the SED, as also demonstrated in
photoionization studies Dannen et al. (2019) plays a crucial role in
determining the force multiplier. Further, close to the threshold luminosities
where line driving can operate as an acceleration mechanism, $\Gamma\sim
1/M_{\rm{max}}$, the mass outflow rate is strongly dependent on the maximum
value of the force multiplier. This work emphasises the need to correctly
model the SED, and hence the radiation transfer, throughout the hydro
simulation, while simultaneously computing the effects of this SED on the bulk
gas properties such as force multiplier and heating and cooling.

Our study further highlights the need for a more self-consistent disc
modeling. We have assumed a constant accretion rate of $1.3M_{\odot}$/yr, and
generate non-stationary winds with $\gtrsim 10\%$ this mass flux. This should
have a non-negligible feedback on the disc, both in terms of changing the
local accretion rate and altering the local radiation flux. Thus far,
simulations have relied on semi-analytic Ganguly & Proga (2020), Kirilov et
al. (2023) or iterative methods Nomura et al. (2020) to account for this wind
driven mass loss. Such large mass loss rates may have important implications
for the accretion disk structure, possibly contributing to the mismatch
between observed spectra and predictions of standard thin disk models (Slone &
Netzer, 2012; Laor & Davis, 2014; Nomura et al., 2020).

Furthermore the contours of $E_{r}$ shown in Figure 4 indicate that the
optical depth in the radial direction near the disk is large when the outflow
is present. This will impact the self-irradiation of the disk and lead to a
different reprocessing geometry than if the wind material was absent. This may
help to explain the relatively large sizes inferred for optical/UV emission
regions from continuum reverberation mapping.

In the future, we will build on these models by incorporating more realistic
treatment of absorption, scattering and reprocessing. Using photoionization
modeling, we can generate absorption opacity tables as a function of
temperature and ionization parameter. The relative importance of radiative
processes can then be accounted for in a self-consistent fashion, rather than
the parameter survey type treatment we have used in this work.

## Acknowledgments

Support for this work was provided by the National Aeronautics and Space
Administration under TCAN grant 80NSSC21K0496. We thank Tim Kallman, Jim
Stone, Yan-Fei Jiang and the entire DAWN TCAN collaboration for fruitful
discussions. Further thanks to Nick Higginbottom, Nico Scepi, Christian
Knigge, James Mattherws and Knox Long for sharing their recent Monte Carlo
work. This work made use of Camber for cloud-based MPI provisioning.

## References

  * Abbott (1982) Abbott D. C., 1982, ApJ, 259, 282
  * Arav & Li (1994) Arav N., Li Z.-Y., 1994, ApJ, 427, 700
  * Arav et al. (1994) Arav N., Li Z.-Y., Begelman M. C., 1994, ApJ, 432, 62
  * Blondin (1994) Blondin J. M., 1994, ApJ, 435, 756
  * Blondin et al. (1990) Blondin J. M., Kallman T. R., Fryxell B. A., Taam R. E., 1990, ApJ, 356, 591
  * Castor et al. (1975) Castor J. I., Abbott D. C., Klein R. I., 1975, ApJ, 195, 157
  * Castro Segura et al. (2022) Castro Segura N., et al., 2022, Nature, 603, 52
  * Dannen et al. (2019) Dannen R. C., Proga D., Kallman T. R., Waters T., 2019, ApJ, 882, 99
  * Dannen et al. (2023) Dannen R., Proga D., Waters T., 2023, arXiv e-prints, p. arXiv:2306.04063
  * Dyda & Proga (2018a) Dyda S., Proga D., 2018a, MNRAS, 475, 3786
  * Dyda & Proga (2018b) Dyda S., Proga D., 2018b, MNRAS, 478, 5006
  * Ganguly & Proga (2020) Ganguly S., Proga D., 2020, ApJ, 890, 54
  * Gayley (1995) Gayley K. G., 1995, ApJ, 454, 410
  * Higginbottom et al. (2014) Higginbottom N., Proga D., Knigge C., Long K. S., Matthews J. H., Sim S. A., 2014, ApJ, 789, 19
  * Higginbottom et al. (2017) Higginbottom N., Proga D., Knigge C., Long K. S., 2017, ApJ, 836, 42
  * Higginbottom et al. (2020) Higginbottom N., Knigge C., Sim S. A., Long K. S., Matthews J. H., Hewitt H. A., Parkinson E. J., Mangham S. W., 2020, MNRAS, 492, 5271
  * Jiang (2021) Jiang Y.-F., 2021, ApJS, 253, 49
  * Jiang (2022) Jiang Y.-F., 2022, ApJS, 263, 4
  * Kirilov et al. (2023) Kirilov A., Dyda S., Reynolds C. S., 2023, MNRAS, 520, 44
  * Laor & Davis (2014) Laor A., Davis S. W., 2014, MNRAS, 438, 3024
  * Lucy & Solomon (1970) Lucy L. B., Solomon P. M., 1970, ApJ, 159, 879
  * Mata Sánchez et al. (2018) Mata Sánchez D., et al., 2018, MNRAS, 481, 2646
  * Matthews et al. (2023) Matthews J. H., et al., 2023, MNRAS,
  * Miller et al. (2016) Miller J. M., et al., 2016, ApJ, 821, L9
  * Muñoz-Darias et al. (2019) Muñoz-Darias T., et al., 2019, ApJ, 879, L4
  * Murray et al. (1995) Murray N., Chiang J., Grossman S. A., Voit G. M., 1995, ApJ, 451, 498
  * Nomura & Ohsuga (2017) Nomura M., Ohsuga K., 2017, MNRAS, 465, 2873
  * Nomura et al. (2016) Nomura M., Ohsuga K., Takahashi H. R., Wada K., Yoshida T., 2016, PASJ, 68, 16
  * Nomura et al. (2020) Nomura M., Ohsuga K., Done C., 2020, MNRAS, 494, 3616
  * Owocki et al. (1988) Owocki S. P., Castor J. I., Rybicki G. B., 1988, ApJ, 335, 914
  * Pereyra et al. (1997) Pereyra N. A., Kallman T. R., Blondin J. M., 1997, ApJ, 477, 368
  * Proga (2003) Proga D., 2003, ApJ, 585, 406
  * Proga & Kallman (2002) Proga D., Kallman T. R., 2002, ApJ, 565, 455
  * Proga & Kallman (2004) Proga D., Kallman T. R., 2004, ApJ, 616, 688
  * Proga et al. (1998) Proga D., Stone J. M., Drew J. E., 1998, MNRAS, 295, 595
  * Proga et al. (1999) Proga D., Stone J. M., Drew J. E., 1999, MNRAS, 310, 476
  * Proga et al. (2000) Proga D., Stone J. M., Kallman T. R., 2000, ApJ, 543, 686
  * Schurch et al. (2009) Schurch N. J., Done C., Proga D., 2009, ApJ, 694, 1
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33
  * Sim et al. (2005) Sim S. A., Drew J. E., Long K. S., 2005, MNRAS, 363, 615
  * Sim et al. (2008) Sim S. A., Long K. S., Miller L., Turner T. J., 2008, MNRAS, 388, 611
  * Sim et al. (2010) Sim S. A., Proga D., Miller L., Long K. S., Turner T. J., 2010, MNRAS, 408, 1396
  * Slone & Netzer (2012) Slone O., Netzer H., 2012, MNRAS, 426, 656
  * Stevens & Kallman (1990) Stevens I. R., Kallman T. R., 1990, ApJ, 365, 321
  * Tomaru et al. (2020) Tomaru R., Done C., Ohsuga K., Odaka H., Takahashi T., 2020, MNRAS, 494, 3413
  * White et al. (2023) White C. J., Mullen P. D., Jiang Y.-F., Davis S. W., Stone J. M., Morozova V., Zhang L., 2023, ApJ, 949, 103
  * de Kool & Begelman (1995) de Kool M., Begelman M. C., 1995, ApJ, 455, 448


solutionSolutionsolutionfile solutionfile[Papazafeiropoulos_WCL2021-0832.R2]

# Asymptotic Analysis of Max-Min Weighted SINR for IRS-Assisted MISO Systems
with Hardware Impairments

Anastasios Papazafeiropoulos, Cunhua Pan, Ahmet Elbir, Van Nguyen, Pandelis
Kourtessis, Symeon Chatzinotas A. Papazafeiropoulos is with the Communications
and Intelligent Systems Research Group, University of Hertfordshire, Hatfield
AL10 9AB, U. K., and with SnT at the University of Luxembourg, Luxembourg. C.
Pan is with the School of Electronic Engineering and Computer Science at Queen
Mary University of London, London E1 4NS, U.K. A. M. Elbir is with the
Department of Electrical and Electronics Engineering, Koc University,
Istanbul, Turkey, and SnT at the University of Luxembourg, Luxembourg. P.
Kourtessis is with the Communications and Intelligent Systems Research Group,
University of Hertfordshire, Hatfield AL10 9AB, U. K. V.-D. Nguyen and S.
Chatzinotas are with the SnT at the University of Luxembourg, Luxembourg.
E-mail<EMAIL_ADDRESS>This work was supported by Luxembourg National
Research Fund (FNR) under the CORE project RISOTTI C20/IS/14773976 and by the
University of Hertfordshire’s 5–year Vice Chancellor’s Research Fellowship.

###### Abstract

We focus on the realistic maximization of the uplink minimum signal-to-
interference-plus-noise ratio (SINR) of a general multiple-input single-output
(MISO) system assisted by an intelligent reflecting surface (IRS) in the large
system limit accounting for HIs. In particular, we introduce the HIs at both
the IRS (IRS-HIs) and the transceiver HIs (AT-HIs), usually neglected despite
their inevitable impact. Specifically, the deterministic equivalent analysis
enables the derivation of the asymptotic weighted maximum-minimum SINR with
HIs by jointly optimizing the HIs-aware receiver, the transmit power, and the
reflect beamforming matrix (RBM). Notably, we obtain the optimal power
allocation and reflect beamforming matrix with low overhead instead of their
frequent necessary computation in conventional MIMO systems based on the
instantaneous channel information. Monte Carlo simulations verify the
analytical results which show the insightful interplay among the key
parameters and the degradation of the performance due to HIs.

###### Index Terms:

Intelligent reflecting surface, hardware impairments, massive MIMO systems,
deterministic equivalents, beyond 5G networks.

##  I Introduction

Intelligent reflecting surfaces (IRSs), consisted of low-cost, passive
reflecting elements with adjustable phase shifts, have been recognized as a
promising solution for enhancing the spectral and energy efficiency of
wireless systems [1]. Notably, a significant amount of research has been
devoted to IRS-aided systems [2, 3, 4, 5, 6]. For example, in [2], a
minimization of the transmit power at the base station (BS) subject to
individual signal-to-interference-plus-noise ratio (SINR) constraints took
place to address the transmit and reflect beamforming (RB) optimization
problem. Also, in [4], the optimum linear precoder (OLP) was studied in the
large number of antennas regime.

However, the majority of works in IRS-aided systems have relied on the highly
unrealistic assumption of ideal hardware, while practical implementations of
IRSs require taking into account the unavoidable residual transceiver hardware
impairments (T-HIs) whose omission may result in misleading design
conclusions. Especially, a cost-attractive implementation of massive multiple-
input multiple-output (mMIMO) systems suggests the use of cheap hardware,
resulting in more severe HIs [7, 8, 9]. One major category of T-HIs, known as
additive transceiver HIs (AT-HIs), can be modeled as additive Gaussian
distributed noise by accounting for the accumulated effect from all individual
HIs such as the in-phase/quadrature-phase imbalance [10, 7, 8, 9]. Another
interesting type of HIs, called IRS-HIs, emerges in IRS-assisted systems
because of the incapability of infinite precision at the IRS phase shifts [11,
12, 13]. Hence, the performance analysis of the IRS-aided systems should
include the impact of both AT-HIs and IRS-HIs.

Recently, the study of HIs on IRSs has attracted significant interest [14, 15,
16, 17, 18, 13]. The authors in [14] focused on the achievable rate while they
considered only single-input and single-output (SISO) channels. In [15], upper
bounds on the channel capacities were obtained while relying on the assumption
of no correlation among the columns of the IRS, assuming a single user
communication, and not providing closed forms regarding these bounds. In [16],
no small-scale fading was assumed and just a single user equipment (UE) was
considered. The latter limiting design setting was also assumed in [17] to
maximize the received signal-to-noise ratio (SNR) while, in [18], the authors
focused on the maximization of the secrecy rate for a finite number of BS
antennas. In this direction, in [13], we studied the impact of HIs on the
achievable rate in a multi-user setting for a finite number of BS antennas
without optimizing the transmit power.

In this paper, we make a substantial leap beyond previous works by accounting
for both AT-HIs and IRS-HIs in terms of the deterministic equivalent (DE)
analysis by formulating a max-min weighted SINR problem in the large antenna
regime. 111There are many differences between our work and existing works,
e.g., [17]. Therein, first, just the SNR was studied for a single UE scenario
while we focus on the SINR accounting for multiple UEs. Second, we consider
the max-min optimization problem, while [17] considered the rate maximization
problem. Third, [17] considered the phase shift design based on instantaneous
CSI, while our work designs the phase shift based on the statistical CSI.
Fourth, we consider the asymptotic case when the number of transmit antennas
is infinite, while [17] is suitable for a limited number of antennas due to
its high complexity when the number of transmit antennas is large.Contrary to
[4], we introduce both AT-HIs and IRS-HIs, and focus on the uplink instead of
the downlink. Therein, an already obtained OLP, based on [19], was applied
which limits the analysis while our methodology, taking HIs into
consideration, is more general in terms of DEs and optimization as the
following analysis reveals, which adds to the novelty of this work. For
example, we obtain the optimal decoder and the optimal power allocation with
HIs by following another approach than [19] that was based on uplink-downlink
duality. Also, we have taken into account the direct channel while its
manipulation was not possible by the analysis in [4]. It should be noted that
the introduction of HIs, increasing the complexity/difficulty, requires
substantial manipulations. Specifically, by considering the general realistic
scenario of correlated Rayleigh fading channels with HIs, we obtain the
optimal HIs-aware linear minimum mean square error (LMMSE) receiver and the
corresponding asymptotic optimal weighted SINR. We achieve optimal power
allocation and an RB design that require only large-scale channel statistics
and do not depend on small-scale fading changing at the order of milliseconds
that would result in prohibitively high overhead. The results allow shedding
light on the impact of HIs on such systems towards their realistic evaluation.

##  II System Model

We consider an IRS-aided multi-user massive MIMO system, where a multi-antenna
BS with $M$ antennas communicates with $K$ single-antenna UEs. To focus on the
impact of hardware distortions, we assume perfect channel state information
(CSI), which allows more direct mathematical manipulations. Hence, the results
play the role of upper bounds of practical scenarios with imperfect CSI. Note
that the CSI could be assumed perfectly known when the coherence intervals are
sufficiently long. The extension to the imperfect CSI scenario, which is of
practical importance, is the topic of future work. In particular, the IRS is
assumed in the line-of-sight (LoS) of the BS and includes $N$ passive
reflecting elements introducing shifts on the phases of the impinging waves.
Also, the IRS is controlled by the BS by means of a backhaul link.

We rely on a block-fading channel model with fixed channels in each time-
frequency coherence block but with independent realizations in each block.
Specifically,
${\mathbf{H}}_{1}=[{\mathbf{h}}_{1,1}\ldots,{\mathbf{h}}_{1,N}]\in\mathbb{C}^{M\times
N}$, ${\mathbf{h}}_{\mathrm{d},k}\in\mathbb{C}^{M\times 1}$, and
${\mathbf{h}}_{2,k}\in\mathbb{C}^{N\times 1}$ express the LoS channel between
the BS and IRS, the direct channel between the BS and UE $k$, and the channel
between the IRS and UE $k$. The vector ${\mathbf{h}}_{1,i}$ for $i=1,\ldots,N$
corresponds to the $i$th column vector of ${\mathbf{H}}_{1}$. Notably, we
consider spatial correlation instead of independent Rayleigh fading assumed in
the majority of previous works, e.g., [2]. Hence, we have

$\displaystyle{\mathbf{h}}_{\mathrm{d},k}$
$\displaystyle=\sqrt{\beta_{\mathrm{d},k}}{\mathbf{R}}_{\mathrm{BS},k}^{1/2}{\mathbf{z}}_{\mathrm{d},k},$
(1) $\displaystyle{\mathbf{h}}_{2,k}$
$\displaystyle=\sqrt{\beta_{2,k}}{\mathbf{R}}_{\mathrm{IRS},k}^{1/2}{\mathbf{z}}_{2,k},$
(2)

where ${\mathbf{R}}_{\mathrm{BS},k}\in\mathbb{C}^{M\times M}$ with
$\tr\left({\mathbf{R}}_{\mathrm{BS},k}\right)=M$ and
${\mathbf{R}}_{\mathrm{IRS},k}\in\mathbb{C}^{N\times N}$ with
$\tr\left({\mathbf{R}}_{\mathrm{IRS},k}\right)=N$ express the deterministic
Hermitian-symmetric positive semi-definite correlation matrices at the BS and
the IRS, respectively. Also, $\beta_{\mathrm{d},k}$ and $\beta_{2,k}$ are the
path-losses of the BS-UE $k$ and IRS-UE $k$ links. Note that the correlation
matrices and the path-losses are assumed to be known by practical methods,
e.g., see [20]. In addition,
${\mathbf{z}}_{\mathrm{d},k}\sim\mathcal{CN}\left({\mathbf{0}},{\bm{\mathrm{I}}}_{M}\right)$
and
${\mathbf{z}}_{2,k}\sim\mathcal{CN}\left({\mathbf{0}},{\bm{\mathrm{I}}}_{N}\right)$
express the respective fast-fading vectors. Moreover, we assume that
${\mathbf{H}}_{1}$ is a full rank channel matrix described as

$\displaystyle[{\mathbf{H}}_{1}]_{m,n}$ $\displaystyle=\sqrt{\beta_{1}}$,

where $\lambda$ and $\beta_{1}$ are the carrier wavelength and the path-loss
between the BS and IRS, while $d_{\mathrm{BS}}$ and $d_{\mathrm{IRS}}$ are the
inter-antenna separation at the BS and inter-element separation at the IRS,
respectively [4]. In addition, $\theta_{1,n}$, $\psi_{1,n}$ express the
elevation and azimuth LoS angles of departure (AoD), respectively at the BS
with respect to IRS element $n$, while $\theta_{2,n}$ and $\psi_{2,n}$ are the
elevation and azimuth LoS angles of arrival (AoA) at the IRS. The design of
${\mathbf{H}}_{1}$ could be realized by several techniques as suggested in
[21]. In addition, the response of the IRS elements is described by the
diagonal matrix
${\bm{\Phi}}=\mathrm{diag}\left(\alpha\exp\left(j\phi_{1}\right),\ldots,\alpha\exp\left(j\phi_{N}\right)\right)\in\mathbb{C}^{N\times
N}$, where $\phi_{n}\in\left[0,2\pi\right],n=1,\ldots,N$ and $\alpha\in(0,1]$
express the phase shifts applied by the IRS elements and the independent
amplitude reflection coefficient, respectively.222Recently, it was shown that
the amplitude and phase responses are intertwined [22, 23], which suggests an
interesting idea for extension of the current work, i.e., to study the impact
of active (additive transceiver distortion) and passive (IRS phase noise) HIs
by accounting for this intertwinement.

####  II-1 IRS-HIs

Since it is not possible to configure the IRS elements with infinite
precision, phase errors are introduced [11]. These IRS-HIs can be described by
means of a random diagonal phase error matrix consisting of $N$ random phase
errors, i.e.,
$\widetilde{{\bm{\Phi}}}=\text{diag}\left(\right)\left(e^{j\tilde{\phi}_{1}},\ldots,e^{j\tilde{\phi}_{N}}\right)\in\mathbb{C}^{N\times
N}$ with $\tilde{\phi}_{i},i=1,\ldots,N$ being the random phase errors of the
IRS phase shifts that are i.i.d. randomly distributed in $[-\pi,\pi)$ and
based on a certain circular distribution. 333The probability density function
(PDF) of $\tilde{\theta}_{i}$ is assumed symmetric with its mean direction
equal is zero, i.e.,
$\arg\left(\mathbb{E}[\mathrm{e}^{j\tilde{\theta}_{i}}]\right)=0$ [11].
Hence, the channel vector between the BS and UE $k$ is written as
${\mathbf{h}}_{k}={\mathbf{h}}_{\mathrm{d},k}+{\mathbf{H}}_{1}{\bm{\Phi}}\tilde{{\bm{\Phi}}}{\mathbf{h}}_{2,k},\in\mathbb{C}^{M\times
1}$, distributed as ${\cal C}{\cal N}\left(0,{\mathbf{R}}_{k}\right)$, where
${\mathbf{R}}_{k}=\beta_{\mathrm{d},k}{\mathbf{R}}_{\mathrm{BS},k}+\beta_{2,k}{\mathbf{H}}_{1}{\bm{\Phi}}\tilde{{\mathbf{R}}}_{\mathrm{IRS},k}{\bm{\Phi}}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{H}}_{1}^{{\scriptscriptstyle\mathsf{H}}}$
with
$\tilde{{\mathbf{R}}}_{\mathrm{IRS},k}=m^{2}{\mathbf{R}}_{\mathrm{IRS},k}+\left(1-m^{2}\right){\bm{\mathrm{I}}}_{N}$
and $m$ denoting its characteristic function (CF) [13, Eq. 12].

Examples of PDFs that could describe the phase noise on IRSs are the uniform
and the Von Mises distributions [11]. The former expresses completely lack of
knowledge and has CF equal to $0$ while the latter has a zero-mean and
concentration parameter $\kappa_{\tilde{\theta}}$, capturing the accuracy of
the estimation. Its CF is
$m=\frac{\mathrm{I}_{1}\\!\left(\kappa_{\tilde{\theta}}\right)}{\mathrm{I}_{0}\\!\left(\kappa_{\tilde{\theta}}\right)}$,
where $\mathrm{I}_{p}\\!\left(\kappa_{\tilde{\theta}}\right)$ is the modified
Bessel function of the first kind and order $p$.

###### Remark 1

If $m=0$ (uniform distribution), we obtain
$\widetilde{{\mathbf{R}}}_{\mathrm{IRS},k}={\bm{\mathrm{I}}}_{N}$, which means
that ${\mathbf{R}}_{k}$ does not depend on the phase shifts and the system
cannot be optimized due to the IRS. Especially, we have
${\mathbf{R}}_{k}=\beta_{\mathrm{d},k}{\mathbf{R}}_{\mathrm{BS},k}+\beta_{2,k}{\mathbf{H}}_{1}{\mathbf{H}}_{1}^{{\scriptscriptstyle\mathsf{H}}}$.
In this case, no knowledge of ${\mathbf{R}}_{\mathrm{IRS},k}$ is required at
the BS. This result is obtained also if
${{\mathbf{R}}}_{\mathrm{IRS},k}={\bm{\mathrm{I}}}_{N}$, which means that in
the case of no IRS correlation, the IRS cannot be optimized if statistical CSI
is considered. However, if the phase errors follow any other circular PDF,
${\mathbf{R}}_{k}$ is phase-dependent and the presence of IRS can be
exploited.

####  II-2 AT-HIs

Disregarding most works in the IRS literature assuming ideal transceiver
hardware, in practice, HIs, remain despite the use of any mitigation
algorithms and affect both the transmit and receive signals. In this
direction, we account for additive HIs at both the transmitter and the
receiver (AT-HIs) being Gaussian distributed with average powers proportional
to the average transmit and receive signals, respectively [10]. The
Gaussianity is a result of the aggregate contribution of many impairments.
Notably, this model is not only analytically tractable but it is also
experimentally validated [10]. For instance, during the uplink, let
$\frac{p_{k}}{M}=\mathbb{E}\\{|x_{k}|^{2}\\}$ be the transmit power from UE
$k$ transmitting signal $x_{k}$. Then, the AT-HIs are described in terms of
conditional Gaussian distributions as

$\displaystyle\delta_{\mathrm{t},k}$ $\displaystyle\sim{\cal C}{\cal
N}\left(0,\Lambda_{k}\right)\\!\\!$ (3)
$\displaystyle\hbox{\boldmath$\delta$}_{\mathrm{r}}$ $\displaystyle\sim{\cal
C}{\cal N}\left({\mathbf{0}},\bm{\Upsilon}\right)\\!,$ (4)

where $\Lambda_{k}=\kappa_{\mathrm{UE}}\frac{p_{k}}{M}$ and
$\bm{\Upsilon}=\kappa_{\mathrm{BS}}\sum_{i=1}^{K}\frac{p_{i}}{M}\mathrm{diag}(|h_{i,1}|^{2},\ldots,$
$|h_{i,M}|^{2})$ with $\kappa_{\mathrm{UE}}$ and $\kappa_{\mathrm{BS}}$
expressing the severity of the residual impairments at the transmitter and
receiver side, respectively. For simplicity, all UEs are assumed with
identical HIs, i.e., $\kappa_{\mathrm{UE}_{i}}=\kappa_{\mathrm{UE}}~{}\forall
i$. The extension to different distortion at each UE is straightforward. Note
that in the case $\kappa_{\mathrm{UE}}=\kappa_{\mathrm{BS}}=0$ and $m=1$, we
result in the ideal scenario with no HIs.

##  III Uplink Data Transmission with HIs

The received complex baseband signal by the BS is

$\displaystyle{\mathbf{y}}=\sum_{i=1}^{K}{\mathbf{h}}_{i}\left(x_{i}+\delta_{\mathrm{t},i}\right)+\hbox{\boldmath$\delta$}_{\mathrm{r}}+{\mathbf{w}},$
(5)

where $\delta_{\mathrm{t},i}$ and $\hbox{\boldmath$\delta$}_{\mathrm{r}}$ are
the transmit and receive distortions given by (3) and (4), respectively. Also,
${\mathbf{w}}\sim\mathcal{CN}\left({\mathbf{0}},\sigma^{2}{\bm{\mathrm{I}}}_{M}\right)$
is the receiver noise. The signal of UE $k$, detected by the combining vector
${\mathbf{v}}_{k}\in\mathbb{C}^{M\times 1}$, can be expressed as
${\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{y}}$.

###### Lemma 1

The instantaneous uplink SINR of an IRS-assisted MIMO system with AT-HIs and
IRS-HIs is given by

$\displaystyle\gamma_{k}=\frac{\frac{p_{k}}{M}{\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{h}}_{k}{\mathbf{h}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{v}}_{k}}{{\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}\\!\left(\sum_{i\neq
k}\frac{p_{i}}{M}{\mathbf{h}}_{i}{\mathbf{h}}_{i}^{{\scriptscriptstyle\mathsf{H}}}+{\mathbf{C}}_{\\!\delta_{\mathrm{t}}}+{\mathbf{C}}_{\\!\delta_{\mathrm{r}}}+\sigma^{2}{\bm{\mathrm{I}}}_{M}\\!\right)\\!{\mathbf{v}}_{k}},$
(6)

where
${\mathbf{C}}_{\\!\delta_{\mathrm{t}}}=\kappa_{\mathrm{UE}}\sum_{i=1}^{K}\frac{p_{i}}{M}{\mathbf{h}}_{i}{\mathbf{h}}_{i}^{{\scriptscriptstyle\mathsf{H}}}$
and
${\mathbf{C}}_{\\!\delta_{\mathrm{r}}}=\kappa_{\mathrm{BS}}\sum_{i=1}^{K}\frac{p_{i}}{M}{\bm{\mathrm{I}}}_{M}\\!\circ\\!{\mathbf{h}}_{i}{\mathbf{h}}_{i}^{{\scriptscriptstyle\mathsf{H}}}$.

###### Proof:

Given that the AT-HIs are Gaussian distributed and uncorrelated with the
transmit signals, we make use of the worst-case uncorrelated additive noise
theorem in [24] to obtain a lower bound of the mutual information
$\mathcal{I}$ between the input $x_{k}$ and output
${\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{y}}$ for a given
channel realization
${\mathbf{H}}=[{\mathbf{h}}_{1},\ldots,{\mathbf{h}}_{K}]\in\mathbb{C}^{M\times
K}$ as
$\mathbb{E}_{{\mathbf{H}}}\\{\mathcal{I}\left(x_{k};{\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{y}}\right)\\}\geq\mathbb{E}_{{\mathbf{H}}}\\{\log_{2}\left(1+\gamma_{k}\right)\\}$,
where
$\mathbb{E}_{{\mathbf{H}}}\\{\mathcal{I}\left(x_{k};{\mathbf{v}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{y}}\right)\\}$
expresses the ergodic achievable SE with
$\mathbb{E}_{{\mathbf{H}}}\\{\cdot\\}$ denoting the expectation with respect
to ${\mathbf{H}}$, and $\gamma_{k}$ is given by (6) with $\circ$ in
${\mathbf{C}}_{\\!\delta_{\mathrm{r}}}$ denoting the Hadamard product. ∎

###  III-A Problem Formulation

The focal point of this work is the max-min weighted uplink SINR under a
weighted sum-power constraint. The optimization problem is described as

$\displaystyle\\!\\!(\mathcal{P}1)~{}$
$\displaystyle\max_{{\mathbf{V}},{\mathbf{p}},{\bm{\Phi}}}\min_{k}\frac{\gamma_{k}\left({\mathbf{V}},{\mathbf{p}},{\bm{\Phi}}\right)}{\eta_{k}}$
(7)
$\displaystyle~{}\mathrm{s.t.}~{}~{}~{}~{}~{}~{}~{}~{}\frac{1}{M}{\bm{\beta}}^{{\scriptscriptstyle\mathsf{T}}}{\mathbf{p}}\leq
p_{\mathrm{max}},~{}p_{k}>0,~{}\|{\mathbf{v}}_{k}\|=1,\forall k$ (8)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}|\phi_{i}|=1,~{}~{}\forall
i\in\\{1,\ldots,N\\},$ (9)

where ${\mathbf{V}}=[{\mathbf{v}}_{1},\ldots,{\mathbf{v}}_{K}]$ and
${\mathbf{p}}=[p_{1},\ldots,p_{K}]^{{\scriptscriptstyle\mathsf{T}}}$ are the
tuple of receive beamforming matrices and the transmit power vector,
respectively. Also, $p_{\mathrm{max}}$ denotes the given power constraint
while
${\bm{\eta}}=[\eta_{1},\ldots,\eta_{K}]^{{\scriptscriptstyle\mathsf{T}}}$ and
${\bm{\beta}}=[\beta_{1},\ldots,\beta_{K}]^{{\scriptscriptstyle\mathsf{T}}}$
with $\eta_{k}$ and $\beta_{k}$ expressing the priority assigned to UE $k$ and
the weight associated with $p_{k}$.

At optimality, the power coefficients are obtained based on the property that
the weighted SINRs for different UEs are identical, i.e.,
$\frac{\gamma_{1}}{\eta_{1}}=\ldots=\frac{\gamma_{K}}{\eta_{K}}=\tau^{\star}$[19].
As a result, the SINR constraint for UE $k$ is written as

$\displaystyle\gamma_{k}\left({\mathbf{V}},{\mathbf{p}},{\bm{\Phi}}\right)\geq\eta_{k}\tau^{\star}~{}~{}~{}~{}~{}\forall
k$ (10)

##  IV Proposed Design

The problem in $(\mathcal{P}1)$ is non-convex and the coupling among the
optimization variables (the active and passive beamforming at the BS and IRS,
respectively, and the power control) raises difficulties to solve. We tackle
them by following the common alternating optimization in two stages. The
notable difference here is that we take HIs and correlated fading into
account. In particular, first, for any given ${\bm{\Phi}}$, we provide the
optimal linear receiver design in terms of the optimal decoders and the
optimal allocated power. Then, we consider the IRS design.

###### Proposition 1

Given the RBM ${\bm{\Phi}}$ and the power vector ${\mathbf{p}}$, the uplink
SINR of an IRS-assisted MIMO system with AT-HIs and IRS-HIs is maximized by
the HIs-aware LMMSE receiver

$\displaystyle{\mathbf{v}}_{k}^{\star}=\\!{\bigg{(}\\!\displaystyle\sum_{i=1}^{K}\frac{p_{i}}{M}{\mathbf{h}}_{i}{\mathbf{h}}_{i}^{{\scriptscriptstyle\mathsf{H}}}\\!+\\!{\mathbf{C}}_{\\!\delta_{\mathrm{t}}}\\!+\\!{\mathbf{C}}_{\\!\delta_{\mathrm{r}}}\\!+\\!\sigma^{2}{\bm{\mathrm{I}}}_{M}\\!\\!\bigg{)}^{\\!\\!\\!-1}{\mathbf{h}}_{k}},$
(11)

and the optimal SINR is obtained as

$\displaystyle\gamma_{k}=\frac{p_{k}}{M}{\mathbf{h}}_{k}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{\Sigma}}{\mathbf{h}}_{k},$
(12)

where ${\mathbf{\Sigma}}\\!=\\!\bigg{(}\displaystyle\sum_{i\neq
k}^{K}\frac{p_{i}}{M}{\mathbf{h}}_{i}{\mathbf{h}}_{i}^{{\scriptscriptstyle\mathsf{H}}}+{\mathbf{C}}_{\\!\delta_{\mathrm{t}}}+{\mathbf{C}}_{\\!\delta_{\mathrm{r}}}+\sigma^{2}{\bm{\mathrm{I}}}_{M}\\!\\!\bigg{)}^{\\!\\!-1}$.

###### Proof:

The SINR $\gamma_{k}$ in (6) can be written as a generalized Rayleigh quotient
that can be maximized according to [24, Lem. B.10] by
${\mathbf{v}}_{k}^{\star}$ given by (11), where the matrix inversion lemma has
been also applied. A straightforward substitution of (11) into (6) results in
(12). ∎

The optimization of the power allocation relies on a large system analysis of
the max-min weighted SINR.444 The application of MMSE-type receivers to IRS-
assisted mMIMO systems includes prohibitively demanding computations such as
the matrix inversion as $M,N,K$ increase. In addition, the corresponding
calculations should take place in every coherence interval. These reasons
indicate the use of the theory of DEs concerning derivations in the asymptotic
limit $M,N,K\to\infty$ while their ratios are kept fixed. Notably, the DE
results are tight approximations even for a medium system size, e.g., $8\times
8$ [25]. In particular, taking into account the same assumptions concerning
the correlation matrices as in [25, Assump. A1-A2], we obtain the DE SINR
$\bar{\gamma}_{k}$ according to
$\gamma_{k}-\bar{\gamma}_{k}\xrightarrow[M\rightarrow\infty]{\mbox{a.s.}}0$.
The notation $\xrightarrow[M\rightarrow\infty]{\mbox{a.s.}}$ denotes almost
sure convergence as $M\\!\rightarrow\\!\infty$. Hence, the DE weighted SINR is
given by $\bar{\gamma}_{k}/\eta_{k}=\bar{\tau}$.

###### Lemma 2

The DE of the optimal SINR of an IRS-assisted mMIMO system, accounting for
HIs, is given by

$\displaystyle\bar{\gamma}_{k}=p_{k}\frac{\delta_{k}}{1+\kappa_{\mathrm{UE}}p_{k}\delta_{k}},$
(13)

where $\delta_{k}=\frac{1}{M}\tr\left({\mathbf{R}}_{k}{\mathbf{T}}\right)$
with ${\mathbf{T}}=\\!\\!\sum_{i\neq
k}^{K}\\!\frac{p_{i}}{M}\big{(}\\!\frac{\left(1\\!+\\!\kappa_{\mathrm{UE}}\right)}{\left(1+\delta_{i}\right)}{\mathbf{R}}_{i}\\!+\\!\sum_{i=1}^{K}\\!\frac{p_{i}}{M}\kappa_{\mathrm{BS}}{\bm{\mathrm{I}}}_{M}\\!\circ\\!{\mathbf{R}}_{i}\\!+\\!\sigma^{2}{\bm{\mathrm{I}}}_{M}\\!\big{)}^{\\!\\!-1}$.

###### Proof:

We define the matrix ${\mathbf{\Sigma}}_{k}\\!=\\!\big{(}\\!\sum_{i\neq
k}^{K}\\!\frac{p_{i}}{M}\left(1+\kappa_{\mathrm{UE}}\right){\mathbf{h}}_{i}{\mathbf{h}}_{i}^{{\scriptscriptstyle\mathsf{H}}}\\!+\\!{\mathbf{C}}_{\\!\delta_{\mathrm{r}}}\\!+\\!\sigma^{2}{\bm{\mathrm{I}}}_{M}\\!\big{)}^{\\!-\\!1}$.
From (12), the optimal SINR becomes

$\displaystyle\gamma_{k}$
$\displaystyle\xrightarrow[M\rightarrow\infty]{\mbox{a.s.}}p_{k}\frac{\frac{1}{M}\tr\left({\mathbf{\Sigma}}_{k}{\mathbf{R}}_{k}\right)}{1+\frac{1}{M}\kappa_{\mathrm{UE}}p_{k}\tr\left({\mathbf{\Sigma}}_{k}{\mathbf{R}}_{k}\right)}$
(14)
$\displaystyle\xrightarrow[M\rightarrow\infty]{\mbox{a.s.}}p_{k}\frac{\frac{1}{M}\tr\left({\mathbf{R}}_{k}{\mathbf{T}}\right)}{1+\frac{1}{M}\kappa_{\mathrm{UE}}p_{k}\tr\left({\mathbf{R}}_{k}{\mathbf{T}}\right)},$
(15)

where in (14), we have applied the matrix inversion lemma [25, Lem. 1], and
have used [25, Lem. 4]. 555We assume that the diagonal matrix inside
${\mathbf{\Sigma}}_{k}$ is considered deterministic with diagonal elements
given by the limits of the individual diagonal elements [8]. Specifically, we
exploit the uniform convergence $\lim\sup_{M}\max_{1\leq i\leq
M}{\left|\left[{\mathbf{h}}_{i}{\mathbf{h}}_{i}^{{\scriptscriptstyle\mathsf{H}}}\right]_{mm}-\left[{\mathbf{R}}_{i}\right]_{mm}\right|}=0$,
and obtain
$\left\|\frac{1}{M}\mathrm{diag}({\mathbf{h}}_{i}{\mathbf{h}}_{i}^{H})-\frac{1}{M}\tr\left(\text{diag}\left(\right)\left({\mathbf{R}}_{i}\right)\right)\\!\right\|\xrightarrow[M\rightarrow\infty]{\mbox{a.s.}}\\!0$
where $[A]_{mm}$ denotes the $m$th diagonal element of matrix ${\mathbf{A}}$.
The last step makes use of [25, Th. 1] with ${\mathbf{T}}$ given in Lemma 2. ∎

Using the combiner in (11), the SINR constraint in (10) is fulfilled with
equality by choosing the optimal power allocation according to the following
proposition. In other words, the following proposition provides a necessary
and sufficient condition for optimality of $\mathcal{P}1$.

###### Proposition 2

Given the RBM ${\bm{\Phi}}$ and the decoder matrix ${\mathbf{V}}$, the optimal
power vector ${\mathbf{p}}^{\star}$, accounting for HIs, is obtained
geometrically fast as the positive solution to the fixed-point equation given
by

$\displaystyle\\!\frac{1}{\bar{\tau}^{\star}}{\mathbf{p}}^{\\!\star}$
$\displaystyle\\!=\left({\bm{\eta}}\\!\circ\\!{\mathbf{u}}\\!\circ\\!\left(\kappa_{\mathrm{UE}}\hbox{\boldmath$\delta$}\\!\circ\\!{\mathbf{1}}+{\mathbf{1}}{\bm{\beta}}^{{\scriptscriptstyle\mathsf{T}}}\right)\right)\\!{\mathbf{p}}^{\\!\star},$
(16)

where $\bar{\tau}^{\star}$ is the deterministic optimal weighted SINR,
$\hbox{\boldmath$\delta$}=\left[\delta_{1},\ldots,\delta_{K}\right]^{{\scriptscriptstyle\mathsf{T}}}$,
and
${\mathbf{u}}=\left[\frac{1}{G_{11}},\ldots,\frac{1}{G_{KK}}\right]^{{\scriptscriptstyle\mathsf{T}}}$.

###### Proof:

We define
$\hbox{\boldmath$\delta$}=\left[\delta_{1},\ldots,\delta_{K}\right]^{{\scriptscriptstyle\mathsf{T}}}$
and
${\mathbf{u}}=\left[\frac{1}{\delta_{1}},\ldots,\frac{1}{\delta_{K}}\right]^{{\scriptscriptstyle\mathsf{T}}}$.
Now, the weighted SINR can be written as

$\displaystyle\frac{\bar{\tau}_{k}}{\eta_{k}}\\!=\\!\frac{p_{k}}{\left({\bm{\eta}}\\!\circ\\!{\mathbf{u}}\\!\circ\\!\left(\kappa_{\mathrm{UE}}\hbox{\boldmath$\delta$}\\!\circ\\!{\mathbf{p}}+\\!{\mathbf{1}}\right)\right)_{k}},$

where $\bar{\tau}_{k}$ is the deterministic weighted SINR. Finally, taking
advantage of (8) and that, at optimality, the weighted SINR for different UEs
is the same, we result in (16), which converges geometrically fast as follows
from the remark after Theorem 1 in [26]. ∎

Given the decoder matrix ${\mathbf{V}}$ and the optimal power vector
${\mathbf{p}}^{\star}$, the design of the IRS RBM ${\bm{\Phi}}$ is obtained by
means of the optimization problem

$\displaystyle\begin{split}(\mathcal{P}2)~{}~{}~{}~{}~{}~{}~{}\max_{{\bm{\Phi}}}~{}~{}~{}&\bar{\tau}^{*}\\\
\mathrm{s.t}~{}~{}~{}&|\phi_{n}|=1,~{}~{}n=1,\dots,N,\end{split}$ (17)

which is a maximization problem with a unit-modulus constraint regarding
$\phi_{n}$ that can be solved by using projected gradient ascent until
converging to a stationary point as in [4].

Let
${\mathbf{s}}^{i}=[\phi_{1}^{i},\ldots,\phi_{N}^{i}]^{{\scriptscriptstyle\mathsf{T}}}$
be the the induced phases at step $i$ and ${\mathbf{q}}^{i}$ be the adopted
ascent direction at step $i$ with
$[{\mathbf{q}}^{i}]_{n}=\partialderivative{\bar{\tau}^{*}}{\phi_{n}^{*}}$(given
by Lemma (3)), the next iteration point is given by

$\displaystyle\tilde{{\mathbf{s}}}^{i+1}$
$\displaystyle={\mathbf{s}}^{i}+\mu{\mathbf{q}}^{i}$ (18)
$\displaystyle{\mathbf{s}}^{i+1}$
$\displaystyle=\exp\left(j\arg\left(\tilde{{\mathbf{s}}}^{i+1}\right)\right),$
(19)

where $\mu$ is the step size, computed at each iteration by means of the
backtracking line search [27]. The projection problem
$\min_{|s_{n}|=1,n=1,\ldots,N}\|{\mathbf{s}}-\tilde{{\mathbf{s}}}\|^{2}$
provides the solution while satisfying the unit-modulus constraint.

###### Lemma 3

The derivative of $\bar{\tau}^{\star}$ with respect to $\phi_{n}^{*}$ is given
by the fixed-point equation (20) at the top of the next page.

###### Proof:

Refer to Appendix A.∎

After establishing the optimal receiver, power allocation, and RBM, we combine
them by using alternate optimization to find a locally optimal solution. Note
that the non-convexity of $\left(\mathcal{P}1\right)$ cannot guarantee any
global optimality.

Since each subproblem achieves an optimal solution, the objective function of
$\left(\mathcal{P}1\right)$ is non-decreasing over iterations. Moreover, the
optimal value of the objective function is bounded from above due to the power
constraint. Hence, the proposed algorithm converges.

The proposed algorithm is quite advantageous compared to algorithms based on
instantaneous CSI since the power allocation and the phase shifts converge to
deterministic values that depend on statistical CSI. Hence, they can be a
priori calculated and stored while can be updated at every several coherence
intervals due to variation of these channel statistics. On the contrary,
instantaneous CSI algorithms would require frequent optimization taking place
at each coherence interval.

$\displaystyle\partialderivative{\bar{\tau}^{\star}}{\phi_{n}^{*}}\\!$
$\displaystyle=\\!\frac{p_{k}}{\eta_{k}M\left(1+\kappa_{\mathrm{UE}}\delta_{k}\right)^{2}}\alpha\bigg{[}{\mathbf{H}}_{1}^{{\scriptscriptstyle\mathsf{H}}}\Big{(}\beta_{2,k}{\mathbf{T}}{\mathbf{H}}_{1}{\bm{\Phi}}\tilde{{\mathbf{R}}}_{\mathrm{IRS},k}-\kappa_{\mathrm{BS}}\sum_{i=1}^{K}\\!\frac{p_{i}\beta_{2,i}}{M}({\bm{\mathrm{I}}}_{M}\circ{\mathbf{T}}{\mathbf{R}}_{k}{\mathbf{T}}){\mathbf{H}}_{1}{\bm{\Phi}}\tilde{{\mathbf{R}}}_{\mathrm{IRS},i}\\!$
$\displaystyle-\left(1\\!+\\!\kappa_{\mathrm{UE}}\right)\sum_{i\neq
k}^{K}\\!\frac{p_{i}\beta_{2,i}}{M\left(1+\delta_{i}\right)}{\mathbf{T}}{\mathbf{R}}_{k}{\mathbf{T}}{\mathbf{H}}_{1}{\bm{\Phi}}\tilde{{\mathbf{R}}}_{\mathrm{IRS},i}\Big{)}\bigg{]}_{n,n}-\left(1\\!+\\!\kappa_{\mathrm{UE}}\right)\sum_{i\neq
k}^{K}\\!\frac{p_{i}\beta_{2,i}\tr({\mathbf{T}}{\mathbf{R}}_{k}{\mathbf{T}}{\mathbf{R}}_{i})}{M\left(1+\delta_{i}\right)^{2}}\partialderivative{\delta_{i}}{\phi^{*}_{n}}.$
(20)

##  V Numerical Results

We consider a uniform linear array (ULA) and a uniform planar array (UPA) for
the configuration of the BS and IRS, respectively. In particular, we have
$d_{\mathrm{BS}}=d_{\mathrm{IRS}}=0.5\lambda$ while $\theta_{1,n}$,
$\psi_{1,n}$ are uniformly distributed between $0$ to $\pi$ and $0$ to $2\pi$,
respectively. Also, $\theta_{2,n}=\pi-\theta_{1,n}$,
$\psi_{2,n}=\pi+\psi_{1,n}$. Moreover, we employ the 3GPP Urban Micro (UMi)
scenario from TR36.814 for a carrier frequency of $2.5$ GHz and noise level
$-80$ dBm, where the path losses for ${\mathbf{h}}_{2,k}$ and
${\mathbf{H}}_{1}$ are generated based on the NLOS and LOS versions,
respectively [4]. Specifically, therein, the overall path loss for the IRS-
assisted link is $\beta_{k}=\beta_{1,k}\beta_{2,k}$, where

$\displaystyle\beta_{1,k}=C_{1}d_{\mathrm{BS}-\mathrm{IRS}}^{-\nu_{1}},~{}~{}\beta_{2,k}=C_{2}d_{\mathrm{IRS}-\mathrm{UE}_{k}}^{-nu_{2}}$
(21)

with $C_{1}=26$ dB, $C_{2}=28$ dB, $\nu_{1}=2.2$, $\nu_{2}=3.67$. The
variables $d_{\mathrm{BS}-\mathrm{IRS}}$ and
$d_{\mathrm{IRS}-\mathrm{UE}_{k}}$ express the distances between the BS and
IRS, and the IRS and UE $k$, respectively. The penetration losses of the IRS-
assisted links are assumed negligible by deploying the IRS higher than the BS.
For $\beta_{\mathrm{d},k}$, we assume the same parameters as for
$\beta_{2,k}$, but we also consider an additional penetration loss equal to
$15~{}\mathrm{dB}$. We use $5$ dBi antennas at the BS and IRS, and
${\mathbf{R}}_{\mathrm{BS},k},{\mathbf{R}}_{\mathrm{IRS},k}$ are generated as
in [28, 4], respectively. The size of each IRS-element dimension is
$\lambda/4$. The “solid” and “dashed” lines correspond to
$p_{\mathrm{max}}=0~{}\mathrm{dB}$ and $p_{\mathrm{max}}=20~{}\mathrm{dB}$,
respectively, while different line symbols correspond to different values of
impairments, given by
$\kappa_{\mathrm{BS}}=\kappa_{\mathrm{UE}}=\\{0,0.05^{2},0.1^{2}\\}$,
respectively. The optimization takes place by choosing arbitrary values for
the phase shifts and the power as in [19]. For simplicity, we assume
$\alpha=1$ and that all data streams have the same priority
$({\bm{\eta}}={\mathbf{1}})$ and power weight associated
(${\bm{\beta}}=\frac{1}{K}{\mathbf{1}}$).

Fig. 1.(a) depicts the minimum uplink user rate
$\log_{2}\left(1+\bar{\tau}\right)$ with respect to the number of IRS elements
for different IRS-HIs and correlation conditions with no AT-HIs. The line
describing the ideal case (no HIs) appears for comparison. In the case of
uniform phase noise $m=0$, the rate takes the lowest value since
${\mathbf{R}}_{k}$ does not depend on the phase shifts. However, when the Von
Mises PDF is assumed (“solid” lines), the IRS can be optimized and the rate
increases. In particular, based on the concentration parameter
$\kappa_{\tilde{\theta}}$, we observe that its decrease results in the
decrease of the rate. In the special case $\kappa_{\tilde{\theta}}=0$ (“star”
symbols), the line coincides with the line describing the uniform
distribution. Moreover, if no IRS correlation is assumed (“$\times$” symbols),
the rate is lower since the rate cannot be maximized due to IRS exploitation
(see Rem. 1). Also, we have added a “dotted” line corresponding to
$\kappa_{\tilde{\theta}}=2$ that describes a “naive” scheme, where the impact
of no HIs is taken into account during the optimization design. Its lower SINR
reveals the robustness of the proposed design.

In Fig. 1.(b), we illustrate the minimum uplink user rate with respect to the
number of IRS elements for different SNR values and varying AT-HIs (no IRS-
HIs). Despite that the rate generally increases with SNR, we now observe its
increase with $N$. In the case of perfect hardware, the rate appears no
ceiling as $N$ increases, but it saturates in practice, where AT-HIs are met.
Moreover, the saturation appears earlier in the case of higher SNR
($20~{}\mathrm{dB}$). Although, the more severe AT-HIs result in higher
degradation, we depict the impact of the BS additive distortion
$\kappa_{\mathrm{BS}}$ as $N$ increases by the “dot” lines in the circle. We
notice that the curves converge to the same value when $N\to\infty$, which
means that the impact of $\kappa_{\mathrm{BS}}$ is negligible at large $N$,
i.e., the larger the IRS, the more beneficial the communication with mMIMO
despite the use of low-quality transceiver hardware. Furthermore, the “star”
lines correspond to Monte Carlo simulations verifying the DE results.

Fig. 1.(c) shows the minimum user rate versus the number of BS antennas for
varying AT-HIs (no IRS-HIs). Especially, we observe that the variations of $M$
and $N$ present similar behavior. In the case of perfect hardware, the
achievable rate increases unboundedly as $M\to\infty$. However, when the AT-
HIs are taken into account, we observe finite limits. Notably, the lower
quality (more severe AT-HIs) results in larger degradation. Moreover, apart
from that, a higher SNR leads to a higher rate, we observe that the
convergence to saturation between different SNRs is different since AT-HIs are
power-dependent. Thus, at $20~{}\mathrm{dB}$, the rate saturates faster, i.e.,
the majority of the multi-antenna gain takes place at low $M$, but still, a
large $M$ contributes to larger multiplexing and inter-user interference
mitigation. Hence, an IRS-assisted mMIMO works better at high SNR values. In
addition, we observe that, at $0~{}\mathrm{dB}$, the convergence requires more
antennas.

$(a)$

$(b)$

$(c)$

Figure 1: Asymptotic weighted max-min user rate of an IRS-assisted mMIMO
system for varying $(a)$ IRS-HIs and correlation (correlated/uncorrelated
Rayleigh fading channels) conditions versus the number of IRS elements $N$
($M=80$, $K=5$); AT-HIs in terms of $\kappa_{\mathrm{BS}}$,
$\kappa_{\mathrm{UE}}$ and transmit power budget $p_{\mathrm{max}}$ versus
$(b)$ the number of IRS elements $N$ ($M=80$, $K=5$); and $(c)$ the number of
BS antennas $M$ ($N=80$, $K=5$).

##  VI Conclusion

This paper provided a thorough investigation of the impact of AT-HIs and IRS-
HIs on an IRS-assisted mMIMO system. Specifically, we obtained the optimal
asymptotic max-min weighted uplink SINR by optimizing the transmit power, the
HIs-aware receiver, and the RBM. Remarkably, this asymptotic SINR, being
dependent only on the large-scale statistics, allows optimizing the transmit
power and the RBM only on every several coherence when these statistics
change. We verified the tightness of the analytical expression by simulations
even for practical system dimensions. Moreover, we shed light on the impact of
HIs on the asymptotic max-min weighted SINR by varying the hardware quality at
both the IRS and the transceiver.

##  Appendix A Proof of Lemma 3

After adding and subtracting $1$ in the numerator of (13), the partial
derivative of $\bar{\tau}^{\star}$ is written as

$\displaystyle\partialderivative{\bar{\tau}^{\star}}{\phi_{n}^{*}}=\frac{p_{k}}{\eta_{k}}\frac{1}{\left(1+\kappa_{\mathrm{UE}}\delta_{k}\right)^{2}}\partialderivative{\delta_{k}}{\phi_{n}^{*}}.$
(22)

We have

$\displaystyle\partialderivative{\delta_{k}}{\phi^{*}_{n}}=\frac{1}{M}\tr\left(\partialderivative{{\mathbf{R}}_{k}}{\phi^{*}_{n}}{\mathbf{T}}+{\mathbf{R}}_{k}\partialderivative{{\mathbf{T}}}{\phi^{*}_{n}}\right)$
(23)
$\displaystyle=\frac{1}{M}\tr\bigg(\partialderivative{{\mathbf{R}}_{k}}{\phi^{*}_{n}}{\mathbf{T}}\bigg{missing})-\frac{1}{M}\tr\bigg({\mathbf{R}}_{k}{\mathbf{T}}\partialderivative{{\mathbf{T}}^{-1}}{\phi^{*}_{n}}{\mathbf{T}}\bigg{missing})$
(24)
$\displaystyle\\!=\\!\frac{1}{M}\tr\\!\bigg{(}\\!\partialderivative{{\mathbf{R}}_{k}}{\phi^{*}_{n}}{\mathbf{T}}\bigg{)}\\!-\\!\frac{1}{M}\kappa_{\mathrm{BS}}\tr\\!\bigg{(}\\!{\mathbf{R}}_{k}{\mathbf{T}}\bigg{(}\\!\sum_{i=1}^{K}\\!\frac{p_{i}}{M}\\!\partialderivative{\big{(}{\bm{\mathrm{I}}}_{M}\\!\circ\\!{\mathbf{R}}_{i}\big{)}}{\phi^{*}_{n}}\\!{\mathbf{T}}\\!\bigg{)}$
$\displaystyle-\\!\frac{1}{M}\\!\tr\\!\bigg{(}\\!{\mathbf{R}}_{k}{\mathbf{T}}\bigg{(}\\!\sum_{i\neq
k}^{K}\\!\frac{p_{i}\left(1\\!+\\!\kappa_{\mathrm{UE}}\right)}{M\left(1+\delta_{i}\right)}\big{(}\\!\partialderivative{{\mathbf{R}}_{i}}{\phi^{*}_{n}}\\!-\\!\frac{\partialderivative{\delta_{i}}{\phi^{*}_{n}}}{\left(1\\!+\\!\delta_{i}\right)}{\mathbf{R}}_{i}\\!\big{)}\\!{\mathbf{T}}\\!\bigg{)}\\!,$
(25)

where, in (24), we used the derivative of the inverse matrix ${\mathbf{T}}$.

Use of [13, Lem. 1] in (25) and substitution into (22) concludes the proof
after several algebraic manipulations.

## References

  * [1] E. Basar _et al._ , “Wireless communications through reconfigurable intelligent surfaces,” _IEEE Access_ , vol. 7, pp. 116 753–116 773, 2019\.
  * [2] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 11, pp. 5394–5409, 2019.
  * [3] A. M. Elbir _et al._ , “Deep channel learning for large intelligent surfaces aided mm-wave massive MIMO systems,” _IEEE Wireless Commun. Lett._ , 2020.
  * [4] A. Kammoun _et al._ , “Asymptotic max-min SINR analysis of reconfigurable intelligent surface assisted MISO systems,” _IEEE Trans. Wireless Commun._ , 2020.
  * [5] A. Papazafeiropoulos _et al._ , “Coverage probability of distributed IRS systems under spatially correlated channels,” _IEEE Wireless Commun. Lett._ , pp. 1–1, 2021.
  * [6] T. Van Chien _et al._ , “Outage probability analysis of IRS-assisted systems under spatially correlated channels,” _IEEE Wireless Commun. Lett._ , 2021.
  * [7] E. Björnson _et al._ , “Massive MIMO systems with non-ideal hardware: Energy efficiency, estimation, and capacity limits,” _IEEE Trans. Inf. Theory_ , vol. 60, no. 11, pp. 7112–7139, 2014.
  * [8] A. Papazafeiropoulos, B. Clerckx, and T. Ratnarajah, “Rate-splitting to mitigate residual transceiver hardware impairments in massive MIMO systems,” _IEEE Trans. Veh. Tech._ , vol. 66, no. 9, pp. 8196–8211, 2017\.
  * [9] A. K. Papazafeiropoulos _et al._ , “Ergodic capacity analysis of AF DH MIMO relay systems with residual transceiver hardware impairments: Conventional and large system limits,” _IEEE Trans. Veh. Tech._ , vol. 66, no. 8, pp. 7010–7025, 2017.
  * [10] C. Studer, M. Wenk, and A. Burg, “MIMO transmission with residual transmit-RF impairments,” in _ITG/IEEE Work. Smart Ant. (WSA)_. IEEE, 2010, pp. 189–196.
  * [11] M.-A. Badiu and J. P. Coon, “Communication through a large reflecting surface with phase errors,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 2, pp. 184–188, 2019.
  * [12] D. Li, “Ergodic capacity of intelligent reflecting surface-assisted communication systems with phase errors,” _IEEE Commun. Lett._ , vol. 24, no. 8, pp. 1646–1650, 2020.
  * [13] A. Papazafeiropoulos _et al._ , “Intelligent reflecting surface-assisted MU-MISO systems with imperfect hardware: Channel estimation, beamforming design,” _arXiv preprint arXiv:2102.05333_ , 2021.
  * [14] Z. Xing and R. Wang, “Achievable rate analyses and phase shift optimizations on intelligent reflecting surface with hardware impairments,” _arXiv preprint arXiv:2005.14411_ , 2020.
  * [15] Y. Liu _et al._ , “Beamforming designs and performance evaluations for intelligent reflecting surface enhanced wireless communication system with hardware impairments,” _arXiv preprint arXiv:2006.00664_ , 2020.
  * [16] S. Zhou _et al._ , “Spectral and energy efficiency of IRS-assisted MISO communication with hardware impairments,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 9, pp. 1366–1369, 2020.
  * [17] H. Shen _et al._ , “Beamforming optimization for IRS-aided communications with transceiver hardware impairments,” _IEEE Trans. Commun._ , vol. 69, no. 2, pp. 1214–1227, 2021.
  * [18] G. Zhou _et al._ , “Secure wireless communication in RIS-aided MISO systems with hardware impairments.”
  * [19] D. W. Cai, T. Q. Quek, and C. W. Tan, “A unified analysis of max-min weighted SINR for MIMO downlink system,” _IEEE Trans. Signal Process._ , vol. 59, no. 8, pp. 3850–3862, 2011.
  * [20] D. Neumann, M. Joham, and W. Utschick, “Covariance matrix estimation in massive MIMO,” _IEEE Signal Process. Lett._ , vol. 25, no. 6, pp. 863–867, 2018.
  * [21] F. Bohagen, P. Orten, and G. E. Oien, “Design of optimal high-rank line-of-sight MIMO channels,” _IEEE Trans. Wireless Commun._ , vol. 6, no. 4, pp. 1420–1425, 2007.
  * [22] G. Gradoni and M. Di Renzo, “End-to-end mutual coupling aware communication model for reconfigurable intelligent surfaces: An electromagnetic-compliant approach based on mutual impedances,” _IEEE Wireless Commun. Lett._ , 2021\.
  * [23] Y. Zhang _et al._ , “Performance analysis of RIS-aided systems with practical phase shift and amplitude response,” _IEEE Trans. Veh. Tech._ , 2021.
  * [24] E. Björnson _et al._ , “Massive MIMO networks: Spectral, energy, and hardware efficiency,” _Foundations and Trends® in Signal Processing_ , vol. 11, no. 3-4, pp. 154–655, 2017.
  * [25] J. Hoydis, S. ten Brink, and M. Debbah, “Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?” _IEEE J. Select. Areas Commun._ , vol. 31, no. 2, pp. 160–171, February 2013.
  * [26] U. Krause, “Concave perron-frobenius theory and applications,” vol. 47, no. 3, pp. 1457–1466, 2001.
  * [27] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [28] E. Björnson and L. Sanguinetti, “Rayleigh fading modeling and channel hardening for reconfigurable intelligent surfaces,” _IEEE Wireless Commun. Lett._ , 2020.

solutionfile