# Prediction of crystal structures and motifs in the Fe-Mg-O system at Earth’s
core pressures

Renhai Wang  Department of Physics, University of Science and Technology of
China, Hefei 230026, China  Department of Physics, Iowa State University,
Ames, Iowa 50011, USA  Yang Sun<EMAIL_ADDRESS>Department of Applied
Physics and Applied Mathematics, Columbia University, New York, NY, 10027, USA
Renata M. Wentzcovitch<EMAIL_ADDRESS>Department of Applied Physics and
Applied Mathematics, Columbia University, New York, NY, 10027, USA  Department
of Earth and Environmental Sciences, Columbia University, New York, NY, 10027,
USA  Lamont–Doherty Earth Observatory, Columbia University, Palisades, NY,
10964, USA  Feng Zheng  Department of Physics, Xiamen University, Xiamen
361005, China  Yimei Fang  Department of Physics, Xiamen University, Xiamen
361005, China  Shunqing Wu  Department of Physics, Xiamen University, Xiamen
361005, China  Zijing Lin  Department of Physics, University of Science and
Technology of China, Hefei 230026, China  Cai-Zhuang Wang  Department of
Physics, Iowa State University, Ames, Iowa 50011, USA  Kai-Ming Ho  Department
of Physics, Iowa State University, Ames, Iowa 50011, USA

(Feb. 4, 2021)

###### Abstract

Fe, Mg, and O are among the most abundant elements in terrestrial planets.
While the behavior of the Fe-O, Mg-O, and Fe-Mg binary systems under pressure
have been investigated, there are still very few studies of the Fe-Mg-O
ternary system at relevant Earth’s core and super-Earth’s mantle pressures.
Here, we use the adaptive genetic algorithm (AGA) to study ternary FexMgyOz
phases in a wide range of stoichiometries at 200 GPa and 350 GPa. We
discovered three dynamically stable phases with stoichiometries FeMg2O4,
Fe2MgO4, and FeMg3O4 with lower enthalpy than any known combination of Fe-Mg-O
high-pressure compounds at 350 GPa. With the discovery of these phases, we
construct the Fe-Mg-O ternary convex hull. We further clarify the composition-
and pressure-dependence of structural motifs with the analysis of the AGA-
found stable and metastable structures. Analysis of binary and ternary stable
phases suggest that O, Mg, or both could stabilize a BCC iron alloy at inner
core pressures.

##  I Introduction

O, Fe, Si, Mg, Al, and Ca (CMAS+F) are the most abundant elements in
terrestrial planets [1]. Among these planets, Earth provides essential general
information, yet it is incompletely deciphered. All CMAS+F elements are
lithophile (rock-loving) elements and are present in the Earth’s rocky mantle
and crust. Fe is the predominant element in the core and is a siderophile
(metal-loving) element as well. Based on current knowledge, this
classification is believed to be valid in the pressure and temperature (PT)
range achieved in Earth’s interior. Seismology and high-pressure data on iron
shows that the Earth’s core is $\sim$ 5-10 wt% [2] less dense than iron at
expected conditions, i.e., $\sim$ 136-364 GPa and $\sim$ 4,000-6,500 K [3, 4,
5, 6, 7, 8]. This indicates the presence of lighter elements in the core
partitioned differently between its solid and liquid regions [9]. Extensive
research has been carried out experimentally and computationally to shed light
on the light elements’ nature. Despite much progress, there is still a great
deal of uncertainty [2, 10, 11, 12]. With every experimental or computational
development, this question is revisited from a different angle. In particular,
the development of materials discovery methods [13, 14, 15, 16, 17, 18, 19,
20] has propelled the exploration of novel chemistries under pressure, which
has fueled the debate on the possible nature of light elements in the core.
Among the CMAS elements, the most likely light element candidates in the core
are Si, and O. O is considered a required element today, a view that has
evolved within the last decade [21, 22, 23, 24]. The volatile elements S, C,
and H are also regarded as likely candidates, but the abundances of C and H on
Earth are still largely unconstrained. Mg, Al, and Ca have been branded
lithophile elements up to core pressures. But the recent computational
discovery of Mg-Fe compounds up to inner core pressures [25] suggests the
possibility of Mg turning siderophile and its presence in the core. The
formation of Fe-O [16] and Fe-Si [15] compounds with variable stoichiometry
have been investigated using materials discovery methods. The theoretical
prediction of pyrite-type FeO2 [16] and its experimental confirmation [26, 27]
has been one of the greatest successes of this approach, which rarely explores
the possibility of ternary compounds [17]. Given the present understanding
that O is a required element in the outer core and should also exist in the
inner core, we explore the possible formation of Fe-Mg-O compounds at typical
core pressures of $\sim$ 350 GPa. The investigation of solid compounds
provides the most critical information. Light elements must be present in both
solid and liquid phases but more abundantly in the liquid phase. It is
energetically more costly to accommodate these elements in the solid phase, a
geochemical definition of incompatible elements. Therefore, the discovery of
thermodynamically stable Fe-Mg-O solids is essential to investigate Mg’s
presence in the core.

Besides being essential for addressing Mg’s presence in the Earth’s core, the
present study of Fe-Mg-O solids has significant ramifications for the mantle
of terrestrial planets larger than Earth whose interiors can reach much higher
pressures and temperatures. The B1-type isomorphous alloy (Mg1-xFex)O,
ferropericlase (x $<$ 0.5) or magnesiumwüstite (x $>$ 0.5) is the second most
abundant phase of the Earth’s lower mantle. Ferropericlase is the
thermodynamically stable form of the Fe-Mg-O ternary compound up to $\sim$ 135
GPa and $\sim$ 4,000 K, i.e., CMB conditions in the Earth. The higher
pressures and temperatures expected in Super-Earths, e.g., $\sim$ 4 TPa and
$\sim$ 9,000 K at the CMB in a 20 M⨁ terrestrial planet [28], raises the
possibility of other compounds and alloy structures with other compositions.
The existence of other stable ternary phases under pressure may induce
decomposition and recombination reactions between Fe-Mg-O, Fe-O, Mg-O, and Fe-
Mg compounds under pressure, similar to what has been observed in the Si-Mg-O
system [17, 18]. Therefore, the current research can also provide the first
glimpses on the essential (Mg1-xFex)O alloy behavior in these planets.

This paper is organized as follows. In the next section, we describe the
computational method used in this work. Sec. III presents the crystal
structure search results, the ternary convex hull, and the predominant
structural motifs in low enthalpy structures. Section IV discusses some
potential geophysical implications of these results. Conclusions are present
in Section V.

##  II Methods

Crystal structures of Fe-Mg-O at high pressure were investigated using the
adaptive genetic algorithm (AGA)[17, 29]. This method integrates auxiliary
interatomic potentials and ab initio calculations adaptively. The auxiliary
interatomic potentials accelerate crystal structure searches in the genetic
algorithm (GA) loop. At the same time, ab initio calculations are used to
adapt the potentials after several GA generations to ensure accuracy. The
structure searches were only constrained by the chemical composition, without
any assumption on the Bravais lattice type, symmetry, atomic basis, or unit
cell dimensions. In our AGA searches, the enthalpy was used as the selection
criteria for optimizing the candidate pool. The candidate structure pool size
in GA search is 128. At each GA generation, 32 new structures are generated
from the parent structure pool via a mating procedure described in [30]. The
structures in the pool were updated by keeping the lowest-energy 128
structures. The structure search with a given auxiliary interatomic potential
sustained 600 consecutive GA generations. Then, 16 structures from the GA
search were randomly selected for ab initio calculations to re-adjust the
interatomic potential parameters for the next round of the GA search. This
sequence of steps was repeated 40 times. In the AGA search in the Fe-Mg-O
system, interatomic potentials based on the embedded-atom method (EAM) [31]
were chosen as the auxiliary classical potential. In EAM, the total energy of
an N-atom system is described by

$E_{total}=\frac{1}{2}\sum_{i,j\left(i\neq
j\right)}^{N}{\varphi\left(r_{ij}\right)+\sum_{i}{F_{i}\left(n_{i}\right)}},$
(1)

where $\varphi\left(r_{ij}\right)$ is the pair term for atoms $i$ and $j$ at a
distance $r_{ij}$. $F_{i}\left(n_{i}\right)$ is the embedded term with
electron density term $n_{i}=\sum_{j\neq i}{\rho_{j}\left(r_{ij}\right)}$ at
the site occupied by atom $i$. The fitting parameters in the EAM formula are
chosen as follows: The parameters for Fe-Fe and Mg-Mg interactions were taken
from the literature [32]. Other pair interactions (O-O, Fe-Mg, Fe-O and Mg-O)
were modeled with the Morse function,

$\varphi\left(r_{ij}\right)=D\left[e^{-2\alpha\left(r_{ij}-r_{0}\right)}-2e^{-\alpha\left(r_{ij}-r_{0}\right)}\right],$
(2)

where $D$, $\alpha$, $r_{0}$ are fitting parameters. The density function for
O atoms are modeled by an exponentially decaying function,

$\rho\left(r_{ij}\right)=\alpha e\ x\
p{\left[-\beta\left(r_{ij}-r_{0}\right)\right]},$ (3)

where $\alpha$ and $\beta$ are fitting parameters. The form proposed by
Benerjea and Smith [33] was used as the embedding function with fitting
parameters $F_{0}$, $\gamma$ as

$F\left(n\right)=F_{0}\left[1-\gamma\ In\ n\right]n^{\gamma}.$ (4)

For Fe and Mg, the parameters of the density function and embedding function
were taken from ref.[32] as well. In the AGA scheme[17], the potential
parameters were adjusted adaptively by fitting to the ab initio energies,
forces, and stresses of selected structures. The fitting process was performed
using the force-matching method with a stochastic simulated annealing
algorithm implemented in the POTFIT code[34, 35].

Ab initio calculations were carried out using the projector augmented wave
(PAW) method[36] within density functional theory (DFT) as implemented in the
VASP code [37, 38]. The exchange and correlation energy are treated without
the spin-polarized generalized gradient approximation (GGA) and parameterized
by the Perdew-Burke-Ernzerhof formula (PBE) [39]. A plane-wave basis was used
with a kinetic energy cutoff of 520 eV, and the convergence criterion for the
total energy was set to $10^{-5}$ eV. Monkhorst-Pack’s sampling scheme [40]
was adopted for Brillouin zone sampling with a k-point grid of $2\pi\times
0.033\text{\AA}^{-1}$, and the unit cell lattice vectors (both the unit cell
shape and size) are fully relaxed under fixed pressure (200 GPa and 350 GPa)
together with the atomic coordinates until the force on each atom is less than
0.01 eV/Å. The phonon dispersions were computed with density functional
perturbation theory (DFPT) implemented in the VASP code and the Phonopy
software [41]. The formation enthalpy (Hf) of compound FexMgyOz was calculated
as

$H_{f}=\frac{H\left({\rm Fe}_{x}{\rm
Mg}_{y}O_{z}\right)-xH\left(Fe\right)-yH\left(Mg\right)-zH\left(O\right)}{x+y+z},$
(5)

where $H\left({\rm Fe}_{x}{\rm Mg}_{y}O_{z}\right)$ is the total enthalpy of
the FexMgyOz alloy. $H\left(Fe\right)$, $H\left(Mg\right)$ and
$H\left(O\right)$ are the enthalpy of the ground state of Fe, Mg, and O at
corresponding pressures, i.e., hcp-Fe, bcc-Mg, and $\zeta$-O2, respectively.

##  III Results and Discussion

###  III.1 AGA search for ternary Fe-Mg-O phases

In Fig. 1(a) we present the AGA results of Fe-Mg-O at 350 GPa. For the sake of
simplicity, all chemical formulae are expressed as Fe/Mg/O reduced ratios. For
example, 123 represents the compound with FeMg2O3. During the structural
search, we select a range of different stoichiometries surrounding 111
composition (i.e., 211, 121, 112, 311, 131, 113, 411, 141, 114, 221, 122, 212,
331, 133, 313, 441, 144, 414, 332, 233, 323, 321, 312, 123, 132, 231, 213,
421, 412, 124, 142, 241, 214, 431, 413, 134, 143, 341 and 314) with 2 or 4
formula units to perform the AGA search (up to 32 atoms per primitive cell).
After the AGA search, we use the following method to determine the stability
of compounds. For a ternary compound $\text{A}_{x}\text{B}_{y}\text{C}_{z}$,
we select three existing compounds $\text{A}_{x1}\text{B}_{y1}\text{C}_{z1}$,
$\text{A}_{x2}\text{B}_{y2}\text{C}_{z2}$, and
$\text{A}_{x3}\text{B}_{y3}\text{C}_{z3}$ on the diagram; these can also be
elementary or binary end-members. If $\text{A}_{x}\text{B}_{y}\text{C}_{z}$
can be written as $a\ \times\ \text{A}_{x1}\text{B}_{y1}\text{C}_{z1}+b\
\times\ \text{A}_{x2}\text{B}_{y2}\text{C}_{z2}+c\ \times\
\text{A}_{x3}\text{B}_{y3}\text{C}_{z3}$ with $a\geq 0$, $b\geq 0$ and $c\geq
0$, we compute its relative enthalpy $\Delta H=H(A_{x}B_{y}C_{z})–a\times
H(A_{x1}B_{y1}C_{z1})–b\times H(A_{x2}B_{y2}C_{z2})–c\times
H(A_{x3}B_{y3}C_{z3})$. If $\Delta H\leq 0$ for all scanned combinations of
existing phases on the diagram, $A_{x}B_{y}C_{z}$ is determined as an
energetic ground state. The energetic ground states form the convex hull as
shown in Fig. 1(a). $H_{d}$ is introduced as the enthalpy above the convex
hull to represent the relative stability on the phase diagram. By definition,
all the ground-state phases have $H_{d}=0$.

The AGA search found three new ternary ground state compounds at 350 GPa:
FeMg2O4, Fe2MgO4, and FeMg3O4. They define the current Fe-Mg-O ternary phase
diagram. The stability of these phases from 200 GPa and 350 GPa is shown in
Fig. 1(b). This stability pressure range is computed by considering the
relative stability of these phases against decomposition into all end-members
(see Supplementary Fig. 1 for the stability range of all ground-state phases).
We will discuss the construction of the phase diagram in the next section. We
also identify low-enthalpy metastable structures such as FeMgO3 with enthalpy
very close to the convex hull ($H_{d}=18$ meV/atom). Here we first analyze
these new ground states and low-enthalpy structures.

Figure 1:  (a) AGA search results for the Fe-Mg-O system at 350 GPa. The dots
represent searched ternary compositions. The color bar corresponds to the
relative enthalpy above the convex hull. The ground-states ($H_{d}=0$) on the
convex hull are connected. The new phases are indicated by the text. 124, 134,
214 and 113 represents FeMg2O4, FeMg3O4, Fe2MgO4 and FeMgO3, respectively. (b)
Stability range of discovered ternary ground states at $T=0K$. The gray bars
indicates the pressure range of decomposition.

Figure 2 shows the atomic structure, phonon dispersion, and electronic density
of states for tetragonal FeMg2O4 with space group I-42d. Fe and Mg atoms are
coordinated with eight oxygen atoms to form a similar MO8 (M for metal)
polyhedra. Unlike a typical cubic polyhedron, this MO8 consists only of
triangular faces. These triangular faces form pentagonal caps and are similar
to the Frank-Kasper polyhedra [42]. The Fe- and Mg-centered MO8 polyhedra pack
in various edge- and face-sharing arrangements. This structure is the same as
the I-42d-type Mg2SiO4 found previously [17, 43]. As shown in Fig. 1(b), this
phase becomes the ground state at 349 GPa. Below this pressure, it decomposes
into MgO and FeO2. As shown in Fig. 2(b), there are no imaginary phonon mode
frequencies, which confirms this phase is dynamically stable. The electronic
density of state in Fig. 2(c) shows a metallic state with somewhat localized
Fe and O states near the Fermi level.

Figure 2:  (a) Atomic structure of I-42d FeMg2O4 and Fe and Mg coordination
polyhedra. Blue is Fe, green is Mg and red is O. Red dashed lines indicate a
pentagonal cap; (b) phonon dispersion; (c) electronic density of states.

Figure 3 shows the atomic structure, phonon dispersion, and electronic density
of states for monoclinic Fe2MgO4 with space group Cc. By inspecting the
structure, we identify the same type of MgO8 polyhedra found in I-42d FeMg2O4.
However, the Fe-O polyhedra are more complicated. It contains two polyhedral
types, one six-fold and one seven-fold coordinated as shown in Fig. 3(a). The
FeO6 is a highly distorted octahedron. The FeO7 also shows the pentagonal cap
similar to the MO8 polyhedron found in I-42d FeMg2O4. As shown in Fig. 1(b),
this phase’s stability pressure range against the decomposition into MgO and
Fe2O3 starts at 325GPa. Phonon calculations confirm its dynamic stability. The
electronic density of states also indicates this phase is metallic in Fig.
3(c).

Figure 3:  (a) Atomic structure Cc Fe2MgO4 and Fe and Mg coordination
polyhedra. Blue is Fe, green is Mg and red is O. Red dashed lines indicate a
pentagonal cap; (b) phonon dispersion; (c) electronic density of states.

The FeMg3O4 Cmmm structure in Fig. 4 shows Fe-O and Mg-O octahedral building
blocks. Visually it is similar to ferropericlase (Fe1-xMgx)O, which has a
NaCl-type (B1) structure. However, unlike the cubic structure and random Fe/Mg
cation site occupancies of ferropericlase, FeMg3O4 is an orthorhombic
structure with ordered Fe/Mg site occupancies. The octahedra in the FeMg3O4
structure are Jahn-Teller distorted because of the orthorhombic symmetry. This
phase becomes stable against decomposition into FeO and MgO at $\sim$ 228 GPa.
Phonon calculations in Fig. 4(b) also confirm its dynamical stability. Unlike
the metallic Fe2MgO4 and FeMg2O4 phases, FeMg3O4 Cmmm is a semiconductor.

Figure 4:  (a) Atomic structure of Cmmm FeMg3O4 and Fe and Mg coordination
polyhedral. Blue is Fe, green is Mg, and red is O; (b) phonon dispersion; (c)
electronic density of states.

Besides the ternary ground states, we also analyze a metastable FeMgO3 Immm
structure with $H_{d}=18$ meV/atom. This enthalpy difference is so small that
the compound may become stable at high temperatures. Structural analysis of
FeMgO3 Immm in Fig. 5 shows an interesting combination of various polyhedra.
Iron shows three different oxygen coordination pulyhedra FeO6, FeO7 and FeO8.
The FeO6 is an octahedron. The FeO7 is a trigonal prism with an extra
rectangular face capping neighbor. The FeO8 is a cube. MgO shows one
coordination polyhedron type, MgO8, not a cube but a triangular prism with two
rectangular face capping oxygens. Such a combination of octahedra and prisms
is similar to the building blocks in the complex Fe2O3 polytypes [44]. The
appearance of cubic polyhedron is consistent with the observation of the CsCl-
type (B2) structure of FeO at Earth’s core conditions [11]. This FeMgO3 Immm
structure shows an intermediate packing between NaCl-type FeO/MgO phase to the
CsCl-type FeO/MgO phases. This phase is dynamically stable and metallic, as
shown by the phonon dispersion and the electronic density of states in Fig. 5.

Figure 5:  (a) Atomic structure of metastable Immm FeMgO3 and Fe and Mg
coordination polyhedra. Blue is Fe, green is Mg and red is O. Red dashed lines
indicate the trigonal prism; (b) phonon dispersion; (c) electronic density of
states.

###  III.2 Construction of the Fe-Mg-O ternary convex hull

In this section, we discuss the construction of the ternary phase diagram and
convex hull. In a binary system, the compositional space is one-dimensional so
that the convex hull is a curve connecting the formation enthalpies of ground-
state phases. In a ternary system, the compositional space is two-dimensional,
and the convex hull consists of surface segments connecting the formation
enthalpies of three stable phases, as shown in Fig. 1. In a discrete
compositional space, these surface segments are triangles. Any new structure
having formation enthalpy below this convex hull surface will be a new ground
state. The convex hull surface needs to be reconstructed after the discovery
of any new stable phase.

For binary references at 350 GPa, the ground-state phases of Fe-O, Fe-Mg, and
Mg-O have been investigated in Refs. [16], [25] and [45], respectively (see
Supplementary Fig. 2 for their crystal structures). Because these crystal
structure searches have already covered the current study’s pressure range, we
do not perform new AGA searches for the binary phases. Still, we re-calculate
the energetics of the previously found crystal structures. Ab initio
calculations confirm these reported relative phase stabilities of Fe-O [16],
Fe-Mg [25], and Mg-O [45]. By exploring an experimental database [46], we find
two previously reported MgFe2O4 stoichiometric compounds [47,48]. However, our
calculations indicate these phases are metastable at Earth’s core pressures.
Based on these ground-state binary phases, we established the Fe-Mg-O ternary
system’s convex-hull shown in Fig. 6. At 200 GPa, all the AGA searched ternary
compounds have relatively higher enthalpy than the elementary or binary
ground-state references. Therefore no Fe-Mg-O stoichiometric phase can be a
ternary ground state at 200 GPa. At 350 GPa three ternary phases become stable
ground states. Detailed energetics and crystallographic information on these
ground-state phases is given in Supplementary Tables S1, Table S2 and Table
S3.

Figure 6:  Ground-state Fe-Mg-O phases and convex hull at (a) 200 GPa and (b)
350 GPa. The blue text refers to stable ternary phases.

###  III.3 Analysis of structural motifs under pressure

Since Fig. 6 suggests a strong effect of pressure on these phases’ stability,
we now investigate how the structural motifs change under pressure. Because
the current calculation does not include temperature effects on phase
relations and at finite temperatures, the ground states may differ. We now
focus on the ground-state phases and metastable phases with formation enthalpy
within 0.8 eV/atom ($\sim$9,000 K) above the convex hull. Phases in this
energy range provide much better statistical information than just ground
states. Here we employ a cluster alignment (CA) method [49, 50, 51] to analyze
the motifs of Fe-centered and Mg-centered clusters in each structure. The
clusters here are defined by a center atom and its first-shell neighbor atoms.
The cluster alignment method identifies the similarity between an as-extracted
cluster and the perfect template cluster. Here we use typical motifs in the
metallic alloys, including FCC-, BCC- and HCP-type clusters as the templates.
We also include common motifs in Fe-O binary compounds which are OCT
(octahedron) and Prism (trigonal prism) [44]. The snapshots of these motifs
are shown in Fig. 7. We first align the atomic cluster against these five
motifs to check its similarity with the CA method. If the cluster does not
match any of the above five polyhedra, it is marked as an ”other” type. In
Fig. 7 we show scatter-type plots of enthalpy differences between these
structures and that of the convex-hull and their respective atomic volumes. We
differentiate the cluster type and O or Fe concentrations with symbols and
colors, respectively. Since we obtain similar results between the analysis of
Fe-centered clusters and Mg-centered clusters, we only show the analysis of
Fe-centered clusters in Fig. 7 and provide the analysis of Mg-centered
clusters in Supplementary Fig. S3.

Figure 7:  Scatter plot of enthalpies above convex-hull (Hd) and volume for
both stable and metastable Fe-Mg-O structures from AGA search. The color bar
in (a) and (b) represent total oxygen concentration as
$O\%=\frac{n(O)}{n(Fe)+n(Mg)+n(O)}\times 100\%$. The color bar in (c) and (d)
represent iron concentration in Fe and Mg as
$Fe\%=\frac{n(Fe)}{n(Fe)+n(Mg)}\times 100\%$. The right panel shows the
template motifs.

At 200 GPa, Figs. 7(a,c) indicate that crystal structures with octahedral type
clusters and 50 mole% oxygen concentration generally have lower energy than
the other structures. While there is no ternary ground state at 200 GPa, these
structures are concentrated within an energy range very close to the convex
hull. Therefore, it is likely that they become ground states at high
temperatures. These structures are all similar to the NaCl-type B1 structure
and have variable Fe/Mg ratios and occupancies, which is essentially
(Mg1-xFex)O, i.e. ferropericlase (xFe $<$ 0.5) or magnesiumwüstite ($x_{Fe}$
$>$ 0.5). This result is consistent with ferropericlase being a dominant phase
not only in the Earth’s lower mantle but at the higher pressures of Super-
Earth’s mantles. It’s worth noting that the lowest enthalpy structure at 200
GPa is the same as the ground-state FeMg3O4 Cmmm at 350 GPa shown in Fig. 4.
Inspecting the higher energy range in Fig. 7(a) and (c), one finds that
oxygen-rich structures generally have lower enthalpies than oxygen-poor
structures. The oxygen-rich structures mainly contain octahedral clusters,
while the oxygen-poor structures can have a greater variety of motifs,
including FCC, BCC, HCP-type clusters. This is mainly because Fe and Mg start
to alloy to form closely packed motifs under unsaturated oxygen conditions.

At 350 GPa, structures with 50% oxygen concentrations still have the lowest
formation enthalpy. A few oxygen-rich phases become more stable and approach
the convex hull energetically compared to 200 GPa. BCC-type clusters start to
appear in these oxygen-rich structures at 350 GPa, indicating that the B2-type
Fe-O clusters are favored at higher pressures over B1-type clusters at lower
pressures. The situation with oxygen-poor structures at 350 GPa is similar to
the one at 200 GPa. We note that at both 200GPa and 350GPa, several motifs
(”other” type) that cannot be classified into the current simple cluster
templates appear. Some of them are due to distortions, while some indeed form
more complex clusters, e.g., the ones in Fig. 2 and Fig. 3.

##  IV Geophysical implications

Our findings on the Fe-Mg-O system at core pressures appear to have some
straightforward geophysical consequences. The Fe-rich side (right corner) of
the ternary phase diagram in Fig. 1 suggests that Fe2Mg and Fe2O can form a
continuous isomorphic solid Fe2(Mg1-xOx) solution. Both end-members are BCC-
like structures at 350 GPa, as shown in Fig. 8. BCC-like Fe2Mg and hcp
$\epsilon$-Fe are likely to inter-alloy and form a eutectic system, with two
coexisting solid phases for some composition-temperature ranges. Small Mg
concentrations might produce hcp-like Fe1-xMgx alloys, but a BCC-like
Fe2+xMg1-x might precipitate and coexist beyond a certain concentration
threshold. The situation is very similar for the Fe2O-Fe system. Therefore,
the Fe-Mg-O system might contain Mg and O dissolved substitutionally in
$\epsilon$-Fe for small Mg and O concentrations, but beyond a certain
concentration threshold BCC-like Fe2+x+y(Mg1/2-xO1/2-y) might precipitate.
BCC-Fe can be stabilized at inner core pressures by alloying with S [52, 53],
and it has been argued, but not confirmed, that BCC iron could be stabilized
at inner core conditions [54]. Therefore, the precipitation of BCC-like
Fe2+x+y(Mg1/2-xO1/2-y) for non-negligible amounts of Mg, O, or both is not a
surprising conclusion.

The ternary phases discovered in the O-rich side (left corner) of the phase
diagram are relevant for the mantle of some Super-Earths. The absence of
stable ternary phases at pressures lower than 228 GPa suggests that stable
phases involving all three elements are solid-solutions of end-member phases
with a small concentration of inter-alloying metals. For example, Fig. 7(a)
shows that at 200 GPa, the low-energy structures are dominated by structures
with octahedral coordination, with more Mg than Fe, and approximately 50% O,
i.e., ferropericlase or B1-type (Mg1-xFex)O. At 350 GPa, the oxygen-rich
ternary phases FeMgO3, Fe2MgO4, FeMg2O4, and FeMg3O4 emerge as ground states
or low-enthalpy phases, besides the B1-type phase. One of them, I-42d FeMg2O4,
has the same structure as I-42d Mg2SiO4, the stable silicate phase predicted
to exist in the mantle of Super-Earths above 500 GPa [17, 43]. Here emerges
the possibility of an I-42d-type Mg2(Si1-xFex)O4 phase, with Fe substitutional
in the Si site, or vice-versa, an unusual type of substitution in the Earth’s
mantle, unless as a coupled Mg-Si substitution. From the chemistry standpoint,
the newly found phases at 350 GPa can all be viewed as combinations of binary
end-members, e.g., FeMgO3 as (MgO)(FeO2), Fe2MgO4 as (MgO)(Fe2O3), FeMg2O4 as
(MgO)2(FeO2), and FeMg3O4 as (MgO)3(FeO). Such stable compositions suggest
other stable stoichiometric phases might be found by exploring combinations of
such end-member compounds, as seen in the Mg-Si-O system, i.e., (MgO)n(SiO2)m
phases [17, 43]. Further AGA searches aiming at these complex compositions are
needed to identify other possible ternary phases in the Mg-Fe-O system.
Finally, O’s greater intermixing with the metallic elements at 350 GPa
suggests that Mg and O abundances might be non-negligible in the Earth’s inner
core. Also, core formation by Fe exsolution from the oxides might be a more
complicated process during Super-Earths’ core formation, or O and Mg might be
more abundant light elements in Super-Earths’ cores.

Figure 8:  Similar BCC-like Fe2O and Fe2Mg ground states at 350GPa.

##  V Conclusion

We use the AGA combined with ab initio calculations to identify high-pressure
structures in the Fe-Mg-O system at 0 K across a wide range of
stoichiometries. This procedure is a crucial preparatory stage for modeling
the system at finite temperatures. At 350 GPa, we identify mechanically stable
phases with FeMg2O4, Fe2MgO4 and FeMg3O4 compositions and one low-enthalpy
phase with FeMgO3 composition. These discoveries lead to the construction of
the ternary phase diagram and convex hull at 350 GPa. While we have not found
any ground-state stoichiometric ternary compound at 200 GPa, the metastable
phases’ analysis indicates that ferropericlase- or magnesiumwüstite-type
phases with 50% oxygen are very close to the convex hull. Oxygen-rich phases
are generally closer to the convex hull than the oxygen-poor phases at all
pressures. Motif analyses show octahedral clusters are energetically favored
at both pressures and BCC-type clusters start to appear in oxygen-rich phases
at 350 GPa. In particular, the nature of iron-rich phases at 350 GPa indicates
that Mg, O, or both simultaneously could stabilize a BCC-type iron alloy at
inner-core pressures.

##  VI acknowledgments

This work was supported primarily by National Science Foundation awards
EAR-1918134 and EAR-1918126 and the Extreme Science and Engineering Discovery
Environment (XSEDE), which is supported by National Science Foundation grant
number ACI-1548562. R.M.W. and Y.S. also acknowledge partial support from the
Department of Energy Theoretical Chemistry Program through grant DOE-
DESC0019759. R.W. and Z.L. were supported by the National Natural Science
Foundation of China (11774324 & 12074362) and the Supercomputing Center of
USTC. Y.F., F.Z., and S.W. were supported by the National Natural Science
Foundation of China (11874307).

## References

  * [1] A. E. Doyle, E. D. Young, B. Klein, B. Zuckerman, and H. E. Schlichting, Science 366, 356 (2019).
  * [2] K. Hirose, S. Labrosse, and J. Hernlund, Annu. Rev. Earth Planet. Sci. 41, 657 (2013).
  * [3] F. Birch, J. Geophys. Res. 57, 227 (1952).
  * [4] F. Birch, J. Geophys. Res. 69, 4377 (1964).
  * [5] L. S. Dubrovinsky, S. K. Saxena, F. Tutti, S. Rekhi, and T. LeBehan, Phys. Rev. Lett. 84, 1720 (2000).
  * [6] J. Shanker, B. P. Singh, and S. K. Srivastava, Phys. Earth Planet. Inter. 147, 333 (2004).
  * [7] A. Dewaele, P. Loubeyre, F. Occelli, M. Mezouar, P. I. Dorogokupets, and M. Torrent, Phys. Rev. Lett. 97, 29 (2006).
  * [8] Y. Kuwayama, G. Morard, Y. Nakajima, K. Hirose, A. Q. R. Baron, S. I. Kawaguchi, T. Tsuchiya, D. Ishikawa, N. Hirao, and Y. Ohishi, Phys. Rev. Lett. 124, 165701 (2020).
  * [9] G. Morard, D. Andrault, D. Antonangeli, and J. Bouchet, Comptes Rendus - Geosci. 346, 130 (2014).
  * [10] D. Alfè, M. J. Gillan, and G. D. Price, Earth Planet. Sci. Lett. 195, 91 (2002).
  * [11] H. Ozawa, F. Takahashi, K. Hirose, Y. Ohishi, and N. Hirao, Science 334, 792 (2011).
  * [12] S. Tateno, K. Hirose, Y. Ohishi, and Y. Tatsumi, Science 330, 359 (2010).
  * [13] M. J. Gillan, D. Alfè, J. Brodholt, L. Vočadlo, and G. D. Price, Reports Prog. Phys. 69, 2365 (2006).
  * [14] F. Zhang and A. R. Oganov, Geophys. Res. Lett. 37, n/a (2010).
  * [15] G. L. Weerasinghe, R. J. Needs, and C. J. Pickard, Phys. Rev. B - Condens. Matter Mater. Phys. 84, 174110 (2011).
  * [16] G. L. Weerasinghe, C. J. Pickard, and R. J. Needs, J. Phys. Condens. Matter 27, 455501 (2015).
  * [17] S. Q. Wu, M. Ji, C. Z. Wang, M. C. Nguyen, X. Zhao, K. Umemoto, R. M. Wentzcovitch, and K. M. Ho, J. Phys. Condens. Matter 26, 035402 (2014).
  * [18] K. Umemoto and R. M. Wentzcovitch, Phys. Rev. Mater. 3, 123601 (2019).
  * [19] L. Zhang, Y. Wang, J. Lv, and Y. Ma, Nat. Rev. Mater. 2, 1 (2017).
  * [20] A. R. Oganov, C. J. Pickard, Q. Zhu, and R. J. Needs, Nat. Rev. Mater. 4, 331 (2019).
  * [21] E. Ohtani and A. E. Ringwood, Earth Planet. Sci. Lett. 71, 85 (1984).
  * [22] D. Alfè, M. J. Gillan, and G. D. Price, J. Chem. Phys. 116, 7127 (2002).
  * [23] H. Huang, Y. Fei, L. Cai, F. Jing, X. Hu, H. Xie, L. Zhang, and Z. Gong, Nature 479, 513 (2011).
  * [24] J. Badro, G. Fiquet, F. Guyot, J. P. Rueff, V. V. Struzhkin, G. Vankó, and G. Monaco, Science 300, 789 (2003).
  * [25] P. Gao, C. Su, S. Shao, S. Wang, P. Liu, S. Liu, and J. Lv, New J. Chem. 43, 17403 (2019).
  * [26] Q. Hu, D. Y. Kim, W. Yang, L. Yang, Y. Meng, L. Zhang, and H.-K. Mao, Nature 534, 241 (2016).
  * [27] M. Nishi, Y. Kuwayama, J. Tsuchiya, and T. Tsuchiya, Nature 547, 205 (2017).
  * [28] A. P. P. Van den Berg, D. A. A. Yuen, K. Umemoto, M. H. G. H. G. Jacobs, and R. M. M. Wentzcovitch, Icarus 317, 412 (2019).
  * [29] X. Zhao, M. C. Nguyen, W. Y. Zhang, C. Z. Wang, M. J. Kramer, D. J. Sellmyer, X. Z. Li, F. Zhang, L. Q. Ke, V. P. Antropov, and K. M. Ho, Phys. Rev. Lett. 112, 045502 (2014).
  * [30] D. M. Deaven and K. M. Ho, Phys. Rev. Lett. 75, 288 (1995).
  * [31] S. M. Foiles, M. I. Baskes, and M. S. Daw, Phys. Rev. B 33, 7983 (1986).
  * [32] X. W. Zhou, R. A. Johnson, and H. N. G. Wadley, Phys. Rev. B 69, 144113 (2004).
  * [33] A. Banerjea and J. R. Smith, Phys. Rev. B 37, 6632 (1988).
  * [34] P. Brommer and F. Gähler, Model. Simul. Mater. Sci. Eng. 15, 295 (2007).
  * [35] P. Brommer and F. Gähler, Philos. Mag. 86, 753 (2006).
  * [36] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * [37] G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996).
  * [38] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * [39] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * [40] H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * [41] A. Togo and I. Tanaka, Scr. Mater. 108, 1 (2015).
  * [42] F. C. Frank and J. S. Kasper, Acta Crystallogr. 11, 184 (1958).
  * [43] K. Umemoto, R. M. Wentzcovitch, S. Wu, M. Ji, C. Z. Wang, and K. M. Ho, Earth Planet. Sci. Lett. 478, 40 (2017).
  * [44] E. Bykova, L. Dubrovinsky, N. Dubrovinskaia, M. Bykov, C. McCammon, S. V. Ovsyannikov, H.-P. Liermann, I. Kupenko, A. I. Chumakov, R. Rüffer, M. Hanfland, and V. Prakapenka, Nat. Commun. 7, 10661 (2016).
  * [45] Q. Zhu, A. R. Oganov, and A. O. Lyakhov, Phys. Chem. Chem. Phys. 15, 7696 (2013).
  * [46] M. Hellenbrandt, Crystallogr. Rev. 10, 17 (2004).
  * [47] D. Andrault and N. Bolfan-Casanova, Phys. Chem. Miner. 28, 211 (2001).
  * [48] S. M. Antao, I. Hassan, W. A. Crichton, and J. B. Parise, Am. Mineral. 90, 1500 (2005).
  * [49] X. W. Fang, C. Z. Wang, Y. X. Yao, Z. J. Ding, and K. M. Ho, Phys. Rev. B 82, 184204 (2010).
  * [50] Y. Sun, F. Zhang, Z. Ye, Y. Zhang, X. Fang, Z. Ding, C. Z. Wang, M. I. Mendelev, R. T. Ott, M. J. Kramer, and K. M. Ho, Sci. Rep. 6, 23734 (2016).
  * [51] S. Ren, Y. Sun, F. Zhang, A. Travesset, C. Z. Wang, and K. M. Ho, ACS Nano 14, 6795 (2020).
  * [52] L. Vočadlo, D. Alfè, M. J. Gillan, I. G. Wood, J. P. Brodholt, and G. D. Price, Nature 424, 536 (2003).
  * [53] P. N. Gavryushkin, Z. I. Popov, K. D. Litasov, A. B. Belonoshko, and A. Gavryushkin, Geophys. Res. Lett. 43, 8435 (2016).
  * [54] R. Hrubiak, Y. Meng, and G. Shen, arXiv 1804.05109, (2018).


# Cascaded Channel Estimation for RIS Assisted mmWave MIMO Transmissions

Yushan Liu, Shun Zhang, _Senior Member, IEEE,_ Feifei Gao, _Fellow, IEEE,_ Jie
Tang, _Senior Member, IEEE,_ and Octavia A. Dobre, _Fellow, IEEE_ Y. Liu and
S. Zhang are with the State Key Laboratory of Integrated Services Networks,
Xidian University, Xi’an 710071, P. R. China (Email:
<EMAIL_ADDRESS>zhangshunsdu@xidian.edu.cn). F. Gao is with
Department of Automation, Tsinghua University, State Key Lab of Intelligent
Technologies and Systems, Tsinghua University, State Key for Information
Science and Technology (TNList) Beijing 100084, P. R. China (Email:
feifeigao@ieee.org). J. Tang is with the School of Electronic and Information
Engineering, South China University of Technology, Guangzhou 510641, China
(E-mail: eejtang@scut.edu.cn). O. A. Dobre is with Faculty of Engineering and
Applied Science, Memorial University, St. John’s NL AIC-5S7, Canada (Email:
odobre@mun.ca).

###### Abstract

Channel estimation is challenging for the reconfigurable intelligence surface
(RIS) assisted millimeter wave (mmWave) communications. Since the number of
coefficients of the cascaded channels in such systems is closely dependent on
the product of the number of base station antennas and the number of RIS
elements, the pilot overhead would be prohibitively high. In this letter, we
propose a cascaded channel estimation framework for an RIS assisted mmWave
multiple-input multiple-output system, where the wideband effect on
transmission model is considered. Then, we transform the wideband channel
estimation into a parameter recovery problem and use a few pilot symbols to
detect the channel parameters by the Newtonized orthogonal matching pursuit
algorithm. Moreover, the Cramer-Rao lower bound on the channel estimation is
introduced. Numerical results show the effectiveness of the proposed channel
estimation scheme.

###### Index Terms:

RIS, mmWave MIMO, wideband effect, channel estimation, NOMP, CRLB.

##  I Introduction

With wide frequency bands, millimeter wave (mmWave) communications can provide
unprecedented gigabits-per-second data rates and satisfy the rapidly growing
transmission speed demand of wireless communications [1]. However, radio
signals in mmWave bands are sensitive to the obstacles and suffer from the
severe path loss. To address this issue, many large-antenna based
technologies, such as massive multiple-input multiple-output (MIMO) [2] and
reconfigurable intelligent surface (RIS) assisted MIMO, have been explored.
Different from radio frequency chain based MIMO systems, RIS is formed as an
artificial planar structure with integrated electronic circuits, and can be
programmed to manipulate an incoming electromagnetic field in a wide variety
of functionalities. It consists of a large number of reconfigurable reflective
elements, which can induce an adjustable independent phase shift on the
incident signal [3]. Recent research results show that RIS assisted MIMO can
achieve similar or even higher communication performance gains with much
smaller hardware cost than massive MIMO.

Generally, RIS usually works in the passive state and has no signal processing
capability [4]. Hence, channel estimation in RIS assisted wireless systems is
more challenging than that in traditional systems. In [5, 6], they examined
the cascaded channel estimation over the RIS aided MIMO systems. In [7], Hu et
al. proposed a novel location information aided channel estimation method,
which substantially reduced the estimation overhead. However, in the large-
scale antenna systems, different antennas at the same sampling time would
receive different time-domain symbols from the same physical path due to the
large propagation delay of electromagnetic waves travelling across the whole
antenna array, which is known as the spatial wideband effect [8]. In such
case, the RIS assisted MIMO channel model in the above researches, which only
considers phase difference and ignores delay difference among the transmission
signals at different RIS elements and BS antennas, are not applicable anymore.
Moreover, the algorithms based on such models, such as for the channel
estimation and parameter recovery, need to be revised.

In this letter, we investigate the wideband channel estimation scheme for an
RIS assisted mmWave MIMO system, which takes the wideband effect into
consideration. The RIS assisted channels are depicted as the functions of
physical parameters, including the angle of arrival/departure (AoA/AoD), the
time delay and the complex gain. The phase and delay differences among the
received signals at different BS antennas and RIS elements are considered.
Then, the frequency response of the received signal is derived and the channel
estimation is formulated as a parameter recovery problem. To achieve the
acquisition of channel parameters, we resort to the Newtonized orthogonal
matching pursuit (NOMP) algorithm. To make the study complete, we also derive
the Cramer-Rao lower bound (CRLB). Finally, the numerical evaluations show the
effectiveness of the proposed scheme.

##  II Wideband Channel Model over RIS Assisted Network

We consider an RIS assisted mmWave MIMO system with a single antenna user, a
BS of $N_{\text{b}}$ antennas and an RIS consisting of $N_{\text{r}}$ passive
elements, as depicted in Fig. 1. The antennas at BS form a uniform linear
array (ULA), and the elements in the RIS form a uniform planar array (UPA)
with $N_{\text{x}}$ and $N_{\text{y}}$ elements along the horizontal and
vertical directions, respectively. Orthogonal frequency division multiplexing
(OFDM) with $N_{\text{c}}$ subcarriers is adopted for combating the multipath
delay spread. Let us denote the transmission bandwidth as $W$, and the
subcarrier spacing is $\Delta\\!f\\!=\\!W\\!/\\!N_{\text{c}}$. The length of
the cyclic prefix is assumed to be longer than that of the maximum multipath
delay plus the maximum delay of the antenna domain.

Figure 1: An RIS assisted mmWave MIMO system.

Suppose that the multipath channel from the user to RIS consists of $P$
incident paths. Denote $\tau_{p,r}^{\text{UR}}$ as the time delay of the
$p$-th path from the user to the $r$-th element of the RIS, where
$p\\!\in\\!\\{0\\!,\cdots\\!,P\\!-\\!1\\}$ and
$r\\!\in\\!\\{0\\!,\cdots\\!,N_{\text{r}}\\!-\\!1\\}$. Note that
$r=iN_{\text{y}}+j$ where
$i\\!\in\\!\\{0\\!,\cdots\\!,N_{\text{x}}\\!-\\!1\\}$ and
$j\\!\in\\!\\{0,\cdots,N_{\text{y}}-1\\}$ are the element indexes along the
horizontal and vertical directions, respectively. Moreover, it can be inferred
that $j\\!=\\!(r)_{N_{\text{y}}}$ and
$i\\!=\\!(r\\!-\\!(r)_{N_{\text{y}}})/N_{\text{y}}$, where
$j\\!=\\!(r)_{N_{\text{y}}}$ denotes the remainder of $r$ divided by
$N_{\text{y}}$. Denote $\phi_{p}$ and $\psi_{p}$ as the elevation and azimuth
AoAs along the $p$-th path at the RIS, respectively. With the RIS structure
and the far-field assumption [8], $\tau_{p,r}^{\text{UR}}$ can be denoted as

$\displaystyle\tau_{p,r}^{\text{UR}}\\!=\\!\tau_{p}\\!+\\!\frac{d\big{(}\frac{r-(r)\\!_{N_{\text{y}}}}{N_{\text{y}}}\sin(\phi_{p})\\!\sin(\psi_{p})\\!\\!+\\!\\!(r)\\!_{N_{\text{y}}}\sin(\phi_{p})\\!\cos(\psi_{p})\big{)}}{c},$
(1)

where $d$ is the antenna spacing and $c$ is the speed of light. Denote the
complex path gain of the $p$-th path as $g_{p}^{\text{UR}}$. Then, the impulse
response of the multipath channel from the user to the $r$-th element at the
RIS can be expressed as

$\displaystyle
h^{\text{UR},\text{T}}_{r}(t)=\sum_{p=0}^{P-1}\underbrace{g^{\text{UR}}_{p}e^{-\jmath
2\pi f_{\text{c}}\tau_{p}}}_{\bar{g}^{\text{UR}}_{p}}e^{-\jmath
2\pi\varpi_{r}(\phi_{p},\psi_{p})}\delta(t-\tau_{p,r}^{\text{UR}}),$ (2)

where $f_{\text{c}}$ is the carrier frequency, $\delta(\cdot)$ is the Delta
function, $\bar{g}^{\text{UR}}_{p}$ is the equivalent complex gain,
$\varpi_{r}(\phi,\psi)\triangleq\big{(}d(\frac{r-(r)_{N_{\text{y}}}}{N_{\text{y}}}\sin(\phi)\\!\sin(\psi)\\!+\\!(r)_{N_{\text{y}}}\sin(\phi)\\!\cos(\psi))\big{)}/\lambda_{\text{c}}$
is the normalized angle at the RIS and $\lambda_{\text{c}}$ is the carrier
wavelength [8].

The direct path from the user to BS may be blocked by possible obstacles such
as buildings and trees [9]. Hence, we mainly focus on the RIS assisted link.
Generally, BS and RIS are considered to be located in an environment with
limited local scattering, which would cause the MIMO link between BS and RIS
to be a light-of-sight (LoS) [10]. Thus, we define $g^{\text{RB}}$ as the
complex path gain, $\theta_{\text{b}}$ as the AoA at BS,
$\phi_{\text{r}}/\psi_{\text{r}}$ as the elevation/azimuth AoDs at RIS and
$\tau_{b,r}^{\text{RB}}=bd\sin(\theta_{\text{b}})/c+\varpi_{r}(\phi_{\text{r}},\psi_{\text{r}})/f_{\text{c}}$
as the time delay from the $r$-th element of the RIS to the $b$-th antenna of
the BS, where $c$ is the speed of light. Then, the impulse response of the LoS
channel from the $r$-th RIS element to the $b$-th BS antenna can be expressed
as

$\displaystyle H_{b,r}^{\text{RB},\text{T}}(t)=g^{\text{RB}}e^{-\jmath 2\pi
f_{\text{c}}b\frac{d\sin(\theta_{\text{b}})}{c}}e^{-\jmath
2\pi\varpi_{r}(\phi_{\text{r}},\psi_{\text{r}})}\delta(t-\tau_{b,r}^{\text{RB}}).$
(3)

##  III NOMP based Wideband Channel Estimation for RIS Assisted MIMO

###  III-A Wideband Effect on Transmission Model

By taking the Fourier transform of (2) and stacking it from RIS’s different
antennas into a $N_{\text{r}}\times 1$ vector, the frequency response between
the RIS and user can be obtained as

$\displaystyle\mathbf{h}^{\text{UR},\text{F}}(f)=\sum_{p=0}^{P-1}\bar{g}^{\text{UR}}_{p}\mathbf{a}_{\text{R}}(f,\phi_{p},\psi_{p})e^{-\jmath
2\pi f\tau_{p}},$ (4)

where $\mathbf{a}_{\text{R}}(f,\phi,\psi)=[1,\cdots,e^{-\jmath
2\pi(1+\frac{f}{f_{\text{c}}})\varpi_{N_{\text{r}-1}}(\phi,\psi)}]^{T}$ is the
array steering vector at the RIS. By a similar mathematical manipulation on
(3), the frequency response between the BS and RIS can be expressed as

$\displaystyle\mathbf{H}^{\text{RB},\text{F}}(f)=g^{\text{RB}}\mathbf{a}_{\text{B}}(f,\theta_{\text{b}})\mathbf{a}^{T}_{\text{R}}(f,\phi_{\text{r}},\psi_{\text{r}}),$
(5)

where $\mathbf{a}_{\text{B}}(f,\theta)=[1,\cdots,e^{-\jmath
2\pi(1+\frac{f}{f_{\text{c}}})(N_{\text{b}}-1)\frac{d\sin(\theta)}{\lambda_{\text{c}}}}]^{T}$
is the array steering vector at the BS.

Before proceeding, let us denote the RIS phase shift vector at the $m$-th OFDM
pilot symbol as
$\boldsymbol{\rho}_{m}\triangleq[e^{\jmath\rho_{m,0}},\cdots,e^{\jmath\rho_{m,N_{\text{r}}-1}}]^{T}$,
where $m\in\\{0,\cdots,M-1\\}$, $M$ is the number of pilot symbols during
channel estimation, and $\rho_{m,r}\in[0,2\pi)$ represents the phase shift of
the $r$-th RIS element at the $m$-th pilot symbol. Then, the operation of RIS
can be described by the diagonal matrix
$\boldsymbol{\Omega}_{m}=\text{diag}\\{\boldsymbol{\rho}_{m}\\}\in\mathbb{C}^{N_{\text{r}}\times
N_{\text{r}}}$. Note that $\boldsymbol{\rho}_{m}$ is assumed to be constant
during each pilot symbol. Then, the frequency response of BS’s received signal
at the $m$-th pilot symbol can be expressed as

$\displaystyle\mathbf{y}_{m}(f)$
$\displaystyle=\sum_{p=0}^{P-1}\underbrace{g^{\text{RB}}\bar{g}^{\text{UR}}_{p}}_{g_{p}}\mathbf{a}_{\text{B}}(f,\theta_{\text{b}})\mathbf{a}^{T}_{\text{R}}(f,\phi_{\text{r}},\psi_{\text{r}})\boldsymbol{\Omega}_{m}\mathbf{a}_{\text{R}}(f,\phi_{p},\psi_{p})$
$\displaystyle\kern 60.0pt\times s_{m}(f)e^{-\jmath 2\pi
f\tau_{p}}+\mathbf{v}_{m}(f),$ (6)

where $g_{p}$ is defined as the cascaded complex gain, $s_{m}(f)$ is the
frequency-domain signal from the user, and $\mathbf{v}_{m}(f)$ is the additive
Gaussian noise with zero mean and covariance matrix
$\sigma_{\text{v}}^{2}\mathbf{I}_{N_{\text{b}}}$.

Assume that $L$ out of $N_{\text{c}}$ subcarriers are exclusively assigned to
the user as pilots whose index set is denoted by
$\mathcal{L}=\\{l_{0},\cdots,l_{L-1}\\}$. Without loss of generality, we
assume that the values of all pilots are 1, i.e., $s_{m}(k\Delta f)=1$,
$m\in\\{0,\cdots,M-1\\}$, $k\in\mathcal{L}$. Then, we collect the received
pilot vectors $\mathbf{y}_{0}{(k\Delta f)},\cdots,\mathbf{y}_{M-1}{(k\Delta
f)}$ into an $MN_{\text{b}}\times 1$ vector at the $k$-th subcarrier as

$\displaystyle\mathbf{y}(k\Delta f)$
$\displaystyle\triangleq[\mathbf{y}_{0}^{T}(k\Delta
f),\cdots,\mathbf{y}_{M-1}^{T}(k\Delta f)]^{T}$
$\displaystyle=\sum_{p=0}^{P-1}\\!g_{p}\\!\Big{(}\\!\underbrace{\left[\begin{array}[]{c}\\!\\!\\!\mathbf{a}^{T}_{\text{R}}(k\Delta
f,\phi_{\text{r}},\psi_{\text{r}})\boldsymbol{\Omega}_{0}\\!\\!\\!\\\
\vdots\\\ \\!\\!\\!\mathbf{a}^{T}_{\text{R}}(k\Delta
f,\phi_{\text{r}},\psi_{\text{r}})\boldsymbol{\Omega}_{M-1}\\!\\!\\!\end{array}\right]}_{\bar{\mathbf{A}}_{k}(\phi_{\text{r}},\psi_{\text{r}})\in\mathbb{C}^{M\times
N_{\text{r}}}}\\!\mathbf{a}_{\text{R}}(k\Delta f,\phi_{p},\psi_{p})\\!\Big{)}$
(10) $\displaystyle\kern 10.0pt\otimes\mathbf{a}_{\text{B}}(k\Delta
f,\theta_{\text{b}})e^{-\jmath 2\pi k\Delta
f\tau_{p}}\\!+\\!\mathbf{v}(k\Delta f),\quad k\in\mathcal{L},$ (11)

where $\bar{\mathbf{A}}_{k}(\phi_{\text{r}},\psi_{\text{r}})$ is defined above
and $\mathbf{v}(k\Delta f)=[\mathbf{v}_{0}^{T}(k\Delta
f),\cdots,\mathbf{v}_{M-1}^{T}(k\Delta f)]^{T}$ is the noise vector. Define
$\mathbf{f}(\theta_{\text{b}},\phi_{\text{r}},\psi_{\text{r}},\phi_{p},\psi_{p},\tau_{p})\in\mathbb{C}^{LMN_{\text{b}}\times
1}$ as

$\displaystyle\mathbf{f}(\theta_{\text{b}},\phi_{\text{r}},\psi_{\text{r}},\phi_{p},\psi_{p},\tau_{p})\\!=\\!\bigg{[}\Big{[}\big{(}\bar{\mathbf{A}}_{l_{0}}(\phi_{\text{r}},\psi_{\text{r}})\mathbf{a}_{\text{R}}(l_{0}\Delta
f,\phi_{p},\psi_{p})\big{)}\\!\otimes\\!\mathbf{a}_{\text{B}}(l_{0}\Delta
f,\theta_{\text{b}})e^{-\jmath 2\pi l_{0}\Delta f\tau_{p}}\Big{]}^{T},\cdots,$
$\displaystyle\kern
88.0pt\Big{[}\big{(}\bar{\mathbf{A}}_{l_{L-1}}(\phi_{\text{r}},\psi_{\text{r}})\mathbf{a}_{\text{R}}(l_{L-1}\Delta
f,\phi_{p},\psi_{p})\big{)}\\!\otimes\\!\mathbf{a}_{\text{B}}(l_{L-1}\Delta
f,\theta_{\text{b}})e^{-\jmath 2\pi l_{L-1}\Delta
f\tau_{p}}\Big{]}^{T}\bigg{]}^{T}.$ (12)

By collecting $\mathbf{y}(k\Delta f)$ at different subcarriers, we obtain

$\displaystyle\mathbf{y}$ $\displaystyle\triangleq[\mathbf{y}^{T}(l_{0}\Delta
f),\cdots,\mathbf{y}^{T}(l_{L-1}\Delta f)]^{T}$
$\displaystyle=\sum_{p=0}^{P-1}g_{p}\mathbf{f}(\theta_{\text{b}},\phi_{\text{r}},\psi_{\text{r}},\phi_{p},\psi_{p},\tau_{p})+\mathbf{v},$
(13)

where $\mathbf{v}\triangleq[\mathbf{v}^{T}(l_{0}\Delta
f),\cdots,\mathbf{v}^{T}(l_{L-1}\Delta
f)]^{T}\in\mathbb{C}^{LMN_{\text{b}}\times 1}$ is the corresponding noise
vector.

We are interested in estimating the cascaded channels at the BS from (13),
which can be transformed into the problem of parameter recovery. Without loss
of generality, we assume that the locations of the BS and RIS are known, which
implies that the angles $\theta_{\text{b}}$, $\phi_{\text{r}}$ and
$\psi_{\text{r}}$ can be completely determined by the geometric positions of
BS and RIS. Then, $\\{\theta_{\text{b}},\phi_{\text{r}},\psi_{\text{r}}\\}$
are constants and can be omitted in
$\bar{\mathbf{A}}_{k}(\phi_{\text{r}},\psi_{\text{r}})$ and
$\mathbf{f}(\theta_{\text{b}},\phi_{\text{r}},\psi_{\text{r}},\phi_{p},\psi_{p},\tau_{p})$.
In order to estimate the RIS assisted channel, we will resort to the NOMP
algorithm to capture the parameter set
$\\{g_{p},\phi_{\text{p}},\psi_{\text{p}},\tau_{p}\\}_{p=0}^{P-1}$.

###  III-B NOMP Algorithm

The NOMP algorithm is divided into five steps: Greedy searching, precise
searching, single refinement, cyclic refinement and cascaded gain updating.

  1. 1.

Greedy Searching: The over-sampled grids along the horizontal AoA, vertical
AoA and delay are respectively divided with sampling rate $\eta_{\phi}$,
$\eta_{\psi}$ and $\eta_{\tau}$, where the range of actual horizontal AoA,
vertical AoA and time delay are respectively $[-\pi/2,\pi/2)$, $[0,\pi/2)$ and
$[0,1/\Delta f)$. Hence, the codeword corresponding to (12) is expressed as

$\displaystyle\mathbf{f}(\bar{\phi},\bar{\psi},\bar{\tau})\\!$
$\displaystyle=\\!\bigg{[}\Big{[}\big{(}\bar{\mathbf{A}}_{l_{0}}\mathbf{a}_{\text{R}}(l_{0}\Delta
f,\bar{\phi},\bar{\psi})\big{)}\\!\otimes\\!\mathbf{a}_{\text{B}}(l_{0}\Delta
f,\theta_{\text{b}})e^{-\jmath 2\pi l_{0}\Delta
f\bar{\tau}}\Big{]}^{T}\\!\\!,\\!\cdots\\!,\\!$ $\displaystyle\kern
50.0pt\Big{[}\big{(}\bar{\mathbf{A}}_{l_{L\\!-\\!1}}\mathbf{a}_{\text{R}}(l_{L\\!-\\!1}\Delta
f,\bar{\phi},\bar{\psi})\big{)}\\!\otimes\\!\mathbf{a}_{\text{B}}(l_{L\\!-\\!1}\Delta
f,\theta_{\text{b}})e^{-\jmath 2\pi l_{L\\!-\\!1}\Delta
f\bar{\tau}}\Big{]}^{T}\bigg{]}^{T}.$ (14)

where
$\bar{\phi}\in\\{-\frac{\pi}{2},\cdots,-\frac{\pi}{2}+\frac{(\eta_{\phi}N_{\text{x}}-1)\pi}{\eta_{\phi}N_{\text{x}}}\\}$,
$\bar{\psi}\in\\{0,\cdots,\frac{(\eta_{\psi}N_{\text{y}}-1)\pi}{2\eta_{\psi}N_{\text{y}}}\\}$
and
$\bar{\tau}\in\\{0,\cdots,\frac{\eta_{\tau}N_{\text{c}}-1}{\eta_{\tau}N_{\text{c}}\Delta
f}\\}$.

At the beginning of the $i$-th iteration, the residual noisy mixture
$\mathbf{y}_{\text{e}}^{i}$ is calculated by

$\displaystyle\mathbf{y}_{\text{e}}^{i}=\mathbf{y}-\sum_{p=0}^{i-1}\hat{g}_{p}\mathbf{f}(\hat{\phi}_{p},\hat{\psi}_{p},\hat{\tau}_{p}),$
(15)

where
$\\{\hat{g}_{p},\hat{\phi}_{p},\hat{\psi}_{p},\hat{\tau}_{p}\\}_{p=0}^{i-1}$
are the estimated parameters in the previous iterations. The course estimated
$\hat{\phi}_{i}$, $\hat{\psi}_{i}$, $\hat{\tau}_{i}$ are obtained by greedily
searching the grid points as

$\displaystyle(\hat{\phi}_{i},\hat{\psi}_{i},\hat{\tau}_{i})=\underset{(\bar{\phi},\bar{\psi},\bar{\tau})}{\arg\max}\frac{|\mathbf{f}^{H}(\bar{\phi},\bar{\psi},\bar{\tau})\mathbf{y}_{\text{e}}^{i}|^{2}}{\|\mathbf{f}(\bar{\phi},\bar{\psi},\bar{\tau})\|^{2}},$
(16)

where $\hat{\phi}_{i}$, $\hat{\psi}_{i}$ and $\hat{\tau}_{i}$ are the
estimated horizontal AoA, vertical AoA and time delay of the $i$-th component
path, respectively.

  2. 2.

Precise Searching: After the simultaneously greedy searching for the three
parameters, we perform a precise searching process near the results of the
greedy searching step. Notice that the implementation of this process is
similar to that of the greedy searching step and the details are omitted here
due to space limitation. Then, the estimated parameters
$\\{\hat{\phi}_{i},\hat{\psi}_{i},\hat{\tau}_{i}\\}$ can be obtained through a
similar operation in (16). Afterwards, the coarse estimation of the cascaded
gain $\hat{g}_{i}$ can be obtained as

$\displaystyle\hat{g}_{i}=\frac{\mathbf{f}^{H}(\hat{\phi}_{i},\hat{\psi}_{i},\hat{\tau}_{i})\mathbf{y}_{\text{e}}^{i}}{\|\mathbf{f}(\hat{\phi}_{i},\hat{\psi}_{i},\hat{\tau}_{i})\|^{2}}.$
(17)

  3. 3.

Single Refinement: We will resort to the extended Newton method and refine
$\hat{\phi}_{i}$, $\hat{\psi}_{i}$, $\hat{\tau}_{i}$ and $\hat{g}_{i}$, and
$R_{\text{s}}$ iterations are executed in this step. The goal of the
refinement step is to minimize the power of the new residual
$\|\mathbf{y}_{\text{e}}^{i}-g\mathbf{f}(\phi,\psi,\tau)\|^{2}$ in the $i$-th
iteration. Hence, the target is to maximize

$\displaystyle\mathcal{N}(\phi,\psi,\tau)\\!=\\!2\mathcal{R}\\{(\mathbf{y}_{\text{e}}^{i})^{H}g\mathbf{f}(\phi,\psi,\tau)\\}\\!-\\!\|g\mathbf{f}(\phi,\psi,\tau)\|^{2}.$
(18)

Then, the refined estimations of $\hat{\phi}_{i}$, $\hat{\psi}_{i}$ and
$\hat{\tau}_{i}$ can be expressed as

$\displaystyle\left[\begin{array}[]{l}\\!\\!\hat{\phi}_{i}\\!\\!\\!\\!\\\
\\!\\!\hat{\psi}_{i}\\!\\!\\!\\!\\\
\\!\\!\hat{\tau}_{i}\\!\\!\\!\\!\end{array}\right]\\!\\!=\\!\\!\left[\begin{array}[]{l}\\!\\!\hat{\phi}_{i}\\!\\!\\!\\!\\\
\\!\\!\hat{\psi}_{i}\\!\\!\\!\\!\\\
\\!\\!\hat{\tau}_{i}\\!\\!\\!\\!\end{array}\right]\\!\\!\\!-\\!\\!\ddot{\boldsymbol{\mathcal{N}}}^{-1}\\!\left(\hat{g}_{i},\hat{\phi}_{i},\hat{\psi}_{i},\hat{\tau}_{i}\right)\dot{\boldsymbol{\mathcal{N}}}\\!\left(\hat{g}_{i},\hat{\phi}_{i},\hat{\psi}_{i},\hat{\tau}_{i}\right),$
(25)

where
$\dot{\boldsymbol{\mathcal{N}}}\left(g,\phi,\psi,\tau\right)=[\frac{\partial\mathcal{N}}{\partial\phi},\frac{\partial\mathcal{N}}{\partial\psi},\frac{\partial\mathcal{N}}{\partial\tau}]^{T}$
is the first-order partial derivative vector, and

$\displaystyle\ddot{\boldsymbol{\mathcal{N}}}(g,\phi,\psi,\tau)=\left[\begin{array}[]{lll}\frac{\partial^{2}\mathcal{N}}{\partial\phi^{2}}&\frac{\partial^{2}\mathcal{N}}{\partial\phi\partial\psi}&\frac{\partial^{2}\mathcal{N}}{\partial\phi\partial\tau}\\\
\frac{\partial^{2}\mathcal{N}}{\partial\psi\partial\phi}&\frac{\partial^{2}\mathcal{N}}{\partial\psi^{2}}&\frac{\partial^{2}\mathcal{N}}{\partial\psi\partial\tau}\\\
\frac{\partial^{2}\mathcal{N}}{\partial\tau\partial\phi}&\frac{\partial^{2}\mathcal{N}}{\partial\tau\partial\psi}&\frac{\partial^{2}\mathcal{N}}{\partial\tau^{2}}\end{array}\right]$
(29)

is the second-order partial derivative matrix. According to (18), we can write
the first-order partial derivatives of $\mathcal{N}(\phi,\psi,\tau)$ as
$\frac{\partial\mathcal{N}}{\partial
x}=2\mathcal{R}\left\\{g\left(\mathbf{y}_{\text{e}}^{i}-g\mathbf{f}\right)^{H}\frac{\partial\mathbf{f}}{\partial
x}\right\\}$, where $x$ can be $\phi$, $\psi$ and $\tau$. The second-order
partial derivative of $\mathcal{N}(\phi,\psi,\tau)$ can be calculated as
$\frac{\partial^{2}\mathcal{N}}{\partial x_{1}\partial
x_{2}}=2\mathcal{R}\left\\{g\left(\mathbf{y}_{\text{e}}^{i}-g\mathbf{f}\right)^{H}\frac{\partial^{2}\mathbf{f}}{\partial
x_{1}\partial x_{2}}-\left|g\right|^{2}\frac{\partial\mathbf{f}^{H}}{\partial
x_{2}}\frac{\partial\mathbf{f}}{\partial x_{1}}\right\\}$, where $x_{1}$ and
$x_{2}$ can be $\phi$, $\psi$ and $\tau$. Due to the space limitation, the
partial derivatives are omitted here.

By carrying out (25), $\\{\hat{\phi}_{i},\hat{\psi}_{i},\hat{\tau}_{i}\\}$ are
refined. Then, the cascaded gain is also updated according to (17).

  4. 4.

Cyclic Refinement: After the single refinement step for the parameters of the
current iteration, $R_{\text{c}}$ iterations of cyclically refinement are
taken into consideration to further perfect the estimations
$\\{\hat{g}_{i},\hat{\phi}_{i},\hat{\psi}_{i},\hat{\tau}_{i}\\}_{p=0}^{i-1}$
of the previous iterations. Similar to the single refinement step, the
extended Newton method is also utilized, and the accurate estimations
$\\{\hat{g}_{i}$, $\hat{\phi}_{i}$,
$\hat{\psi}_{i},\hat{\tau}_{i}\\}_{p=0}^{i}$ can be obtained. The derivation
is omitted due to space limitation.

  5. 5.

Cascaded Gain Updating: Based on
$\\{\hat{\phi}_{p},\hat{\psi}_{p},\hat{\tau}_{p}\\}_{p=0}^{i}$ obtained in the
previous iterations, we can further update the cascaded gains through the
least square algorithm as

$\displaystyle[\hat{g}_{0},\cdots,\hat{g}_{i}]^{T}=(\mathbf{F}_{i}^{H}\mathbf{F}_{i})^{-1}\mathbf{F}_{i}^{H}\mathbf{y},$
(30)

where
$\mathbf{F}_{i}=[\mathbf{f}(\hat{\phi}_{0},\hat{\psi}_{0},\hat{\tau}_{0}),\cdots,\mathbf{f}(\hat{\phi}_{i},\hat{\psi}_{i},\hat{\tau}_{i})]$.

As the NOMP runs, the power of the residual $\mathbf{y}_{\text{e}}^{i}$
decreases after each iteration. If the extracted parameters are accurate
enough, the power of the residual is reduced to the power of the noise in the
end, i.e., $\|\mathbf{y}^{i}_{\text{e}}\|^{2}\approx\|\mathbf{v}\|^{2}$.
Hence, the NOMP algorithm terminates when
$|\mathbf{f}^{H}(\bar{\phi},\bar{\psi},\bar{\tau})\mathbf{y}_{\text{e}}^{i}|^{2}<\epsilon$
for all possible $(\bar{\phi},\bar{\psi},\bar{\tau})$. The stopping criterion
threshold $\epsilon$ is chosen from the false alarm rate
$P_{\text{fa}}\\!=\\!P\\{\|\mathbf{v}\|^{2}_{\infty}\\!>\\!\epsilon\\}\\!=\\!1\\!-\\!(1\\!-\\!\exp(-\\!\epsilon/\sigma^{2}_{\text{v}}))^{LMN_{\text{b}}}$
as
$\epsilon\\!=\\!-\sigma_{\text{v}}^{2}\ln(1\\!-\\!(1\\!-\\!P_{\text{fa}})^{1/LMN_{\text{b}}})$
[11].

###  III-C Performance Analysis

To evaluate the performance of the NOMP algorithm, we derive the CRLB. Let us
define the $5P\times 1$ unknown parameter vector
$\boldsymbol{\xi}=[\boldsymbol{\xi}^{T}_{1},\cdots,\boldsymbol{\xi}^{T}_{P}]^{T}$,
where $\boldsymbol{\xi}_{p}=[|g_{p}|,\angle
g_{p},\phi_{p},\psi_{p},\tau_{p}]^{T}\in\mathbb{R}^{5\times 1}$ denotes the
unknown parameters of the $p$-th path, and $|g_{p}|$, $\angle g_{p}$
separately represent the amplitude and phase of the complex-valued path gain.
Then, for an unbiased estimator, the estimation variance is bounded by CRLB,
which is the inverse of the $5P\times 5P$ Fisher information matrix (FIM)
$\boldsymbol{\mathcal{F}}(\boldsymbol{\xi})$ [12] defined as

$\displaystyle[\boldsymbol{\mathcal{F}}(\boldsymbol{\xi})]_{i,j}=\mathbb{E}_{\mathbf{y}\mid\boldsymbol{\xi}}\left\\{\frac{\partial\ln
p(\mathbf{y}\mid\boldsymbol{\xi})}{\partial\xi_{i}}\frac{\partial\ln
p(\mathbf{y}\mid\boldsymbol{\xi})}{\partial\xi_{j}}\right\\},$ (31)

where $p(\mathbf{y}\mid\boldsymbol{\xi})$ is the likelihood function of
$\mathbf{y}$ conditioned on $\boldsymbol{\xi}$, and the expectation is taken
over the noise distribution. The $5P\times 5P$ FIM can be sliced into $P^{2}$
submatrices as

$\displaystyle\boldsymbol{\mathcal{F}}(\boldsymbol{\xi})=\left[\begin{array}[]{ccc}\boldsymbol{\mathcal{F}}(\boldsymbol{\xi}_{1},\boldsymbol{\xi}_{1})&\cdots&\boldsymbol{\mathcal{F}}(\boldsymbol{\xi}_{1},\boldsymbol{\xi}_{P})\\\
\vdots&\ddots&\vdots\\\
\boldsymbol{\mathcal{F}}(\boldsymbol{\xi}_{P},\boldsymbol{\xi}_{1})&\cdots&\boldsymbol{\mathcal{F}}(\boldsymbol{\xi}_{P},\boldsymbol{\xi}_{P})\end{array}\right],$
(35)

where $\boldsymbol{\mathcal{F}}(\boldsymbol{\xi}_{i},\boldsymbol{\xi}_{j})$ is
a $5\times 5$ matrix. Generally, a closed-form CRLB analysis for multipath
channel estimation is hard to obtain [13]. Moreover, if two or more paths have
an extremely close angle and delay, the rank deficiency of
$\boldsymbol{\mathcal{F}}(\boldsymbol{\xi})$ appears, which causes the
determinant of $\boldsymbol{\mathcal{F}}(\boldsymbol{\xi})$ to be close to
$0$. However, the number of paths to the BS is very small and these paths are
characterized by separation in delay in wideband mmWave systems. Hence, we
assume that $\boldsymbol{\mathcal{F}}(\boldsymbol{\xi})$ is nonsingular and
can be transformed into a block diagonal matrix [12], and each submatrix on
the diagonal of $\boldsymbol{\mathcal{F}}(\boldsymbol{\xi})$ can be written as

$\displaystyle\boldsymbol{\mathcal{F}}(\\!\boldsymbol{\xi}\\!_{p}\\!,\\!\boldsymbol{\xi}\\!_{p}\\!)=\\!\\!\left[\begin{array}[]{lllll}\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!|g\\!_{p}\\!|\\!,\\!|g\\!_{p}\\!|\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!|g\\!_{p}\\!|,\\!\angle\\!g\\!_{p}\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!|g_{p}\\!|,\\!\phi\\!_{p}\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!|g_{p}\\!|,\\!\psi\\!_{p}\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!|g_{p}\\!|,\\!\tau\\!_{p}\\!)\\!\\!\\!\\\
\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\angle\\!g\\!_{p}\\!,\\!|g\\!_{p}\\!|\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\angle\\!g\\!_{p}\\!,\\!\angle\\!g\\!_{p}\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\angle\\!g\\!_{p}\\!,\\!\phi\\!_{p}\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\angle\\!g\\!_{p}\\!,\\!\psi\\!_{p}\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\angle\\!g\\!_{p}\\!,\\!\tau\\!_{p}\\!)\\!\\!\\!\\\
\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\phi\\!_{p}\\!,\\!|g\\!_{p}\\!|\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\phi\\!_{p}\\!,\\!\angle\\!g\\!_{p}\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\phi\\!_{p}\\!,\\!\phi\\!_{p}\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\phi\\!_{p}\\!,\\!\psi\\!_{p}\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\phi\\!_{p}\\!,\\!\tau\\!_{p}\\!)\\!\\!\\!\\\
\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\psi\\!_{p}\\!,\\!|g\\!_{p}\\!|\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\psi\\!_{p}\\!,\\!\angle\\!g\\!_{p}\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\psi\\!_{p}\\!,\\!\phi\\!_{p}\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\psi\\!_{p}\\!,\\!\psi\\!_{p}\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\psi\\!_{p}\\!,\\!\tau\\!_{p}\\!)\\!\\!\\!\\\
\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\tau_{p}\\!,\\!|g_{p}\\!|\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\tau\\!_{p}\\!,\\!\angle\\!g\\!_{p}\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\tau\\!_{p}\\!,\\!\phi\\!_{p}\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\tau\\!_{p}\\!,\\!\psi\\!_{p}\\!)\\!\\!\\!&\\!\\!\\!\boldsymbol{\mathcal{F}}(\\!\tau\\!_{p}\\!,\\!\tau\\!_{p}\\!)\\!\\!\\!\end{array}\right].$
(41)

Considering that the inverse of a block diagonal matrix is a block diagonal
matrix with the inverse of the original blocks on its diagonal, we can
separately calculate the inverse matrix of each submatrix. Thus,
$\boldsymbol{\mathcal{F}}(\boldsymbol{\xi}_{p},\boldsymbol{\xi}_{p})$ can be
derived as in Appendix and the FIM
$\boldsymbol{\mathcal{F}}(\boldsymbol{\xi})$ can be obtained.

Let us define the channel statement information (CSI) as

$\displaystyle\widetilde{H}_{k,m,b}=\sum_{p=0}^{P-1}|g_{p}|e^{\jmath\angle
g_{p}}\breve{A}_{k,m,b,p}e^{-\jmath 2\pi k\Delta f\tau_{p}}.$ (42)

From [13], the lower bound of the CSI can be obtained via the transformation
vector $\partial\widetilde{H}_{k,m,b}/\partial\boldsymbol{\xi}$ as

$\displaystyle\text{LB}_{k,m,b}(\boldsymbol{\xi})=\left(\frac{\partial\widetilde{H}_{k,m,b}}{\partial\boldsymbol{\xi}}\right)^{H}\boldsymbol{\mathcal{F}}^{-1}(\boldsymbol{\xi})\frac{\partial\widetilde{H}_{k,m,b}}{\partial\boldsymbol{\xi}},$
(43)

where the entries of $\partial\widetilde{H}_{k,m,b}/\partial\boldsymbol{\xi}$
are given as

$\displaystyle\frac{d\widetilde{H}_{k,m,b}}{d|g_{p}|}\\!\\!$
$\displaystyle=\\!\\!e^{\jmath\angle g_{p}}\breve{A}_{k,m,b,p}e^{-\jmath 2\pi
k\Delta f\tau_{p}},\quad\frac{d\widetilde{H}_{k,m,b}}{d\angle
g_{p}}\\!\\!=\\!\\!\jmath|g_{p}|e^{\jmath\angle
g_{p}}\breve{A}_{k,m,b,p}e^{-\jmath 2\pi k\Delta f\tau_{p}},$
$\displaystyle\frac{d\widetilde{H}_{k,m,b}}{d\phi_{p}}\\!\\!$
$\displaystyle=\\!\\!|g_{p}|e^{\jmath\angle g_{p}}e^{-\jmath 2\pi k\Delta
f\tau_{p}}\frac{\partial\breve{A}_{k,m,b,p}}{\partial\phi_{p}},\quad\frac{d\widetilde{H}_{k,m,b}}{d\psi_{p}}\\!\\!=\\!\\!|g_{p}|e^{\jmath\angle
g_{p}}e^{-\jmath 2\pi k\Delta
f\tau_{p}}\frac{\partial\breve{A}_{k,m,b,p}}{\partial\psi_{p}},$
$\displaystyle\frac{d\widetilde{H}_{k,m,b}}{d\tau_{p}}\\!\\!$
$\displaystyle=\\!\\!-\jmath 2\pi k\Delta f|g_{p}|e^{\jmath\angle
g_{p}}\breve{A}_{k,m,b,p}e^{-\jmath 2\pi k\Delta f\tau_{p}}.$ (44)

Finally, the CSI error (43) can be constructed by solving the inverse of
$\boldsymbol{\mathcal{F}}(\boldsymbol{\xi})$, where the inversion process of
$\boldsymbol{\mathcal{F}}(\boldsymbol{\xi})$ is omitted due to space
limitations.

##  IV Simulation Results

In this section, we evaluate the performance of the proposed scheme for the
wideband channel estimation over RIS assisted mmWave MIMO. The system
parameter settings are as follows. The carrier frequency is $f_{\text{c}}=28$
GHz, the bandwidth is $W=600$ MHz, and the number of subcarriers is
$N_{\text{c}}=512$. The RIS assisted mmWave channel parameter settings are:
the number of antennas at BS is $N_{\text{b}}=64$, the antenna spacing is
$d=\lambda_{\text{c}}/2$, the number of paths is $P=5$, and $\tau_{p}$,
$\phi_{p}$ and $\psi_{p}$ are uniformly distributed within $[0,32/W)$,
$[-\pi/2,\pi/2)$ and $[-\pi/2,\pi/2)$, respectively. Moreover, during the
channel estimation, the under-sampling rates in the greedy searching are set
as $\eta_{\phi}=2$, $\eta_{\psi}=2$ and $\eta_{\tau}=2$, the over-sampling
rates in the precise searching are set as $\eta_{\phi}=4$, $\eta_{\psi}=4$ and
$\eta_{\tau}=4$, and the iteration number of the single refinement and cyclic
refinement are set as $R_{\text{s}}=R_{\text{c}}=5$.

The signal-to-noise ratio (SNR) is expressed as
$10\log_{10}\sigma_{\text{t}}^{2}/\sigma_{\text{v}}^{2}$, where
$\sigma_{\text{t}}^{2}$ is the average power of the effective signal. Here, we
use the normalized mean square error (NMSE) for the channel parameters as
performance metric, which is defined as
$\text{NMSE}_{\mathbf{x}}\\!\\!\\!=\\!\\!\\!\mathbb{E}\left\\{\frac{\|\hat{\mathbf{x}}-\mathbf{x}\|^{2}}{\|\mathbf{x}\|^{2}}\right\\}$
with $\hat{x}$ being the estimate of $x$ and the $p$-th element of $x$
consisting of $\\{g_{p},\phi_{p},\psi_{p},\tau_{p}\\}$. Note that
$\mathbf{g}\\!\\!=\\!\\![g_{0},\\!\cdots\\!,g_{P-1}]^{T}$,
$\boldsymbol{\phi}\\!\\!=\\!\\![\phi_{0},\\!\cdots\\!,\phi_{P-1}]^{T}$,
$\boldsymbol{\psi}\\!\\!=\\!\\![\psi_{0},\\!\cdots\\!,\psi_{P-1}]^{T}$,
$\boldsymbol{\tau}\\!\\!=\\!\\![\tau_{0},\\!\cdots\\!,\tau_{P-1}]^{T}$, and
the NMSE of ‘Angle $\phi/\psi$’ in Fig. 2 and Fig. 3 denotes the average of
the angle $\boldsymbol{\phi}$’s NMSE and $\boldsymbol{\psi}$’s NMSE.

Fig. 2 plots the NMSEs of channel parameters estimates with respect to SNR,
where the number of pilot symbols is $M=5$, $10$, $15$, respectively. Notice
that the number of pilot carriers and elements at RIS are fixed to $L=12$ and
$N_{\text{r}}=16\times 16$, respectively. It can be seen from Fig. 2 that all
NMSEs decrease as the SNR increases, and the proposed scheme achieves a good
estimation performance even at low SNR, which demonstrates its effectiveness.
Moreover, for the same SNR, the NMSE decreases with the increase of $M$, which
shows that more pilot symbols can enhance the parameter recovery performance.

Figure 2: NMSEs of channel parameters versus SNR at different $M$.

Fig. 3 plots the NMSEs of channel parameters estimates versus SNR for
different $L$ and $N_{\text{r}}$, where $M=16$. It can be seen from Fig. 3
that all NMSEs of the parameters decrease as the SNR increases. At the same
SNR, the performance of parameter estimation becomes better as the number of
RIS elements or pilot subcarriers increases as expected.

Figure 3: NMSEs of channel parameters versus SNR at different $L$ and
$N_{\text{r}}$.

In Fig. 4, we employ the NMSE of estimated cascaded channels, i.e.,
$\text{NMSE}_{\widetilde{\mathbf{H}}}\\!=\\!\mathbb{E}\Big{\\{}\\!\frac{\|\widehat{\widetilde{\mathbf{H}}}\\!-\\!\widetilde{\mathbf{H}}\|^{2}}{\|\widetilde{\mathbf{H}}\|^{2}}\\!\Big{\\}}$,
as performance indicator with $\widehat{\widetilde{H}}_{k,m,b}$ as the
estimate of $\widetilde{H}_{k,m,b}$, where
$\widetilde{\mathbf{H}}\\!=\\![\widetilde{H}_{0,0,0},\cdots,\widetilde{H}_{L\\!-\\!1,M\\!-\\!1,N_{\text{b}}\\!-\\!1}]^{T}$.
Fig. 4 compares the cascaded channel NMSEs of the OMP method [14], the NOMP
method and CRLB at different SNRs, where $M\\!=\\!16$, $L\\!=\\!12$ and
$N_{\text{r}}\\!=\\!16\\!\times\\!16$. It can be seen from Fig. 4 that the
NMSE of the proposed NOMP method is always better than that of the OMP method,
and the gap between the NOMP and OMP methods increases with the SNR increase.
Moreover, it can be checked that the performance of the OMP method is limited
by the sampling rate of searching and the extended Newton method can further
enhance the estimation performance.

Figure 4: Comparison of cascaded channels’ NMSEs of OMP, NOMP and CRLB versus
SNR.

##  V Conclusion

In this letter, we have presented a wideband channel estimation scheme for an
RIS assisted mmWave MIMO system. We have firstly described the two individual
channels and the frequency response of the received signal, where the wideband
effect is considered. Then, we have expressed channel estimation as a
parameter recovery problem and utilized the NOMP algorithm to estimate the
channel with a few OFDM pilot symbols. Moreover, we have derived the CRLB for
channel estimates. Simulation results have verified that the proposed channel
estimation scheme achieves promising NMSE performances, especially at low SNR
and with a few pilot symbols.

From (41), we can derive
$\boldsymbol{\mathcal{F}}(\boldsymbol{\xi}_{p},\boldsymbol{\xi}_{p})$ as

$\displaystyle\boldsymbol{\mathcal{F}}(\boldsymbol{\xi}_{p},\boldsymbol{\xi}_{p})=$
$\displaystyle\frac{2}{\sigma_{\text{v}}^{2}}\\!\\!\left[\begin{array}[]{ccccc}\\!\\!\\!\\!1\\!\\!\\!&\\!\\!\\!|g_{p}|\\!\\!\\!&\\!\\!\\!|g_{p}|\\!\\!\\!&\\!\\!\\!|g_{p}|\\!\\!\\!&\\!\\!\\!2\pi\Delta
f|g_{p}|\\!\\!\\!\\!\\\
\\!\\!\\!\\!|g_{p}|\\!\\!\\!&\\!\\!\\!|g_{p}|^{2}\\!\\!\\!&\\!\\!\\!|g_{p}|^{2}\\!\\!\\!&\\!\\!\\!|g_{p}|^{2}\\!\\!\\!&\\!\\!\\!2\pi\Delta
f|g_{p}|^{2}\\!\\!\\!\\!\\\
\\!\\!\\!\\!|g_{p}|\\!\\!\\!&\\!\\!\\!|g_{p}|^{2}\\!\\!\\!&\\!\\!\\!|g_{p}|^{2}\\!\\!\\!&\\!\\!\\!|g_{p}|^{2}\\!\\!\\!&\\!\\!\\!2\pi\Delta
f|g_{p}|^{2}\\!\\!\\!\\!\\\
\\!\\!\\!\\!|g_{p}|\\!\\!\\!&\\!\\!\\!|g_{p}|^{2}\\!\\!\\!&\\!\\!\\!|g_{p}|^{2}\\!\\!\\!&\\!\\!\\!|g_{p}|^{2}\\!\\!\\!&\\!\\!\\!2\pi\Delta
f|g_{p}|^{2}\\!\\!\\!\\!\\\ \\!\\!\\!\\!2\pi\Delta
f|g_{p}|\\!\\!\\!&\\!\\!\\!2\pi\Delta
f|g_{p}|^{2}\\!\\!\\!&\\!\\!\\!2\pi\Delta
f|g_{p}|^{2}\\!\\!\\!&\\!\\!\\!2\pi\Delta
f|g_{p}|^{2}\\!\\!\\!&\\!\\!\\!(2\pi\Delta
f|g_{p}|)^{2}\\!\\!\\!\\!\end{array}\right]$ (50)
$\displaystyle\odot\\!\\!\\!\left[\begin{array}[]{ccccc}\\!\\!\\!\\!\\!\Re\\{[\mathbf{B}_{p}]\\!_{1\\!,\\!1}\\}\\!\\!\\!&\\!\\!\\!-\Im\\{[\mathbf{B}_{p}]\\!_{1\\!,\\!1}\\}\\!\\!\\!&\\!\\!\\!\Re\\{[\mathbf{B}_{p}]\\!_{1\\!,\\!2}\\}\\!\\!\\!&\\!\\!\\!\Re\\{[\mathbf{B}_{p}]\\!_{1\\!,\\!3}\\}\\!\\!\\!&\\!\\!\\!\Im\\{[\widetilde{\mathbf{B}}_{p}]\\!_{1\\!,\\!1}\\}\\!\\!\\!\\!\\!\\!\\\
\\!\\!\\!\\!\\!\Im\\{[\mathbf{B}_{p}]\\!_{1\\!,\\!1}\\}\\!\\!\\!\\!&\\!\\!\\!\\!\Re\\{[\mathbf{B}_{p}]\\!_{1\\!,\\!1}\\}\\!\\!\\!\\!&\\!\\!\\!\\!\Im\\{[\mathbf{B}_{p}]\\!_{1\\!,\\!2}\\}\\!\\!\\!\\!&\\!\\!\\!\\!\Im\\{[\mathbf{B}_{p}]\\!_{1\\!,\\!3}\\}\\!\\!\\!\\!&\\!\\!\\!\\!-\Re\\{[\widetilde{\mathbf{B}}_{p}]\\!_{1\\!,\\!1}\\}\\!\\!\\!\\!\\!\\!\\\
\\!\\!\\!\\!\\!\Re\\{[\mathbf{B}_{p}]\\!_{2\\!,\\!1}\\}\\!\\!\\!\\!&\\!\\!\\!\\!-\Im\\{[\mathbf{B}_{p}]\\!_{2\\!,\\!1}\\}\\!\\!\\!\\!&\\!\\!\\!\\!\Re\\{[\mathbf{B}_{p}]\\!_{2\\!,\\!2}\\}\\!\\!\\!\\!&\\!\\!\\!\\!\Re\\{[\mathbf{B}_{p}]\\!_{2\\!,\\!3}\\}\\!\\!\\!\\!&\\!\\!\\!\\!\Im\\{[\widetilde{\mathbf{B}}_{p}]\\!_{2\\!,\\!1}\\}\\!\\!\\!\\!\\!\\!\\\
\\!\\!\\!\\!\\!\Re\\{[\mathbf{B}_{p}]\\!_{3\\!,\\!1}\\}\\!\\!\\!\\!&\\!\\!\\!\\!-\Im\\{[\mathbf{B}_{p}]\\!_{3\\!,\\!1}\\}\\!\\!\\!\\!&\\!\\!\\!\\!\Re\\{[\mathbf{B}_{p}]\\!_{3\\!,\\!2}\\}\\!\\!\\!\\!&\\!\\!\\!\\!\Re\\{[\mathbf{B}_{p}]\\!_{3\\!,\\!3}\\}\\!\\!\\!\\!&\\!\\!\\!\\!\Im\\{[\widetilde{\mathbf{B}}_{p}]\\!_{3\\!,\\!1}\\}\\!\\!\\!\\!\\!\\!\\\
\\!\\!\\!\\!\\!-\Re\\{[\widetilde{\mathbf{B}}_{p}]\\!_{1\\!,\\!1}\\}\\!\\!\\!&\\!\\!\\!-\Re\\{[\widetilde{\mathbf{B}}_{p}]\\!_{1\\!,\\!1}\\}\\!\\!\\!&\\!\\!\\!-\Im\\{[\widetilde{\mathbf{B}}_{p}]\\!_{1\\!,\\!2}\\}\\!\\!\\!&\\!\\!\\!-\Im\\{[\widetilde{\mathbf{B}}_{p}]\\!_{1\\!,\\!3}\\}\\!\\!\\!&\\!\\!\\!\Re\\{C_{p}\\}\\!\\!\\!\\!\\!\\!\end{array}\right]$
(56)

where $\breve{\mathbf{a}}_{k,m,b}\triangleq[\mathbf{a}_{\text{B}}(k\Delta
f,\theta_{\text{b}})]_{b}\mathbf{a}^{T}_{\text{R}}(k\Delta
f,\phi_{\text{r}},\psi_{\text{r}})\boldsymbol{\Omega}_{m}$,
$\breve{A}_{k,m,b,p}\triangleq\breve{\mathbf{a}}_{k,m,b}\mathbf{a}_{\text{R}}(k\Delta
f,\phi_{p},\psi_{p})$,
$\mathbf{B}_{p}\triangleq\sum_{k,m,b}\mathbf{D}_{k,m,b,p}$,
$\widetilde{\mathbf{B}}_{p}\triangleq\sum_{k,m,b}k\mathbf{D}_{k,m,b,p}$,
$C_{p}\triangleq\sum_{k,m,b}k^{2}\breve{A}_{k,m,b,p}^{*}\breve{A}_{k,m,b,p}$,

$\displaystyle\mathbf{D}_{k,m,b,p}\\!\\!\triangleq\\!\\!\left[\begin{array}[]{ccc}\\!\\!\breve{A}_{k,m,b,p}^{*}\\!\\!\\\
\\!\\!\frac{\partial\breve{A}_{k,m,b,p}^{*}}{\partial\phi_{p}}\\!\\!\\\
\\!\\!\frac{\partial\breve{A}_{k,m,b,p}^{*}}{\partial\psi_{p}}\\!\\!\end{array}\right]\\!\\!\left[\begin{array}[]{ccc}\\!\\!\breve{A}_{k,m,b,p}\\!\\!&\\!\\!\frac{\partial\breve{A}_{k,m,b,p}}{\partial\phi_{p}}\\!\\!&\\!\\!\frac{\partial\breve{A}_{k,m,b,p}}{\partial\psi_{p}}\end{array}\\!\\!\right],$
(61)
$\displaystyle\left[\begin{array}[]{cccc}\\!\\!\\!\\!\frac{\partial\breve{A}_{k\\!,m\\!,b\\!,p}}{\partial\phi_{p}}\\!\\!\\!\\!\\!\\\
\\!\\!\\!\\!\frac{\partial\breve{A}_{k\\!,m\\!,b\\!,p}}{\partial\psi_{p}}\\!\\!\\!\\!\\!\\\
\\!\\!\\!\\!\frac{\partial\breve{A}^{*}_{k\\!,m\\!,b\\!,p}}{\partial\phi_{p}}\\!\\!\\!\\!\\!\\\
\\!\\!\\!\\!\frac{\partial\breve{A}^{*}_{k\\!,m\\!,b\\!,p}}{\partial\psi_{p}}\\!\\!\\!\\!\\!\end{array}\right]\\!\\!\\!=\\!\\!\\!\left[\begin{array}[]{cccc}\\!\\!\\!\\!-\\!\jmath
2\pi(1\\!\\!+\\!\\!\frac{k\Delta
f}{f_{c}})\breve{\mathbf{a}}_{k\\!,m\\!,b}\big{(}{\mathbf{a}}_{\text{R}}(k\Delta
f\\!,\phi\\!_{p}\\!,\psi\\!_{p})\\!\odot\\!\dot{\boldsymbol{\varpi}}\\!_{\phi\\!_{p}}\big{)}\\!\\!\\!\\!\\!\\\
\\!\\!\\!\\!-\\!\jmath 2\pi(1\\!\\!+\\!\\!\frac{k\Delta
f}{f_{c}})\breve{\mathbf{a}}_{k\\!,m\\!,b}\big{(}{\mathbf{a}}_{\text{R}}(k\Delta
f\\!,\phi\\!_{p}\\!,\psi\\!_{p})\\!\odot\\!\dot{\boldsymbol{\varpi}}\\!_{\psi\\!_{p}}\big{)}\\!\\!\\!\\!\\!\\\
\\!\\!\\!\\!\jmath 2\pi(1\\!\\!+\\!\\!\frac{k\Delta
f}{f_{c}})\breve{\mathbf{a}}_{k\\!,m\\!,b}^{*}\big{(}{\mathbf{a}}^{*}_{\text{R}}(k\Delta
f\\!,\phi\\!_{p}\\!,\psi\\!_{p})\\!\odot\\!\dot{\boldsymbol{\varpi}}\\!_{\phi\\!_{p}}\big{)}\\!\\!\\!\\!\\!\\\
\\!\\!\\!\\!\jmath 2\pi(1\\!\\!+\\!\\!\frac{k\Delta
f}{f_{c}})\breve{\mathbf{a}}_{k\\!,m\\!,b}^{*}\big{(}{\mathbf{a}}^{*}_{\text{R}}(k\Delta
f\\!,\phi\\!_{p}\\!,\psi\\!_{p})\\!\odot\\!\dot{\boldsymbol{\varpi}}\\!_{\psi\\!_{p}}\big{)}\\!\\!\\!\\!\\!\end{array}\right],$
(70)
$\displaystyle\dot{\boldsymbol{\varpi}}\\!_{\phi\\!_{p}}\\!\\!\\!=\\!\\!\Big{[}\\!0,\\!\cdots\\!,\\!\frac{d\\!\big{(}\\!(N\\!_{\text{x}}\\!\\!-\\!\\!1)\\!\cos(\\!\phi\\!_{p}\\!)\\!\sin(\\!\psi\\!_{p}\\!)\\!\\!+\\!\\!(N\\!_{\text{y}}\\!\\!-\\!\\!1)\\!\cos(\\!\phi\\!_{p}\\!)\\!\cos(\\!\psi\\!_{p}\\!)\\!\big{)}}{\lambda_{\text{c}}}\\!\Big{]}^{T},$
(71)
$\displaystyle\dot{\boldsymbol{\varpi}}\\!_{\psi\\!_{p}}\\!\\!\\!=\\!\\!\Big{[}\\!0,\\!\cdots\\!,\\!\frac{d\big{(}\\!(N\\!_{\text{x}}\\!\\!-\\!\\!1)\\!\sin(\\!\phi\\!_{p}\\!)\\!\cos(\\!\psi\\!_{p}\\!)\\!\\!-\\!\\!(N\\!_{\text{y}}\\!\\!-\\!\\!1)\\!\sin(\\!\phi\\!_{p}\\!)\\!\sin(\\!\psi\\!_{p}\\!)\\!\big{)}}{\lambda_{\text{c}}}\\!\Big{]}^{T}.$
(72)

## References

  * [1] M. Nemati et. al., “RIS-assisted coverage enhancement in millimeter-wave cellular networks,” _IEEE Access_ , vol. 8, pp. 188171-188185, 2020.
  * [2] J. Ma et. al., “Sparse Bayesian learning for the time-varying massive MIMO channels: Acquisition and tracking,” _IEEE Trans. Commun._ , vol. 67, no. 3, pp. 1925-1938, Mar. 2019.
  * [3] E. Basar et. al., “Present and future of reconfigurable intelligent surface-empowered communications,” _arXiv:2105.00671_ , 2021. [Online]. Available: https://arxiv.org/abs/2105.00671.
  * [4] M. Nemati et. al., “Modeling RIS empowered outdoor-to-indoor communication in mmWave cellular networks,” _arXiv:2101.00736_ , 2021. [Online]. Available: https://arxiv.org/abs/2101.00736.
  * [5] H. Liu et. al., “Matrix-calibration-based cascaded channel estimation for reconfigurable intelligent surface assisted multiuser MIMO,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 11, pp. 2621-2636, Nov. 2020.
  * [6] J. He et. al., “Channel estimation for RIS-aided mmWave MIMO channels,” _arXiv:2002.06453_ , 2020. [Online]. Available: https://arxiv.org/abs/2002.06453.
  * [7] X. Hu et. al., “Location information aided multiple intelligent reflecting surface systems,” _IEEE Trans. Commun._ , vol. 68, no. 12, pp. 7948-7962, Dec. 2020.
  * [8] B. Wang et. al., “Beam squint and channel estimation for wideband mmWave massive MIMO-OFDM systems,” _IEEE Trans. Signal Process._ , vol. 67, no. 23, pp. 5893-5908, Dec. 2019.
  * [9] S. Zhang et. al., “Deep learning optimized sparse antenna activation for reconfigurable intelligent surface assisted communication,” _arXiv:2009.01607_ , 2020. [Online]. Available: https://arxiv.org/abs/2009.01607.
  * [10] B. Deepak et. al., “Channel estimation in reconfigurable intelligent surface assisted mmWave MIMO systems,” _arXiv:2011.00900_ , 2020. [Online]. Available: https://arxiv.org/abs/2011.00900.
  * [11] Y. Han et. al., “Efficient downlink channel reconstruction for FDD multi-antenna systems,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 6, pp. 3161-3176, Jun. 2019.
  * [12] W. -T. Shih et. al., “Fast antenna and beam switching method for mmWave handsets with hand blockage,” _arXiv:2103.08151_ , 2021. [Online]. Available: https://arxiv.org/abs/2103.08151.
  * [13] D. Fan et. al., “Angle domain channel estimation in hybrid millimeter wave massive MIMO systems,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 12, pp. 8165-8179, Dec. 2018.
  * [14] J. Lee et. al., “Channel estimation via orthogonal matching pursuit for hybrid MIMO systems in millimeter wave communications,” _IEEE Trans. Commun._ , vol. 64, no. 6, pp. 2370-2386, Jun. 2016.