# Optimizations and Heuristics to improve Compression in Columnar Database
Systems

Jayanth

Stanford University, CA, US<EMAIL_ADDRESS>
###### Abstract

In-memory columnar databases have become mainstream over the last decade and
have vastly improved the fast processing of large volumes of data through
multi-core parallelism and in-memory compression thereby eliminating the usual
bottlenecks associated with disk-based databases. For scenarios, where the
data volume grows into terabytes and petabytes, keeping all the data in memory
is exorbitantly expensive. Hence, the data is compressed efficiently using
different algorithms to exploit the multi-core parallelization technologies
for query processing. Several compression methods are studied for compressing
the column array, post Dictionary Encoding. In this paper, we will present two
novel optimizations in compression techniques - Block Size Optimized Cluster
Encoding and Block Size Optimized Indirect Encoding - which perform better
than their predecessors. In the end, we also propose heuristics to choose the
best encoding amongst common compression schemes.

###### keywords:

Database Compression; Database Optimization; Encoding Schemes; In-memory
Columnar Databases

##  1 Introduction

Traditional database systems have relied on compression [2] mechanisms to
improve performance significantly. Compression not just reduces the size of
the data but also improves I/O performance. Columnar databases have gained
momentum over the last decade [1]. In this paper, we discuss various state-of-
art compression schemes in columnar database systems. We, further, propose
optimizations and heuristics to arrive at the right compression scheme. A
columnar database/column-store is one in which each attribute is stored in a
separate column such that successive values are stored consecutively on the
disk. By virtue of storing successive values close, such databases present
significant lossless compression ability in the form of dictionary encoding
[3][4] .

For any given column, distinct values appearing in the column are stored in a
dictionary as an array with keys being the indexes of the array. The values
comprising the dictionary are sorted. Each value in the sorted dictionary is
identified by its value ID which is given by its position in the value array.
The actual data in the column is stored as an array of value IDs. The sorted
dictionary and value ID array together represent the actual column.

Dictionary compression presents the advantage of executing search on the
integer values. Since the dictionary is sorted by values, comparison and range
selection queries can simply be performed on the value IDs. This yields better
performance as opposed to when the actual values stored in the columns are
strings. For example, the query select c from T where name > James can
internally be replaced by finding all entries for which the value ID is
greater than the value ID of the greatest value in the dictionary that is less
or equal to James. Furthermore in the context of in-memory columnar stores,
dictionary compressed data can be loaded faster into the CPU cache. As the
bottleneck is often data transport between memory and cache, performance gain
exceeds the computing time needed for decompression.

##  2 Compression Schemes

Given dictionary encoding wherein a column gets encoded as a dictionary and
value ID array, the size of the value ID array readily grows with large data
volumes spanning billions of records. For in-memory systems, it is critical to
further compress the value ID array using advanced compression schemes. These
compression mechanisms are explained further [4].

Figure 1: Compression Schemes [3]

###  2.1 Prefix Encoding

Real-world databases often comprise of columns that contain predominantly one
value and remaining values are mostly unique. In order to apply prefix
encoding, the value ID array is sorted such that this predominant value
appears at the start of the array. The sequence is then replaced by storing
the value once, along with the number of occurrences.

###  2.2 Run Length Encoding

For a sorted value ID array, this encoding replaces the value sequences with a
single instance of the value followed by the number of occurrences.

###  2.3 Sparse Encoding

For value ID array containing predominantly a single value V, sparse encoding
removes V and introduces a bit vector to specify all the positions at which V
occurs in the original array.

###  2.4 Cluster Encoding

Cluster encoding divides the value ID array into N fixed size blocks/clusters
(typically 1024 elements). If a cluster contains unique value, it is replaced
by the single instance of that value else it remains uncompressed.
Additionally, a bit vector of length N indicates the clusters replaced by
single value (1 if compressed else 0).

###  2.5 Indirect Encoding

Similar to cluster encoding, indirect encoding also addresses value ID array,
divided into fixed size blocks. For each block containing few distinct values,
an additional dictionary is used to encode the values in that block. There is
a level of indirections introduced by local dictionaries for blocks in
addition to a global dictionary.

##  3 Optimisations

While these advanced compression schemes prove useful, there are some basic
flaws in their mechanisms. Both cluster and indirect encoding consider fixed
size blocks. This does not provide optimal compression because it is highly
unlikely that all elements in the cluster have the same value. Particularly,
for in-memory systems it gets critical to fine tune the size of the blocks.
The optimizations, discussed further, derive optimal block size based on the
value ID array instead of using a fixed block size b ( 1024 elements,
currently). This won’t affect the query processor as instead of a fixed value
1024, it will use the block size specific to that column stored along with the
value id array.

###  3.1 Block Size Optimized Cluster Encoding

Let N be the total size of the value ID array and b be the block size. “size”
here represents number of elements constituting the entity. Without any
compression, the size of the value ID array is given by N * size(valueID).
Clustering encoding reduces the size to :

$(N-(S*b))*size(valueID)+S*size(valueID)+log(\frac{N}{b})$ (1)

where $S$ is the number of clusters, which contain only occurrences of a
single value. This equation can be re-written as:

$(N-S*(b-1))*size(valueID)$ (2)

considering log terms of lower complexity than rest, i.e. size of bit-vector
is negligible. Thus, the optimal block size $b^{*}$ can be given
mathematically as:

$b^{*}=argmin_{b\in\cup_{i=1}^{log(N)}\\{2^{i}\\}}(N-S*(b-1))$ (3)

Since, N is fixed for a dataset i.e. the count of records, it can removed to
rewrite the optimization formultation as:

$b^{*}=argmax_{b\in\cup_{i=1}^{log(N)}\\{2^{i}\\}}S*(b-1)$ (4)

The above optimization algorithm can be solved iteratively in complexity of
$O(Nlog(N))$. Further, optimization can be done to look for b till $\sqrt{N}$
and also, finding the value of S can be done dynamically from Run Length
Encoding of the value ID array. Let, the value ID array is represented as
$\cup_{i}\\{(v_{i},c_{i})\\}$, where $v_{i}$ is the ith value ID and $c_{i}$
is its count. Let, $S_{i}$ denote the clustered blocks formed from the start
of $(i-1)th$ value ID $v_{i-1}$ till the $ith$ value ID $v_{i}$ and $r_{i}$
denote the number of value IDs unclustered till start of $ith$ value ID. Thus,
$S_{i}$ and $r_{i}$ can be written as following recurrence relations, which
can be exploited to solve using dynamic programming:

$\displaystyle S_{1}$ $\displaystyle=0$ (5) $\displaystyle r_{1}$
$\displaystyle=0$ (6) $\displaystyle S_{i}$
$\displaystyle=\frac{c_{i-1}-(|b-r_{i-1}|\mod b)}{b}$ (7) $\displaystyle
r_{i}$ $\displaystyle=(\sum_{j=1}^{i-1}c_{j})\mod b$ (8)

Finally, the total number of clustered blocks $S$, required to solve the
optimization problem can be derived as:

$S=\sum_{i=1}^{..}S_{i}$ (9)

⬇

def findClusteredBlocks(valCountArray,blockSize):

S_i=[];R_i=[];S=0;R_i.append(0)

for i in range(0,valCountArray.__len__()):

R_i.append((R_i[i] + valCountArray[i][1]) % blockSize)

S_i.append(max((valCountArray[i][1] - ((blockSize - R_i[i]) %
blockSize))/blockSize, 0))

S +=S_i[i]

return S

def findOptimalBlockSize(valArray):

valCountArray = convertValCount(valArray)

Fmax = 0; OptimalBlockSize = 1

for i in range(1, int(numpy.log2(valArray.__len__()))):

b = 2 ** i

S = findClusteredBlocks(valCountArray, b)

F = S * (b - 1)

print i, b, S, F

if F > Fmax:

Fmax = F; OptimalBlockSize = b

return OptimalBlockSize

###  3.2 Block Size Optimized Indirect Encoding

In Indirect encoding, if a block contains only a few distinct values, an
additional dictionary is used to encode the values in that block. The
implementation also needs to store the information which blocks are encoded
with an additional dictionary and the links to the additional dictionaries. If
there are very many distinct values in a column, compression of the value ID
array is not used. Indirect encoding can be applied efficiently, if data
blocks hold a few distinct values as Indirect encoding saves space by mapping
large global value IDs to small local value IDs.

The few distinct values criteria points towards usage of shanon entropy
function, which has minima if all values are same while maxima, when all
values are distinct and thus, have uniform distribution. However, each block
should have as few distinct values and thus, on average the entropy of blocks
should be minimized. Thus, the requirement for the blocks to have as few
distinct values as possible can be used to formulate the following
optimization problem for finding the optimal block size $b^{*}$:

$b^{*}=argmin_{b\in\cup_{i=1}^{log(N)}\\{2^{i}\\}}E[H_{b}(^{b}X_{i})]$ (10)

where $E$ is the expectation or average of the b-ary entropy $H_{b}(X_{i})$ of
the blocks ${}^{b}X_{i}$ defined by block size $b$. The reason for using b-ary
entropy is to distinguish the entropy for differnt block sizes as the maximum
different values the block can hold is $b$. Here, the b-ary entropy
$H_{b}(X_{i})$ will be computed as:

$H_{b}(X_{i})=E[-log_{b}(P(X_{i})]$ (11)

The above optimization problem can be solved using dynamic programming in
complexity of $O(Nlog(N))$ by merging the frequency distribution at previous
step into subsequent step.

⬇

def findOptimalBlockSizeEntropy(valArray):

Fmin = float("inf"); OptimalBlockSize = 1

for i in range(1, int(numpy.log2(valArray.__len__()))):

b = 2 ** i; data = []; S = 0

for i in range(valArray.__len__()):

data.append(valArray[i])

if (i + 1) % b == 0:

S += H(data, b)

data = []

F = S / math.ceil(valArray.__len__() / float(b))

if F < Fmin:

Fmin = F

OptimalBlockSize = b

return OptimalBlockSize

##  4 Encoding decision heuristics

The choice between different encoding decisions is important, as wrong
encoding sometimes may even lead to increase in data size, and thus, bad query
performance. Before proposing any algorithm or decision tree to decide over
encoding decision, lets distinguish them on some important characteristics:

Table 1: Heuristics Criteria Encoding | Blocks | Sorting | Sparsity | Prefix
---|---|---|---|---
Prefix | NO | * | YES | NO
Run-Length | NO | YES | NO | NO
Cluster | YES | * | NO | YES
Sparse | NO | NO | YES | YES
Indirect | YES | NO | NO | YES

Here, the columns in above table refer to the following:

  1. 1.

Blocks \- This column refers to the requirement of chunking the valueID array
into blocks and thus, maximum block repetitions provide optimal compression.
This holds as per the encoding algorithm only valid for Cluster and Indirect
encoding. However, both encoding had an inherent limitation of using fixed
block size as 1024 elements, which has been removed in section 3.

  2. 2.

Sorting \- This column refers to the requirement of sorting the valueID array
for good compression. Prefix encoding requires sorting the table in a way that
puts the most frequent value at the beginning of the column. Cluster encoding,
benefits from sorting but then, Run Length encoding is more optimal encoding
choice.

  3. 3.

Sparsity \- This column refers to the requirement of sparsity in the valueID
array, i.e. repetition of a single value. Apart from Sparse encoding, Prefix
encoding also benefits from high sparsity, if the sparse value ID is the
prefix itself.

  4. 4.

Prefix Encoding \- This column refers to the possibility and benefit of
applying prefix encoding before the encoding choice for more compression.
Cluster encoding, indirect encoding and sparse encoding can be combined with
prefix encoding, after sorting the frequent value ID as prefix.

Based on these observations, a sample algorithmic flowchart for encoding
decision can be proposed as figure 2.

Figure 2: Flowchart for encoding decision on a column

The first condition checks whether the column is candidate key, i.e. all
values are distinct. If all values are distinct, none of the above
compressions will work. If all values are distinct, apply no compression. In
special case, when the column is sequentially increasing or decreasing as in
Time-Series table, linear Run Length Encoding should be applied. The next
condition checks whether the column is sorted, i.e. all the value IDs are in
order. If the value ID array is sorted, then Run Length Encoding should be
applied as it gives the optimal compression for sorted value id array. Then,
the heuristics finds the frequency of each distinct value ID in the value ID
array using a group by distinct count query. If the count of 1st value ID in
value ID array is greater than 2 and that values repeats itself in the
beginning, then Prefix Encoding should be applied. The following condition
checks whether the sparsity in the value ID array is greater than x, a tune-
able parameter (x>1). The sparsity is quantified as the ratio of the maximum
of the frequency of value ID array and the average of the frequencies of value
ID array i.e. how much more is the most sparse value ID repeating over average
repetition. If the sparsity in the value ID array is significantly large, then
Sparse Encoding should be applied. The next condition checks whether there is
sufficient repetition in the value ID array i.e. the average repetition is
greater than y, a tune-able parameter (y>1). The average repetition is
quantified as the average of the frequencies of value ID array. The final
condition checks whether the clustered value IDs after cluster encoding lead
to high compression. The ratio of the size of cluster of compressed value IDs
using cluster encoding and the total value ID array size is compared with z, a
tune-able parameter (0<z<1). Here, S is the number of the clustered blocks and
b is the optimal block size, both values calculated as in in section 3 after
applying the algorithm for calculating optimal block size in cluster encoding.
If the clustered value IDs after cluster encoding lead to high compression,
then Cluster Encoding should be applied. If the condition to apply cluster
encoding is not met but there is sufficient repetition, then Indirect Encoding
should be applied.

##  5 Conclusions

Optimal block size in cluster and indirect encoding perform better than the
conventional techniques which uses fixed block size. Furthermore, the
heuristics based on certain tunable parameters can hint towards the right
compression scheme. While dictionary encoding ensures good compression,
optimizations and heuristics to choose the right advanced compression
technique can further reduce the storage footprint.

## References

  * [1] D. Abadi, S. Madden, and M. Ferreira.  Integrating compression and execution in column-oriented database systems.  In Proceedings of the 2006 ACM SIGMOD international conference on Management of data, pages 671–682. ACM, 2006.
  * [2] G. V. Cormack.  Data compression on a database system.  Communications of the ACM, 28(12):1336–1342, 1985.
  * [3] H. Plattner.  A course in in-memory data management.  Springer, 2013.
  * [4] H. Plattner and A. Zeier.  In-memory data management: technology and applications.  Springer Science & Business Media, 2012.


# Majorana corner states in square and Kagome quantum spin liquids

Haoran Wang  Department of Physics and Astronomy, University of Manchester,
Manchester M13 9PL, UK  Alessandro Principi
<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Manchester, Manchester M13 9PL, UK

###### Abstract

Quantum spin liquids hosting Majorana excitations have recently experienced
renewed interest for potential applications to topological quantum
computation. Performing logical operations with reduced poisoning requires to
localize such quasiparticles at specific point of the system, with energies
that are well defined and inside the bulk energy gap. These are two defining
features of second order topological insulators (SOTIs). Here, we show two
spin models that support quantum spin liquid phases characterised by Majorana
excitations and that behave as SOTIs, one of which is analytically solvable
thanks to a theorem by Lieb. We show that, depending on the values of spin
couplings, it is possible to localize either fermions or Majorana particles at
their corners.

Introduction—Quantum spin liquids (QSLs) are intriguing states of matter in
which spins never freeze due to their high degree of entanglement Anderson
(1973); Zhou  _et al._ (2017); Savary and Balents (2016); Wen (2004).
Recently, they have experienced rekindled interest because of two main
factors. On the one hand, the development of exactly-solvable spin-lattice
models, for example Kitaev’s Kitaev (2006), which host QSL phases whose
excitations are neither fermions nor bosons but Majorana particles. On the
other hand, the discovery of new materials compatible with such models
Banerjee  _et al._ (2017).

Majorana particles are one of the cornerstones of research in topological
quantum computation, since they can be used to design quantum gates that are
resilient to external noise Kitaev (2003); Sarma  _et al._ (2015); Hoffman
_et al._ (2016); Lian  _et al._ (2018). In order to use them to perform
logical operations, it is however crucial to find ways to localize them at
specific points of the system, with well-defined energies inside an energy
gap. The former property is missing in the original Kitaev QSL model. There,
1D Majorana channels can be formed at the edges, in the presence of a magnetic
field Kitaev (2006). (In principle, it is also possible to localize Majorana
particles at vortices of the $\mathbb{Z}_{2}$ gauge field Kitaev (2006);
Knolle (2016). However, generating and controlling such nanoscale objects in
experiments would be extremely challenging.) In this paper we show two spin
models, one of which analytically solvable, which support topologically-
protected Majorana corner states with energies in the middle of the bulk band
gap.

Mid-gap corner states emerge in two-dimensional (2D) “second-order”
topological insulators (SOTIs) Benalcazar  _et al._ (2017a, b); Schindler  _et
al._ (2018a); Langbehn  _et al._ (2017); Geier  _et al._ (2018); Ezawa
(2018a); Song  _et al._ (2017); Ezawa (2018b, c). $d$-dimensional SOTIs are
insulating both in the bulk and at the surfaces, but exhibit gapless
$d-2$-states protected by a variety of crystalline symmetries Benalcazar  _et
al._ (2017b) (among others, mirror reflection, twofold rotation, or inversion
symmetry). SOTIs have been realized in various experiments, most notably in
artificial settings such as mechanical Serra-Garcia  _et al._ (2018),
electrical Imhof  _et al._ (2018), microwave Peterson  _et al._ (2018), and
photonic El Hassan  _et al._ (2019) devices, but they have also been shown to
occur in Nature Schindler  _et al._ (2018b). A similar phenomenology has been
predicted to occur in non-Hermitian systems Luo and Zhang (2019), topological
superconductors Liu  _et al._ (2018); Zhu (2018); Laubscher  _et al._ (2019);
Wang  _et al._ (2018); Kheirkhah  _et al._ (2020); Yan (2019) and QSLs Dwivedi
_et al._ (2018). In the latter two, corner states can be Majorana particles.
However, to the best of our knowledge, no analytically solvable QSL has been
shown to exhibit Majorana corner states 111The model of Ref. Dwivedi  _et al._
(2018) can be solved exactly by mapping spins into Majorana particles which
are free to propagate on top a quenched $\mathbb{Z}_{2}$ gauge potential.
However, the ground-state flux configuration is not analytically known and
must be found with numerical techniques..

Here we study two frustrated spin-$3/2$ systems which exhibit QSL phases and
corner states at low temperature. The models we study are similar to that of
Ref. Dwivedi  _et al._ (2018) but are defined on Kagome Chua  _et al._ (2011)
and square Yao  _et al._ (2009) lattices, rather than on the Shastry-
Sutherland one. Following Kitaev’s construction Kitaev (2006), we
fractionalize the spin degrees of freedom into two itinerant Majorana
particles and a quesched $\mathbb{Z}_{2}$ gauge potential. Notably, the ground
state of the square lattice we study is analytically known Lieb (1994). We
show that both models can support stable fermionic or Majorana corner states,
depending on the values of the spin couplings. The analytical solvability of
the square-lattice model lends credibility to the possibility of finding
Majorana corner states in QSLs.

The model—We consider two lattices, square and Kagome, at whose sites are
located spin-$3/2$ variables. The Hamiltonian of such systems can be written
in terms of $4\times 4$ matrices which operate on the four spin polarizations
at each site. The basis of $4\times 4$ matrices can be chosen to be composed
by: the identity, five Gamma matrices ${\hat{\Gamma}}^{a}$ ($a=1,\ldots,5$)
which can be represented as symmetric bilinear combinations of the spin-$3/2$
operators ${\hat{S}}^{\alpha}$ ($\alpha=x,y,z$) as Chua  _et al._ (2011); Yao
_et al._ (2009)

$\displaystyle{\hat{\Gamma}}^{1}=\frac{\\{{\hat{S}}^{y},{\hat{S}}^{z}\\}}{\sqrt{3}},~{}{\hat{\Gamma}}^{2}=\frac{\\{{\hat{S}}^{z},{\hat{S}}^{x}\\}}{\sqrt{3}},~{}{\hat{\Gamma}}^{3}=\frac{\\{{\hat{S}}^{x},{\hat{S}}^{y}\\}}{\sqrt{3}},$
$\displaystyle{\hat{\Gamma}}^{4}=\frac{1}{\sqrt{3}}\big{[}({\hat{S}}^{x})^{2}-({\hat{S}}^{y})^{2}\big{]},~{}{\hat{\Gamma}}^{5}=({\hat{S}}^{z})^{2}-\frac{5}{4},$
(1)

and the ten bilinears
${\hat{\Gamma}}^{ab}=[{\hat{\Gamma}}^{a},{\hat{\Gamma}}^{b}]/(2i)$. Here,
$[\cdot,\cdot]$ and $\\{\cdot,\cdot\\}$ are the matrix commutator and
anticommutator, respectively. In passing we note that the chosen Gamma
matrices satisfy the Clifford algebra
$\\{{\hat{\Gamma}}^{a},{\hat{\Gamma}}^{b}\\}=2\delta^{ab}$.

We now define the Gamma-matrix model Chua  _et al._ (2011); Yao  _et al._
(2009)

$\displaystyle{\hat{\cal H}}$ $\displaystyle=$ $\displaystyle
J_{1}\sum_{\langle i,j\rangle\in{\cal
P}_{1}}{\hat{\Gamma}}_{i}^{1}{\hat{\Gamma}}_{j}^{2}+J_{2}\sum_{\langle
i,j\rangle\in{\cal P}_{2}}{\hat{\Gamma}}_{i}^{3}{\hat{\Gamma}}_{j}^{4}$
$\displaystyle+$ $\displaystyle J_{1}^{\prime}\sum_{\langle i,j\rangle\in{\cal
P}_{1}}{\hat{\Gamma}}_{i}^{15}{\hat{\Gamma}}_{j}^{25}+J_{2}^{\prime}\sum_{\langle
i,j\rangle\in{\cal
P}_{2}}{\hat{\Gamma}}_{i}^{35}{\hat{\Gamma}}_{j}^{45}+J_{5}\sum_{i}{\hat{\Gamma}}_{i}^{5},$

where $i$ and $j$ label the lattice sites, while ${\cal P}_{1}$ and ${\cal
P}_{2}$ are the collections of plaquettes of types 1 and 2, respectively (see
Fig. 1 for the definition). In Eq. (Majorana corner states in square and
Kagome quantum spin liquids), the coupling between spins depends on the type
of plaquette to which the link $\langle i,j\rangle$ belongs (links are taken
in the counterclockwise direction in plaquettes 1 and 2). For the Kagome
lattice, plaquettes 1 and 2 coincide with upward and downwards triangles,
respectively. Conversely, the plaquettes of the square lattice are taken to
alternate between type 1 and 2 along every other diagonal, while all others
are of type 3. In the Kagome lattice, hexagons are type-3 plaquettes. Type-3
plaquettes do not appear explicitly in the Hamiltonian (Majorana corner states
in square and Kagome quantum spin liquids) to avoid double counting the bonds

\begin{overpic}[width=84.55576pt]{fig1a} \put(0.0,-10.0){(a)} \end{overpic} | \begin{overpic}[width=127.9171pt]{fig1b} \put(0.0,-10.0){(b)} \end{overpic}
---|---

\begin{overpic}[width=52.03558pt]{fig1c} \put(0.0,-10.0){(c)} \end{overpic} | \begin{overpic}[width=52.03558pt]{fig1d} \put(0.0,-10.0){(d)} \end{overpic} | \begin{overpic}[width=52.03558pt]{fig1e} \put(0.0,-10.0){(e)} \end{overpic} | \begin{overpic}[width=52.03558pt]{fig1f} \put(0.0,-10.0){(f)} \end{overpic}
---|---|---|---

Figure 1: A pictorial representation of the two models studied in this paper.
Panel (a) The square model. Panel (b) The Kagome model. In both cases we
identify three types of plaquettes and two types of bonds. Solid (dashed)
arrows correspond to couplings of the type
${\hat{\Gamma}}_{i}^{1}{\hat{\Gamma}}_{j}^{2}$ and
${\hat{\Gamma}}_{i}^{15}{\hat{\Gamma}}_{j}^{25}$
(${\hat{\Gamma}}_{i}^{3}{\hat{\Gamma}}_{j}^{4}$ and
${\hat{\Gamma}}_{i}^{35}{\hat{\Gamma}}_{j}^{45}$). The direction of the arrow
is from site $i$ to $j$. The shaded region represents the unit cell used in
square-lattice calculations (atoms in it are labeled with $1,2,3$ and $4$).
Panels (c)-(f) The four distinct flux patterns in the Kagome unit cell. Panels
(c) and (e) require doubling of the unit cell.

It is possible to define a set of flux operators ${\hat{W}}_{p}$, one per
plaquette $p$, which commute with the Hamiltonian and among themselves Chua
_et al._ (2011); Yao  _et al._ (2009). These are ${\hat{W}}_{p}=\prod_{i\in
p}i{\hat{\Gamma}}^{12}_{i}$, if $p$ is of type 1, ${\hat{W}}_{p}=\prod_{i\in
p}i{\hat{\Gamma}}^{34}_{i}$, if $p$ is of type 2, and
${\hat{W}}_{p}=\prod_{\langle i,j\rangle\in
p}i{\hat{\Gamma}}^{23}_{i}{\hat{\Gamma}}^{14}_{j}$, if $p$ is of type 3. All
products run over links taken in the counterclockwise direction around a given
plaquette. In the type-3 flux operator, the first link is shared with a
plaquette of type 1. The presence of such a large number of constants of
motion is at the root of the solvability of the models.

We now introduce six Majorana operators Chua  _et al._ (2011); Yao  _et al._
(2009)
${\hat{\xi}}_{i}^{1},{\hat{\xi}}_{i}^{2},{\hat{\xi}}_{i}^{3},{\hat{\xi}}_{i}^{4},{\hat{c}}_{i},$
and ${\hat{d}}_{i}$ at each site $i$, which satisfy the anticommutation
relations
$\\{{\hat{\xi}}_{i}^{a},{\hat{\xi}}_{j}^{b}\\}=2\delta_{ij}\delta^{ab}$,
$\\{{\hat{c}}_{i},{\hat{c}}_{j}\\}=2\delta_{ij}$,
$\\{{\hat{d}}_{i},{\hat{d}}_{j}\\}=2\delta_{ij}$, and
$\\{{\hat{c}}_{i},{\hat{\xi}}_{j}^{a}\\}=\\{{\hat{d}}_{i},{\hat{\xi}}_{j}^{a}\\}=0$.
The Gamma matrices are then expressed as
${\hat{\Gamma}}_{i}^{a}=i{\hat{\xi}}_{i}^{a}{\hat{c}}_{i}$,
${\hat{\Gamma}}_{i}^{a5}=i{\hat{\xi}}_{i}^{a}{\hat{d}}_{i}$ and
${\hat{\Gamma}}_{i}^{5}=i{\hat{c}}_{i}{\hat{d}}_{i}$, where $a=1,2,3,4$. This
representation enlarges the Hilbert space, introducing non-physical states
that must be projected out at the end of the calculation Chua  _et al._
(2011); Yao  _et al._ (2009); Kitaev (2006). To define the projection
operator, we note that
${\hat{D}}_{i}=-{\hat{\Gamma}}_{i}^{1}{\hat{\Gamma}}_{i}^{2}{\hat{\Gamma}}_{i}^{3}{\hat{\Gamma}}_{i}^{4}{\hat{\Gamma}}_{i}^{5}=1$.
In the Majorana representation, however, the eigenvalues of
${\hat{D}}_{i}=-i{\hat{\xi}}_{i}^{1}{\hat{\xi}}_{i}^{2}{\hat{\xi}}_{i}^{3}{\hat{\xi}}_{i}^{4}{\hat{c}}_{i}{\hat{d}}_{i}$
are $\pm 1$. For any physical state $|\Psi\rangle_{\rm phys}$, it must then be
${\hat{D}}_{i}|\Psi\rangle_{\rm phys}=|\Psi\rangle_{\rm phys}$, and therefore
one can define the projection operator onto the physical Hilbert subspace as
Chua  _et al._ (2011); Yao  _et al._ (2009)
${\hat{P}}=\prod_{i}(1+{\hat{D}}_{i})/2$.

In the Majorana representation, Eq. (Majorana corner states in square and
Kagome quantum spin liquids) becomes

${\hat{\cal H}_{\rm M}}=i\sum_{\langle i,j\rangle\in{\cal
P}_{\alpha}}[J_{\alpha}{\hat{u}}_{ij}^{\alpha}{\hat{c}}_{i}{\hat{c}}_{j}+J_{\alpha}^{\prime}{\hat{u}}_{ij}^{\alpha}{\hat{d}}_{i}{\hat{d}}_{j}]+iJ_{5}\sum_{i}{\hat{c}}_{i}{\hat{d}}_{i},$
(3)

where $\alpha=1,2$,
${\hat{u}}_{ij}^{1}=-i{\hat{\xi}}_{i}^{1}{\hat{\xi}}_{j}^{2}$, and
${\hat{u}}_{ij}^{2}=-i{\hat{\xi}}_{i}^{3}{\hat{\xi}}_{j}^{4}$. It can be shown
Chua  _et al._ (2011); Yao  _et al._ (2009); Kitaev (2006) that all
${\hat{u}}_{ij}^{\alpha}$ commute with the Hamiltonian (3). Consequently, the
full Hilbert space can be divided into sectors, each obtained by replacing
${\hat{u}}_{ij}^{\alpha}$ with the eigenvalues $u_{ij}^{\alpha}=\pm 1$. In
each sector, Eq. (3) describes free Majorana particles (${\hat{c}}$ and
${\hat{d}}$) propagating on top of a quenched $\mathbb{Z}_{2}$ gauge potential
$u_{ij}^{\alpha}$. The eigenvalues of flux operators are $W_{p}=\prod_{\langle
i,j\rangle\in p}iu_{ij}^{\alpha}\equiv e^{i\phi_{p}}$. Here, $\phi_{p}$ is the
flux through a given plaquette. For Kagome lattices, $\phi_{p}=\pm\pi/2$, if
$p$ is of type 1 or 2, while $\phi_{p}=0$ or $\pi$ if $p$ is of type 3.
Because of the antisymmetry of the Hamiltonian, a total of four distinct flux
patterns exist for any unit cell [see Fig. 1(c)-(f)]. The total flux through a
unit cell is either $0$ or $\pi$, and in the $\pi$ case doubling the unit cell
is required. In the square lattice, the flux through any plaquette can only be
either $\phi_{p}=0$ or $\pi$.

Not all sectors describe different physical systems. In fact, the variables
commuting with the original spin Hamiltonian are the fluxes ${\hat{W}}_{p}$
and not the gauge potential ${\hat{u}}_{ij}^{\alpha}$. By performing the gauge
transformation ${\hat{c}}_{i}\to\Lambda_{i}{\hat{c}}_{i}$,
${\hat{d}}_{i}\to\Lambda_{i}{\hat{d}}_{i}$ and
$u_{ij}\to\Lambda_{i}u_{ij}^{\alpha}\Lambda_{j}$, with $\Lambda_{i}=\pm 1$,
both Hamiltonian and fluxes remain invariant. Hence, $2^{N}$ configurations of
$u_{ij}^{\alpha}$ describe the same physical state Chua  _et al._ (2011); Yao
_et al._ (2009); Kitaev (2006). One can choose to work with any configuration,
depending on convenience: upon projection with ${\hat{P}}$, the physical state
becomes a superposition of states of all equivalent Hilbert space sectors. We
now set $J_{5}=0$ and analyze the two lattices separately.

Square lattice—We rewrite the Hamiltonian (3) as

$\displaystyle{\hat{\cal H}_{\rm M}}$ $\displaystyle=$ $\displaystyle
i\sum_{\ell,m}(-)^{\ell+m}\big{[}u_{\ell,m}^{x}({\tilde{J}}_{m}{\hat{c}}_{\ell,m}{\hat{c}}_{\ell,m+1}+{\tilde{J}}_{m}^{\prime}{\hat{d}}_{\ell,m}{\hat{d}}_{\ell,m+1})$
(4) $\displaystyle-$ $\displaystyle
u_{\ell,m}^{y}({\tilde{J}}_{\ell}{\hat{c}}_{\ell,m}{\hat{c}}_{\ell+1,m}+{\tilde{J}}_{\ell}^{\prime}{\hat{d}}_{\ell,m}{\hat{d}}_{\ell+1,m})\big{]},$

where $\ell$ and $m$ denote the row and column in the lattice of Fig. 1(a),
respectively, and $2{\tilde{J}}_{\ell}=(J_{1}+J_{2})+(-)^{\ell}(J_{1}-J_{2})$
(${\tilde{J}}_{\ell}^{\prime}$ is analogously defined, with $J_{1}^{\prime}$
and $J_{2}^{\prime}$ in lieu of $J_{1}$ and $J_{2}$). In Eq. (4),
$u_{\ell,m}^{x}$ [$u_{\ell,m}^{y}$] is the value of the $\mathbb{Z}_{2}$ gauge
field between sites $i=(\ell,m)$ and $j=(\ell+1,m)$ [$i=(\ell,m)$ and
$j=(\ell,m+1)$] along the $x$ [$y$] direction. Thanks to Lieb’s theorem Lieb
(1994), the ground state of the square lattice with equal hopping amplitudes
is known to contain one flux quantum per plaquette ($\phi_{p}=\pi$)
222Strictly speaking, Lieb’s theorem holds for periodic systems. Following
Kitaev Kitaev (2006), we will assume that the ground state of a large open
lattice coincides with that of the associated periodic system.. When the
difference between hopping amplitudes is much smaller than the two-flux
excitation energy, we expect the ground state configuration to remain
unchanged.

We observe that $-{\hat{\cal H}_{\rm M}}$ has the same flux pattern as
${\hat{\cal H}_{\rm M}}$, since the latter is defined modulo $2\pi$.
Therefore, $-{\hat{\cal H}_{\rm M}}={\hat{\cal G}}{\hat{\cal H}_{\rm
M}}{\hat{\cal G}}^{-1}$, where ${\hat{\cal G}}$ is a gauge transformation that
inverts the signs of all $u_{\ell,m}^{x(y)}$. This implies that the
eigenvalues of ${\hat{\cal H}_{\rm M}}$ must be symmetric about zero, since if
the state $|\psi\rangle$ has energy $E$, then the state ${\hat{\cal
G}}|\psi\rangle$ has energy $-E$. We also note that, if ${\hat{\cal R}}$ is a
$90^{\circ}$ rotation, there exist a gauge operation ${\hat{\cal G}}$ such
that ${\hat{\cal R}}{\hat{\cal H}_{\rm M}}{\hat{\cal R}}^{-1}={\hat{\cal
G}}^{-1}{\hat{\cal H}_{\rm M}}{\hat{\cal G}}$. Therefore ${\hat{\cal
G}}{\hat{\cal R}}$ is a symmetry of the system. It can be easily verified that
$({\hat{\cal G}}{\hat{\cal R}})^{4}=1$, and therefore the eigenstates of (4)
can be one, two or four-fold degenerate.

To progress further, we introduce the fermion operator ${\hat{f}}_{\ell,m}$,
such that
${\hat{c}}_{\ell,m}=i^{\ell+m}{\hat{f}}_{\ell,m}^{\dagger}+(-i)^{\ell+m}{\hat{f}}_{\ell,m}$
and
${\hat{d}}_{\ell,m}=i^{\ell+m+1}{\hat{f}}_{\ell,m}^{\dagger}+(-i)^{\ell+m+1}{\hat{f}}_{\ell,m}$.
Plugging this expressions into Eq. (4) we find

$\displaystyle{\hat{\cal H}_{\rm
M}}=\sum_{\ell,m}\big{\\{}u_{\ell,m}^{x}\big{[}(-)^{\ell+m}t_{m}{\hat{f}}_{\ell,m}^{\dagger}{\hat{f}}_{\ell,m+1}-\Delta_{m}{\hat{f}}_{\ell,m}^{\dagger}{\hat{f}}_{\ell,m+1}^{\dagger}\big{]}$
$\displaystyle-
u_{\ell,m}^{y}\big{[}(-)^{\ell+m}t_{\ell}{\hat{f}}_{\ell,m}^{\dagger}{\hat{f}}_{\ell+1,m}-\Delta_{\ell}{\hat{f}}_{\ell,m}^{\dagger}{\hat{f}}_{\ell+1,m}^{\dagger}\big{]}\big{\\}}+{\rm
h.c.},$ (5)

where $t_{\ell}={\tilde{J}}_{\ell}+{\tilde{J}}_{\ell}^{\prime}$ and
$\Delta_{\ell}={\tilde{J}}_{\ell}-{\tilde{J}}_{\ell}^{\prime}$. The
Hamiltonian (Majorana corner states in square and Kagome quantum spin liquids)
is in form identical to that of spinless electrons paired by p-wave
superconductivity. Here, however, the pairing can have a vortex structure
commensurate to the lattice (see more below). From now on, we choose the
gauge-potential configuration $u_{\ell,m}^{x}=(-1)^{\ell+m}$ and
$u_{\ell,m}^{y}=+1$, that corresponds to having a $\pi$-flux in each
plaquette.

\begin{overpic}[width=130.08731pt]{fig2a} \put(0.0,-10.0){(a)} \end{overpic} |  \begin{overpic}[width=54.2025pt]{fig2b} \put(0.0,-10.0){(b)} \end{overpic} \begin{overpic}[width=54.2025pt]{fig2c} \put(0.0,-10.0){(c)} \end{overpic} \begin{overpic}[width=54.2025pt]{fig2d} \end{overpic}
---|---
Figure 2:  Panel (a) The eigenvalues of Hamiltonian (6) for both $J_{2}=J_{1}$
(inner bands) and $J_{2}>J_{1}$ (outer bands). Panels (b) and (c) The two
corner states of Majorana particles obtained by numerically diagonalizing the
$c$-part of Hamiltonian (4) for $J_{1}=0.3$ and $J_{2}=1$. Due to the
structure of Eq. (4), $d$-particles would exhibit the same states for
$J_{2}^{\prime}>J_{1}^{\prime}$.

Up until now, we have left the couplings $J_{1},J_{2},J_{1}^{\prime}$ and
$J_{2}^{\prime}$ unspecified. Now we study two interesting situations. In the
first case we set $J_{1}^{\prime}=J_{1}$ and $J_{2}^{\prime}=J_{2}$. In turn,
$t_{\ell}=(J_{1}+J_{2})+(-)^{\ell}(J_{1}-J_{2})$ ($t_{m}$ is analogously
defined), while $\Delta_{\ell}=\Delta_{m}=0$. The resulting Hamiltonian is in
form identical to that of a fermionic SOTI Benalcazar  _et al._ (2017a), and
therefore exhibits a phase transition at $J_{1}=J_{2}$. For a system with
periodic boundary conditions we can rewrite ${\hat{\cal H}}_{\rm
M}=\sum_{{\bm{k}},\alpha,\beta}{\hat{F}}_{{\bm{k}},\alpha}^{\dagger}H^{(0)}_{{\bm{k}},\alpha\beta}{\hat{F}}_{{\bm{k}},\beta}$,
where
${\hat{F}}_{{\bm{k}}}={}^{t}({\hat{f}}_{{\bm{k}},1},{\hat{f}}_{{\bm{k}},2},{\hat{f}}_{{\bm{k}},3},{\hat{f}}_{{\bm{k}},4})$
[the unit cell, with the four sites $1=1,\ldots,4$ is shown in Fig. 1(a)].
Here, t denotes transposition, while ($a$ is the side of the unit cell)

$\displaystyle H^{(0)}_{\bm{k}}$ $\displaystyle=$
$\displaystyle\big{[}t_{0}+t_{1}\cos(k_{x}a)\big{]}\openone\otimes\sigma^{x}+\big{\\{}t_{1}\sin(k_{y}a)\sigma^{x}$
(6) $\displaystyle+$
$\displaystyle\big{[}t_{0}-t_{1}\cos(k_{y}a)\big{]}\sigma^{y}+t_{1}\sin(k_{x}a)\sigma^{z}\big{\\}}\otimes\sigma^{y}.$

Here, $\openone$ is the $2\times 2$ identity matrix and $\sigma^{a}$
($a=x,y,z$) are Pauli matrices. Fig. 2(a) shows the eigenvalues of Eq. (6) for
two values of $J_{1}$ and $J_{2}$. When $J_{1}=J_{2}$, the two bands are
doubly-degenerate and touch linearly at the points $(\pm\pi/a,0)$. A gap opens
when $J_{2}\neq J_{1}$. When such system is made finite, it is seen to
transition from trivial insulator to SOTI Benalcazar  _et al._ (2017a) at
$J_{1}=J_{2}$. In Fig. 2(b)-(c), we show the lowest eigenstates of Eq.
(Majorana corner states in square and Kagome quantum spin liquids), defined on
a square lattice with four unit cells along each edge, for $J_{2}>J_{1}$.
These are zero-energy corner states, each of which is doubly degenerate and
exhibits a periodic modulation along the edges reminiscent of SSH edge states.
In fact the weights of corner states alternate between finite and zero going
along each of the two neighboring edges. Edge states also appear at finite
energy. Unlike corner states, these are four-fold degenerate.

More intriguing is the configuration in which $J_{1}^{\prime}=J_{2}$ and
$J_{2}^{\prime}=J_{1}$, which imply $t_{\ell}=t_{m}=J_{1}+J_{2}$ and
$\Delta_{\ell}=(-)^{\ell}(J_{1}-J_{2})$. In this case, both the hopping part
and the p-wave pairing in Eq. (Majorana corner states in square and Kagome
quantum spin liquids) are characterized by a $\pi$-flux structure. The spin
system is therefore mapped into spinless fermions under the simultaneous
effect of $\mathbb{Z}_{2}$ fluxes and p-wave pairing exhibiting a $\pi$-flux
structure (i.e. an Abrikosov lattice of “half vortices” commensurate to the
square lattice). For a system with periodic boundary conditions, the
Hamiltonian is

$\displaystyle{\hat{\cal H}}_{\rm
M}=\sum_{\bm{k}}\left({\hat{F}}_{{\bm{k}}}^{\dagger},{\hat{F}}_{-{\bm{k}}}\right)\left(\\!\begin{array}[]{cc}H^{(0)}_{\bm{k}}&\Delta_{\bm{k}}\vspace{0.1cm}\\\
\Delta^{\dagger}_{\bm{k}}&-{}^{t}H^{(0)}_{-{\bm{k}}}\end{array}\\!\right)\left(\\!\\!\begin{array}[]{c}{\hat{F}}_{{\bm{k}}}\vspace{0.1cm}\\\
{\hat{F}}_{-{\bm{k}}}^{\dagger}\end{array}\\!\\!\right),$ (11)

where

$\displaystyle\Delta_{\bm{k}}$ $\displaystyle=$
$\displaystyle-i\big{[}\Delta_{0}+\Delta_{1}\cos(k_{x}a)\big{]}\openone\otimes\sigma^{y}+i\big{\\{}\Delta_{1}\sin(k_{y}a)\sigma^{x}$
(12) $\displaystyle+$
$\displaystyle\big{[}\Delta_{0}-\Delta_{1}\cos(k_{y}a)\big{]}\sigma^{y}+\Delta_{1}\sin(k_{x}a)\sigma^{z}\big{\\}}\otimes\sigma^{x}.$

Since all hopping amplitudes are equal, if the pairing is artificially made to
vanish the band structure of Eq. (11) exhibits no gap, exactly as in Fig.
2(a). The pairing opens a gap and its vortex structure enables the
localization of Majorana particles at the corners of finite systems
[analogously to Figs. 2(b)-(c)]. Majorana particles are present for all
$J_{1}$ and $J_{2}$, except when $J_{1}=J_{2}$, for which the gap that
stabilizes corner states closes.

All these findings are made clearer by observing that Eq. (4) describes two
copies of the same “Majorana SOTI” Benalcazar  _et al._ (2017a), one for the
${\hat{c}}$\- and one for ${\hat{d}}-$particles. When $J_{5}=0$, the two are
independent and exhibit the SOTI transition at $J_{1}=J_{2}$ and
$J_{1}^{\prime}=J_{2}^{\prime}$, respectively. When $J_{1}^{\prime}=J_{1}$ and
$J_{2}^{\prime}=J_{2}$, both ${\hat{c}}$\- and ${\hat{d}}$-Majorana particles
localize at the corners for $J_{2}>J_{1}$ (when the number of plaquettes per
side is odd). Therefore, as shown above, corner states have fermionic
statistics.

On the contrary, when $J_{1}^{\prime}=J_{2}$ and $J_{2}^{\prime}=J_{1}$,
${\hat{c}}$\- and ${\hat{d}}$-Majorana particles localize at the corners on
opposite sides of the transition point $J_{1}=J_{2}$. Therefore, for
$J_{2}>J_{1}$ ($J_{2}<J_{1}$) only ${\hat{c}}$\- (${\hat{d}}$-)Majorana
particles localize at the corners. Hence, the system is a Majorana SOTI for
all values of $J_{1}\neq J_{2}$. Note that this renders $J_{5}$ largely
ineffective 333See also the Supplemental Online Material. Such coupling pairs
Majorana particles of different kinds at the same site. However, for any value
of $J_{1}\neq J_{2}$ only one type of particle appears at the corners.

Kagome lattice—The ground state of the Kagome lattice is not analytically
known. For this reason, we explore all four configurations shown in Fig.
1(c)-(f). It was indeed shown in Ref. Chua  _et al._ (2011) that these have
very similar energies, and which of them is the ground state can depend on the
size of the lattice. For the sake of presentation, we consider the triangular
lattice shown in Fig. 3. We diagonalize only the ${\hat{c}}$-part of Eq. (3).
In fact, when $J_{5}=0$ the ${\hat{c}}$\- and ${\hat{d}}$-particles are
completely independent and the behavior of the latter can be deduced by that
of the former by replacing $J_{1}\to J_{1}^{\prime}$ and $J_{2}\to
J_{2}^{\prime}$.

\begin{overpic}[width=67.21056pt]{fig3a} \put(0.0,-10.0){(a)} \end{overpic} | \begin{overpic}[width=67.21056pt]{fig3b} \put(0.0,-10.0){(b)} \end{overpic} | \begin{overpic}[width=67.21056pt]{fig3c} \put(0.0,-10.0){(c)} \end{overpic}
---|---|---
Figure 3:  Edge and corner states for flux pattern Fig 1(e). The color scheme
is the same as in Fig. 2(b)-(c). Panel (a) A typical type-one finite energy
edge state, obtained with $J_{1}=0.1$ and $J_{2}=1$. Panel (b) one of the
zero-energy corner states, obtained with the same parameters as Panel (a).
Panel (c) Type-two edge state with $J_{1}=1$ and $J_{2}=0.1$.

In all four cases we find edge states as in Fig. 3(a) coexisting with the
zero-energy corner states of Fig. 3(b) whenever $J_{1}<J_{2}$. We name such
states “type one”, to distinguish them from those of Fig. 3(c) (see below). A
typical type-one edge state is localized along a series of lines parallel to
the edge, and, similar to those of the SSH model, its weights alternate
between zero and nonzero values as one moves farther into the bulk. Similarly,
the weight of type-one corner states alternates between zero and nonzero
values along the two neighboring edges.

Contrary to all other configurations, that of Fig. 1(e) exhibits a second type
of edge states, which we name “type two”. These edge states are considerably
more complicated than the type-one discussed above. Firstly, type-two edge
states occur when one of $J_{1},J_{2}$ is significantly larger than the other
(at $J_{1}=J_{2}$ the system is gapless). Secondly, they are localized along a
line of triangles characterized by strong (either $J_{1}$ or $J_{2}$) bonds.

Conclusion—In this manuscript we have studied two spin-$3/2$ models, defined
on Kagome Chua  _et al._ (2011) and square Yao  _et al._ (2009) lattices, that
support QSL phases at low temperature. As in Kitaev’s work Kitaev (2006),
spins are fractionalized by means of Majorana particles, in this case six of
them. Four of them give rise to a $\mathbb{Z}_{2}$ gauge potential on top of
which the remaining two propagate Chua  _et al._ (2011); Yao  _et al._ (2009).
In the case of the square lattice, the ground-state $\mathbb{Z}_{2}$ flux
configuration is exactly known thanks to a theorem by Lieb Lieb (1994), and
therefore this problem can be solved analytically. On the contrary, since the
Kagome lattice is not bipartite, its ground-state flux configuration is
unknown. To obviate this problem, we have studied the four possible gauge-
inequivalent flux configurations commensurate to the unit cell.

We find that both models can support topologically-protected corner states
when magnetic couplings are made unequal. Such states are fermions when both
free particles resulting from the spin fractionalization are able to localize
at the corners. On the contrary, when only one of the two species can
localize, the corner states are Majorana particles. Notably, in the latter
case such quasiparticles are protected against perturbation that locally mix
the two kinds of free Majorana particles. We observe that the ground-state
energy is lowered by making the magnetic couplings unequal. Therefore, we
speculate that real systems could undergo lattice distortions in order to
lower their magnetic energy and simultaneously localize fermionic or Majorana
states at their corners.

In the Kagome lattice, corner states coexist with edge states, a distinctive
trait of the chiral QSL realized in such system Chua  _et al._ (2011). We find
that all four flux configurations exhibit the same type of edge and corner
states, but one of them supports also a second set of edge states which have
no analog in the other three.

We conclude by noting that corner Majorana states have been predicted to occur
in a similar spin-$3/2$ model defined on a Shastry-Sutherland lattice Dwivedi
_et al._ (2018). The latter is akin to our Kagome model: it exhibits both
chiral and gapped QSL phases but the ground state is not exactly known since
Lieb’s theorem does not apply. On the contrary, the ground state flux
configuration of the square-lattice model we study is exactly known. We
believe that this fact lends credibility to the possibility of localizing
fermionic or Majorana particles at corners of QSLs. Furthermore, we hope that
the relative simplicity of the lattices studied here and their relative
abundance in nature will stimulate the search of material realizations of
second-order QSL supporting topologically-protected corner states.

Acknowledgments—A.P. acknowledges support from the European Commission under
the EU Horizon 2020 MSCA-RISE-2019 programme (project 873028 HYDROTRONICS).

##  Appendix A Protection of edge states against $J_{5}$

\begin{overpic}[width=106.23808pt]{figS1a} \put(0.0,-10.0){(a)} \end{overpic} | \begin{overpic}[width=106.23808pt]{figS1b} \put(0.0,-10.0){(b)} \end{overpic}
---|---
Figure 4:  Energy spectrum of the square model as a function of $J_{5}$. In
both figures, $J_{1}=0.3$ and $J_{2}=1$. Panel (a) Spectrum with
$J_{1}^{\prime}=J_{1}$ and $J_{2}^{\prime}=J_{2}$. Panel (b) Spectrum with
$J_{1}^{\prime}=J_{2}$ and $J_{2}^{\prime}=J_{1}$. It is evident that Majorana
corner states [Panel (b)] are protected against the inter-species coupling
$J_{5}$ up to $|J_{5}|\sim 1.3$.

In this Supplemental Material we show that Majorana corner states are robust
against $J_{5}$ by computing the energy spectrum of the square model for both
particle species as a function of $J_{5}$.

In the case where $J_{1}^{\prime}=J_{1}=0.3$ and $J_{2}^{\prime}=J_{2}=1$
(fermionic corner states), we find that the energies of corner states split
away from zero linearly as $J_{5}$ increases in strength, as shown in Fig.
4(a).

On the contrary, when $J_{1}^{\prime}=J_{2}=1$ and $J_{2}^{\prime}=J_{1}=0.3$
(Majorana corner states), the corner states remain at zero energy for
$|J_{5}|\lesssim 1.3$. We conclude that on-site mixing the of the two species
is forbidden for small-to-moderate inter-species couplings as expected.

## References

  * Anderson (1973) P. Anderson, Materials Research Bulletin 8, 153 (1973).
  * Zhou  _et al._ (2017) Y. Zhou, K. Kanoda, and T.-K. Ng, Rev. Mod. Phys. 89, 025003 (2017).
  * Savary and Balents (2016) L. Savary and L. Balents, Reports on Progress in Physics 80, 016502 (2016).
  * Wen (2004) X. G. Wen,  _Quantum Field Theory of Many-Body Systems: From the Origin of Sound to an Origin of Light and Electrons_ (Oxford University Press, New York, 2004).
  * Kitaev (2006) A. Kitaev, Annals of Physics 321, 2 (2006).
  * Banerjee  _et al._ (2017) A. Banerjee, J. Yan, J. Knolle, C. A. Bridges, M. B. Stone, M. D. Lumsden, D. G. Mandrus, D. A. Tennant, R. Moessner, and S. E. Nagler, Science 356, 1055 (2017).
  * Kitaev (2003) A. Kitaev, Annals of Physics 303, 2 (2003).
  * Sarma  _et al._ (2015) S. D. Sarma, M. Freedman, and C. Nayak, npj Quantum Information 1, 15001 (2015).
  * Hoffman  _et al._ (2016) S. Hoffman, C. Schrade, J. Klinovaja, and D. Loss, Phys. Rev. B 94, 045316 (2016).
  * Lian  _et al._ (2018) B. Lian, X.-Q. Sun, A. Vaezi, X.-L. Qi, and S.-C. Zhang, Proceedings of the National Academy of Sciences 115, 10938 (2018).
  * Knolle (2016) J. Knolle,  _Dynamics of a Quantum Spin Liquid_ (Springer, Heidelberg, 2016).
  * Benalcazar  _et al._ (2017a) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Science 357, 61 (2017a).
  * Benalcazar  _et al._ (2017b) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Phys. Rev. B 96, 245115 (2017b).
  * Schindler  _et al._ (2018a) F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. P. Parkin, B. A. Bernevig, and T. Neupert, Science Advances 4 (2018a), 10.1126/sciadv.aat0346.
  * Langbehn  _et al._ (2017) J. Langbehn, Y. Peng, L. Trifunovic, F. von Oppen, and P. W. Brouwer, Phys. Rev. Lett. 119, 246401 (2017).
  * Geier  _et al._ (2018) M. Geier, L. Trifunovic, M. Hoskam, and P. W. Brouwer, Phys. Rev. B 97, 205135 (2018).
  * Ezawa (2018a) M. Ezawa, Phys. Rev. Lett. 120, 026801 (2018a).
  * Song  _et al._ (2017) Z. Song, Z. Fang, and C. Fang, Phys. Rev. Lett. 119, 246402 (2017).
  * Ezawa (2018b) M. Ezawa, Phys. Rev. B 97, 155305 (2018b).
  * Ezawa (2018c) M. Ezawa, Phys. Rev. B 97, 241402 (2018c).
  * Serra-Garcia  _et al._ (2018) M. Serra-Garcia, V. Peri, R. Süsstrunk, O. R. Bilal, T. Larsen, L. G. Villanueva, and S. D. Huber, Nature 555, 342 (2018).
  * Imhof  _et al._ (2018) S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, and R. Thomale, Nature Physics 14, 925 (2018).
  * Peterson  _et al._ (2018) C. W. Peterson, W. A. Benalcazar, T. L. Hughes, and G. Bahl, Nature 555, 346 (2018).
  * El Hassan  _et al._ (2019) A. El Hassan, F. K. Kunst, A. Moritz, G. Andler, E. J. Bergholtz, and M. Bourennane, Nature Photonics 13, 697 (2019).
  * Schindler  _et al._ (2018b) F. Schindler, Z. Wang, M. G. Vergniory, A. M. Cook, A. Murani, S. Sengupta, A. Y. Kasumov, R. Deblock, S. Jeon, I. Drozdov, H. Bouchiat, S. Guéron, A. Yazdani, B. A. Bernevig, and T. Neupert, Nature Physics 14, 918 (2018b).
  * Luo and Zhang (2019) X.-W. Luo and C. Zhang, Phys. Rev. Lett. 123, 073601 (2019).
  * Liu  _et al._ (2018) T. Liu, J. J. He, and F. Nori, Phys. Rev. B 98, 245413 (2018).
  * Zhu (2018) X. Zhu, Phys. Rev. B 97, 205134 (2018).
  * Laubscher  _et al._ (2019) K. Laubscher, D. Loss, and J. Klinovaja, Phys. Rev. Research 1, 032017 (2019).
  * Wang  _et al._ (2018) Q. Wang, C.-C. Liu, Y.-M. Lu, and F. Zhang, Phys. Rev. Lett. 121, 186801 (2018).
  * Kheirkhah  _et al._ (2020) M. Kheirkhah, Y. Nagai, C. Chen, and F. Marsiglio, Phys. Rev. B 101, 104502 (2020).
  * Yan (2019) Z. Yan, Phys. Rev. Lett. 123, 177001 (2019).
  * Dwivedi  _et al._ (2018) V. Dwivedi, C. Hickey, T. Eschmann, and S. Trebst, Phys. Rev. B 98, 054432 (2018).
  * Note (1) The model of Ref. Dwivedi  _et al._ (2018) can be solved exactly by mapping spins into Majorana particles which are free to propagate on top a quenched $\mathbb{Z}_{2}$ gauge potential. However, the ground-state flux configuration is not analytically known and must be found with numerical techniques.
  * Chua  _et al._ (2011) V. Chua, H. Yao, and G. A. Fiete, Phys. Rev. B 83, 180412 (2011).
  * Yao  _et al._ (2009) H. Yao, S.-C. Zhang, and S. A. Kivelson, Phys. Rev. Lett. 102, 217202 (2009).
  * Lieb (1994) E. H. Lieb, Phys. Rev. Lett. 73, 2158 (1994).
  * Note (2) Strictly speaking, Lieb’s theorem holds for periodic systems. Following Kitaev Kitaev (2006), we will assume that the ground state of a large open lattice coincides with that of the associated periodic system.
  * Note (3) See also the Supplemental Online Material.


$\gamma^{(2)}(\bm{\phi}_{2})=L_{\widetilde{S1},\overline{12},\overline{2U}}T_{S}\mathbf{v}_{2}^{H}\mathbf{G}_{22}^{H}\mathbf{G}_{22}\mathbf{v}_{2}+2\sqrt{L_{\widetilde{S1},\overline{1U}}L_{\widetilde{S1},\overline{12},\overline{2U}}}T_{S}\mathfrak{R}\left\\{\mathbf{v}_{2}^{H}\mathbf{G}_{22}^{H}\mathbf{g}_{22}\right\\}+\gamma^{(0)},$

with
$\mathbf{G}_{22}\triangleq\bar{\mathbf{H}}_{12}^{H}\mathrm{diag}(\mathbf{a}_{D,12})$
and $\mathbf{g}_{12}\triangleq\bar{\mathbf{h}}_{1U}$. By (3), (4), and Lemma
13, we have
$\bm{\phi}_{1}^{\star}=\Lambda\big{(}-\angle(\mathrm{diag}(\mathbf{a}^{H}_{D,12})\mathbf{a}_{A,12})$
$-\angle(\mathbf{a}_{D,12}^{H}\mathbf{a}_{D,1U})\mathbf{1}_{T_{2}}\big{)}\overset{(h)}{=}\Lambda(-\bm{\Delta}_{\overline{12},\overline{2U}}$
$-\angle(r_{\overline{12},\overline{1U}})\mathbf{1}_{T_{2}})$, where $(h)$ is
due to (10) and (11).

Therefore, we complete the proof of Theorem 2.

## Appendix D: Proof of Lemma 2

By (14), (15), (22), and (23), we have:

$\displaystyle\gamma^{(l)}(\phi_{l,t},\bm{\phi}_{l,-t})\overset{}{=}\left(\sum_{k\in\mathcal{T}_{l},k\neq
t\atop i\in\mathcal{T}_{l},i\neq
t}A_{l,k,i}e^{j(\phi_{l,k}-\phi_{l,i})}+A_{l,t,t}\right)+\left(e^{j\phi_{l,t}}\sum_{k\in\mathcal{T}_{l},k\neq
t}A_{l,t,k}e^{-j\phi_{l,k}}+e^{-j\phi_{l,t}}\sum_{k\in\mathcal{T}_{l},k\neq
t}A_{l,k,t}e^{j\phi_{l,k}}\right)$
$\displaystyle\qquad\qquad\qquad\quad+2\mathfrak{R}\left\\{\left(e^{j\phi_{l,t}}b_{l,t}+\sum_{k\in\mathcal{T}_{l},k\neq
t}e^{j\phi_{l,k}}b_{l,k}\right)\right\\}+\gamma^{(0)}$
$\displaystyle\overset{(a)}{=}2\mathfrak{R}\left\\{e^{j\phi_{l,t}}\left(\sum_{k\in\mathcal{T}_{l},k\neq
t}A_{l,t,k}e^{-j\phi_{l,k}}+b_{l,t}\right)+\sum_{k\in\mathcal{T}_{l},k\neq
t}e^{j\phi_{l,k}}b_{l,k}\right\\}+\sum_{k\in\mathcal{T}_{l},k\neq t\atop
i\in\mathcal{T}_{l},i\neq
t}A_{l,k,i}e^{j(\phi_{l,k}-\phi_{l,i})}+A_{l,t,t}+\gamma^{(0)},$

where $(a)$ is due to $\mathbf{A}_{l}=\mathbf{A}_{l}^{H}$. As
$\bm{\phi}_{l,-t}$ is given, to solve Problem 2 is equivalent to solve:

$\displaystyle\max_{\phi_{l,t}\in[0,2\pi)}\mathfrak{R}\left\\{e^{j\phi_{l,t}}\left(\sum_{k\in\mathcal{T}_{l},k\neq
t}A_{l,t,k}e^{-j\phi_{l,k}}+b_{l,t}\right)\right\\}.$

Besides, we have:

$\displaystyle\mathfrak{R}\left\\{e^{j\phi_{l,t}}\left(\sum_{k\in\mathcal{T}_{l},k\neq
t}A_{l,t,k}e^{-j\phi_{l,k}}+b_{l,t}\right)\right\\}\overset{(b)}{\leq}\left|e^{j\phi_{l,t}}\left(\sum_{k\in\mathcal{T}_{l},k\neq
t}A_{l,t,k}e^{-j\phi_{l,k}}+b_{l,t}\right)\right|\overset{}{=}\left|\sum_{k\in\mathcal{T}_{l},k\neq
t}A_{l,t,k}e^{-j\phi_{l,k}}+b_{l,t}\right|,$

where $(b)$ is due to $\mathfrak{R}\left\\{x\right\\}\leq|x|$,
$x\in\mathbb{C}$ and holds with equality if and only if
$\phi_{l,t}=-\angle\left(\sum_{k\in\mathcal{T}_{l},k\neq
t}A_{l,t,k}e^{-j\phi_{l,k}}+b_{l,t}\right)$. By noting that
$\phi_{l,t}\in[0,2\pi)$, the unique optimal solution of the equivalent problem
is
$\phi_{l,t}^{\dagger}=\Lambda\left(-\angle\left(\sum_{k\in\mathcal{T}_{l},k\neq
t}A_{l,t,k}e^{-j\phi_{l,k}}+b_{l,t}\right)\right)$. Therefore, we complete the
proof of Lemma 2.

## Appendix E: Proof of Theorem 3

First, we prove statement (i). Consider $K_{S2}=K_{1U}=0$. We first have:

$\displaystyle\left|\mathbf{v}_{1}^{H}\mathrm{diag}(\mathbf{a}^{H}_{D,12})\mathbf{a}_{A,S1}\right|$
$\displaystyle\overset{(a)}{=}\left|\sum_{t\in\mathcal{T}_{1}}e^{j\left(\phi_{1,t}+\Delta_{\overline{S1},\overline{12},t}\right)}\right|\overset{(b)}{\leq}T_{1},$
(70)
$\displaystyle\left|\mathbf{v}_{2}^{H}\mathrm{diag}(\mathbf{h}^{H}_{2U})\mathbf{a}_{A,12}\right|$
$\displaystyle\overset{(c)}{=}\left|\sum_{t\in\mathcal{T}_{2}}e^{j\left(\phi_{2,t}+\Delta_{\overline{12},\overline{2U},t}\right)}\right|\overset{(d)}{\leq}T_{2},$
(71)

where $(a)$ and $(c)$ are due to (10), and $(b)$ and $(d)$ are due to the
triangle inequality. Note that $(b)$ and $(d)$ hold with equality if and only
if:

$\displaystyle\bm{\phi}_{1}$
$\displaystyle=-\bm{\Delta}_{\overline{S1},\overline{12}}+\psi_{1}\mathbf{1}_{T_{1}},\psi_{1}\in\mathbb{R},$
(72) $\displaystyle\bm{\phi}_{2}$
$\displaystyle=-\bm{\Delta}_{\overline{12},\overline{2U}}+\psi_{2}\mathbf{1}_{T_{2}},\psi_{2}\in\mathbb{R}.$
(73)

Then, we have:

$\displaystyle\gamma^{(3)}(\bm{\phi}_{1},\bm{\phi}_{2})$
$\displaystyle\overset{(e)}{=}L_{\overline{S1},\overline{12},\widetilde{2U}}T_{S}T_{2}\left|\mathbf{v}_{1}^{H}\mathrm{diag}(\mathbf{a}^{H}_{D,12})\mathbf{a}_{A,S1}\right|^{2}+L_{\widetilde{S1},\overline{12},\overline{2U}}T_{S}T_{1}\left|\mathbf{v}_{2}^{H}\mathrm{diag}(\mathbf{h}_{2U}^{H})\mathbf{a}_{A,12}\right|^{2}$
$\displaystyle\quad+L_{\overline{S1},\overline{12},\overline{2U}}T_{S}\left|\mathbf{v}_{2}^{H}\mathrm{diag}(\mathbf{h}_{2U}^{H})\mathbf{a}_{A,12}\mathbf{a}^{H}_{D,12}\mathrm{diag}(\mathbf{a}_{A,S1})\mathbf{v}_{1}^{*}\right|^{2}$
$\displaystyle\quad+2\sqrt{L_{\overline{SU}}L_{\overline{S1},\overline{12},\overline{2U}}}\mathfrak{R}\left\\{\mathbf{v}_{2}^{H}\mathrm{diag}(\mathbf{h}_{2U}^{H})\mathbf{a}_{A,12}\mathbf{a}^{H}_{D,12}\mathrm{diag}(\mathbf{a}_{A,S1})\mathbf{v}_{1}^{*}r_{\overline{S1},\overline{SU}}\right\\}+\gamma^{(0)}$
$\displaystyle\stackrel{{\scriptstyle(f)}}{{\leq}}L_{\overline{S1},\overline{12},\widetilde{2U}}T_{S}T_{2}\left|\mathbf{v}_{1}^{H}\mathrm{diag}(\mathbf{a}^{H}_{D,12})\mathbf{a}_{A,S1}\right|^{2}+L_{\widetilde{S1},\overline{12},\overline{2U}}T_{S}T_{1}\left|\mathbf{v}_{2}^{H}\mathrm{diag}(\mathbf{h}_{2U}^{H})\mathbf{a}_{A,12}\right|^{2}$
$\displaystyle\quad+L_{\overline{S1},\overline{12},\overline{2U}}T_{S}\left|\mathbf{v}_{2}^{H}\mathrm{diag}(\mathbf{h}_{2U}^{H})\mathbf{a}_{A,12}\right|^{2}\left|\mathbf{v}_{1}^{H}\mathrm{diag}(\mathbf{a}^{H}_{D,12})\mathbf{a}_{A,S1}\right|^{2}$
$\displaystyle\quad+2\sqrt{L_{\overline{SU}}L_{\overline{S1},\overline{12},\overline{2U}}}\left|\mathbf{v}_{2}^{H}\mathrm{diag}(\mathbf{h}_{2U}^{H})\mathbf{a}_{A,12}\right|\left|\mathbf{v}_{1}^{H}\mathrm{diag}(\mathbf{a}^{H}_{D,12})\mathbf{a}_{A,S1}\right|\left|r_{\overline{S1},\overline{SU}}\right|+\gamma^{(0)}$
$\displaystyle\stackrel{{\scriptstyle(g)}}{{\leq}}L_{\overline{S1},\overline{12},\widetilde{2U}}T_{S}T_{2}T_{1}^{2}+L_{\widetilde{S1},\overline{12},\overline{2U}}T_{S}T_{1}T_{2}^{2}+L_{\overline{S1},\overline{12},\overline{2U}}T_{S}T_{1}^{2}T_{2}^{2}$
$\displaystyle\quad+2\sqrt{L_{\overline{SU}}L_{\overline{S1},\overline{12},\overline{2U}}}T_{1}T_{2}\left|r_{\overline{S1},\overline{SU}}\right|+\gamma^{(0)},$
(74)

where $(e)$ is due to (16), $(f)$ is due to
$\mathbf{a}^{H}_{D,12}\mathrm{diag}(\mathbf{a}_{A,S1})\mathbf{v}_{1}^{*}=\mathbf{v}_{1}^{H}\mathrm{diag}(\mathbf{a}^{H}_{D,12})\mathbf{a}_{A,S1}$
and $\mathfrak{R}\\{x\\}\leq|x|,x\in\mathbb{C}$, and $(g)$ is due to (70) and
(71). Note that $(f)$ holds with equality if and only if

$\displaystyle\angle\left(\mathbf{v}_{2}^{H}\mathrm{diag}(\mathbf{h}_{2U}^{H})\mathbf{a}_{A,12}\right)+\angle\left(\mathbf{v}_{1}^{H}\mathrm{diag}(\mathbf{a}^{H}_{D,12})\mathbf{a}_{A,S1}\right)+\angle\left(r_{\overline{S1},\overline{SU}}\right)=0.$
(75)

Besides, note that $(g)$ holds with equality if and only if (70) and (71) hold
with equality, i.e., (72) and (73) hold. Thus, the upper bound in (74) is
achieved if and only if (72), (73), and (75) hold simultaneously. By
substituting (72) and (73) into (75), we have
$\psi_{1}+\psi_{2}+\angle\left(r_{\overline{S1},\overline{SU}}\right)=0$. By
noting that $\bm{\phi}_{l,t}\in[0,2\pi),t\in\mathcal{T}_{l},l\in\mathcal{L}$
and by choosing
$\psi_{1}=-\frac{1}{2}\angle(r_{\overline{S1},\overline{SU}})+\psi$ and
$\psi_{2}=-\frac{1}{2}\angle(r_{\overline{S1},\overline{SU}})-\psi$ for all
$\psi\in\mathbb{R}$, we can show that any
$(\bm{\phi}_{1}^{\star},\bm{\phi}_{2}^{\star})$ with
$\bm{\phi}_{1}^{\star}=\Lambda\left(-\bm{\Delta}_{\overline{S1},\overline{12}}+\psi\mathbf{1}_{T_{1}}-\frac{1}{2}\angle\left(r_{\overline{S1},\overline{SU}}\right)\mathbf{1}_{T_{1}}\right)$
and
$\bm{\phi}_{2}^{\star}=\Lambda\left(-\bm{\Delta}_{\overline{12},\overline{2U}}-\psi\mathbf{1}_{T_{2}}-\frac{1}{2}\angle\left(r_{\overline{S1},\overline{SU}}\right)\mathbf{1}_{T_{2}}\right)$,
$\psi\in\mathbb{R}$ is an optimal solution of Problem 1.

Next, we prove statement (ii). Consider $K_{SU}=K_{12}=0$. We first have:

$\displaystyle\left|\mathbf{v}_{l}^{H}\mathrm{diag}(\mathbf{h}^{H}_{lU})\mathbf{a}_{A,Sl}\right|\overset{(h)}{=}\left|\sum_{t\in\mathcal{T}_{l}}e^{j\left(\phi_{l,t}+\Delta_{\overline{Sl},\overline{lU},t}\right)}\right|\overset{(i)}{\leq}T_{l},l\in\mathcal{L},$
(76)

where $(h)$ is due to (10), $(i)$ is due to the triangle inequality and holds
with equality if and only if:

$\displaystyle\bm{\phi}_{l}=\Lambda\left(-\bm{\Delta}_{\overline{Sl},\overline{lU}}+\psi_{l}\right),\psi_{l}\in\mathbb{R}.$
(77)

Then, we have:

$\displaystyle\gamma^{(3)}(\bm{\phi}_{1},\bm{\phi}_{2})\overset{(j)}{=}$
$\displaystyle\sum_{l=1}^{2}L_{\overline{Sl},\overline{lU}}T_{S}\left|\mathbf{v}_{l}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{lU}^{H})\mathbf{a}_{A,Sl}\right|^{2}$
$\displaystyle+2\sqrt{L_{\overline{S1},\overline{1U}}L_{\overline{S2},\overline{2U}}}\mathfrak{R}\left\\{\mathbf{v}_{2}^{H}\mathrm{diag}\left(\bar{\mathbf{h}}_{2U}^{H}\right)\bar{\mathbf{H}}_{S2}\bar{\mathbf{H}}_{S1}^{H}\mathrm{diag}\left(\bar{\mathbf{h}}_{1U}\right)\mathbf{v}_{1}\right\\}+\gamma^{(0)}$
$\displaystyle\overset{(k)}{\leq}$
$\displaystyle\sum_{l=1}^{2}L_{\overline{Sl},\overline{lU}}T_{S}\left|\mathbf{v}_{l}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{lU}^{H})\mathbf{a}_{A,Sl}\right|^{2}+\gamma^{(0)}$
$\displaystyle+2\sqrt{L_{\overline{S1},\overline{1U}}L_{\overline{S2},\overline{2U}}}\left|\mathbf{v}_{2}^{H}\mathrm{diag}\left(\bar{\mathbf{h}}_{2U}^{H}\right)\mathbf{a}_{A,S2}\right|\left|\mathbf{v}_{1}^{H}\mathrm{diag}\left(\bar{\mathbf{h}}_{1U}^{H}\right)\mathbf{a}_{A,S1}\right|\left|r_{\overline{S1},\overline{S2}}\right|$
$\displaystyle\overset{(l)}{\leq}$
$\displaystyle\sum_{l=1}^{2}L_{\overline{Sl},\overline{lU}}T_{S}T_{l}^{2}+2\sqrt{L_{\overline{S1},\overline{1U}}L_{\overline{S2},\overline{2U}}}T_{1}T_{2}\left|r_{\overline{S1},\overline{S2}}\right|+\gamma^{(0)},$
(78)

where $(j)$ is due to (16), $(k)$ is due to
$\mathfrak{R}\\{x\\}\leq|x|,x\in\mathbb{C}$, and $(l)$ is due to (76). Note
that $(k)$ holds with equality if and only if:

$\displaystyle\angle\left(\mathbf{v}_{1}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{1U}^{H})\mathbf{a}_{A,S1}\right)-\angle\left(\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\mathbf{a}_{A,S2}\right)+\angle\left(r_{\overline{S1},\overline{S2}}\right)=0.$
(79)

Besides, note that $(l)$ holds with equality if and only if (76) holds with
equality, i.e., (77) holds. Thus, the upper bound in (78) is achieved if and
only if (77) and (79) hold simultaneously. By substituting (77) into (79), we
have $\psi_{1}-\psi_{2}+\angle\left(r_{\overline{S1},\overline{S2}}\right)=0$.
By noting that
$\bm{\phi}_{l,t}\in[0,2\pi),t\in\mathcal{T}_{l},l\in\mathcal{L}$ and by
choosing $\psi_{1}=\psi$ and
$\psi_{2}=\psi+\angle(r_{\overline{S1},\overline{S2}})$ for all
$\psi\in\mathbb{R}$, we can show that any
$(\bm{\phi}_{1}^{\star},\bm{\phi}_{2}^{\star})$ with
$\phi^{\star}_{1}=\Lambda(-\bm{\Delta}_{\overline{S1},\overline{1U}}+\psi\mathbf{1}_{T_{1}})$,
$\phi^{\star}_{2}=\Lambda(-\bm{\Delta}_{\overline{S2},\overline{2U}}+\psi\mathbf{1}_{T_{2}}+\angle(r_{\overline{S1},\overline{S2}})\mathbf{1}_{T_{2}})$,
$\psi\in\mathbb{R}$ is an optimal solution of Problem 1.

Therefore, we complete the proof of Theorem 3.

## Appendix F: Proof of Lemma 3

By (16), (25), and (26), we have:

$\displaystyle\gamma^{(3)}(\phi_{l,t},\bm{\phi}_{l,-t},\phi_{-l})\overset{}{=}\left(\sum_{k\in\mathcal{T}_{l},k\neq
t\atop i\in\mathcal{T}_{l},i\neq
t}C_{l,k,i}e^{j(\phi_{l,k}-\phi_{l,i})}+C_{l,t,t}\right)+\left(e^{j\phi_{l,t}}\sum_{k\in\mathcal{T}_{l},k\neq
t}C_{l,t,k}e^{-j\phi_{l,k}}+e^{-j\phi_{l,t}}\sum_{k\in\mathcal{T}_{l},k\neq
t}C_{l,k,t}e^{j\phi_{l,k}}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\quad\quad+2\mathfrak{R}\left\\{\left(e^{j\phi_{l,t}}d_{l,t}+\sum_{k\in\mathcal{T}_{l},k\neq
t}e^{j\phi_{l,k}}d_{l,k}\right)\right\\}+\gamma^{(0)}$
$\displaystyle\overset{(a)}{=}2\mathfrak{R}\left\\{e^{j\phi_{l,t}}\left(\sum_{k\in\mathcal{T}_{l},k\neq
t}C_{l,t,k}e^{-j\phi_{l,k}}+d_{l,t}\right)+\sum_{k\in\mathcal{T}_{l},k\neq
t}e^{j\phi_{l,k}}d_{l,k}\right\\}+\sum_{k\in\mathcal{T}_{l},k\neq t\atop
i\in\mathcal{T}_{l},i\neq
t}C_{l,k,i}e^{j(\phi_{l,k}-\phi_{l,i})}+C_{l,t,t}+\gamma^{(0)},$

where $(a)$ is due to $\mathbf{C}_{l}=\mathbf{C}_{l}^{H}$. As
$\left(\bm{\phi}_{l,-t},\phi_{-l}\right)$ is given, to solve Problem 4 is
equivalent to solve:

$\displaystyle\max_{\phi_{l,t}\in[0,2\pi)}\mathfrak{R}\left\\{e^{j\phi_{l,t}}\left(\sum_{k\in\mathcal{T}_{l},k\neq
t}C_{l,t,k}e^{-j\phi_{l,k}}+d_{l,t}\right)\right\\}.$

Besides, we have:

$\displaystyle\mathfrak{R}\left\\{e^{j\phi_{l,t}}\left(\sum_{k\in\mathcal{T}_{l},k\neq
t}C_{l,t,k}e^{-j\phi_{l,k}}+d_{l,t}\right)\right\\}\overset{(b)}{\leq}\left|\sum_{k\in\mathcal{T}_{l},k\neq
t}C_{l,t,k}e^{-j\phi_{l,k}}+d_{l,t}\right|,$

where $(b)$ is due to $\mathfrak{R}\\{x\\}\leq|x|$, $x\in\mathbb{C}$ and holds
with equality if and only if
$\phi_{l,t}=-\angle\left(\sum_{k\in\mathcal{T}_{l},k\neq
t}C_{l,t,k}e^{-j\phi_{l,k}}+d_{l,t}\right)$. By noting that
$\phi_{l,t}\in[0,2\pi)$, the unique optimal solution of the equivalent problem
is
$\phi_{l,t}^{\dagger}=\Lambda\left(-\angle\left(\sum_{k\in\mathcal{T}_{l},k\neq
t}C_{l,t,k}e^{-j\phi_{l,k}}+d_{l,t}\right)\right)$. Therefore, we complete the
proof of Lemma 3.

## Appendix G: Proof of Theorem 4

Before proving Theorem 4, we first show:

$\displaystyle\|\mathbf{v}_{l}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{lU}^{H})\bar{\mathbf{H}}_{Sl}\|_{2}$
$\displaystyle=\|\mathbf{v}_{l}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{lU}^{H})\mathbf{a}_{A,Sl}\mathbf{a}^{H}_{D,Sl}\|_{2}$
$\displaystyle\overset{(a)}{\leq}|\mathbf{v}_{l}^{H}\mathrm{diag}(\mathbf{a}_{D,lU}^{H})\mathbf{a}_{A,Sl}|\|\mathbf{a}^{H}_{D,Sl}\|_{2}\overset{(b)}{\leq}\sqrt{T_{S}}T_{l},$
(80)
$\displaystyle|\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{12}\mathrm{diag}(\mathbf{a}_{A,S1})\mathbf{v}_{1}^{*}|$
$\displaystyle=\|\mathbf{v}_{2}^{H}\mathrm{diag}(\mathbf{a}_{D,2U}^{H})\mathbf{a}_{A,12}\mathbf{a}^{H}_{D,12}\mathrm{diag}(\mathbf{a}_{A,S1})\mathbf{v}_{1}^{*}\|_{2}$
$\displaystyle\overset{(c)}{\leq}|\mathbf{v}_{2}^{H}\mathrm{diag}(\mathbf{a}_{D,2U}^{H})\mathbf{a}_{A,12}||\mathbf{a}^{H}_{D,12}\mathrm{diag}(\mathbf{a}_{A,S1})\mathbf{v}_{1}^{*}|\overset{(d)}{\leq}T_{1}T_{2},$
(81)
$\displaystyle\|\mathbf{v}_{1}^{H}\mathrm{diag}(\mathbf{a}^{H}_{D,12})\bar{\mathbf{H}}_{S1}\|_{2}$
$\displaystyle=\|\mathbf{v}_{1}^{H}\mathrm{diag}(\mathbf{a}^{H}_{D,12})\mathbf{a}_{A,S1}\mathbf{a}^{H}_{D,S1}\|_{2}$
$\displaystyle\overset{(e)}{\leq}|\mathbf{v}_{1}^{H}\mathrm{diag}(\mathbf{a}^{H}_{D,12})\mathbf{a}_{A,S1}|\|\mathbf{a}^{H}_{D,S1}\|_{2}\overset{(f)}{\leq}\sqrt{T_{S}}T_{1},$
(82)
$\displaystyle\|\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U})\bar{\mathbf{H}}_{12}\|_{2}$
$\displaystyle=\|\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U})\mathbf{a}_{A,12}\mathbf{a}^{H}_{D,12}\|_{2}$
$\displaystyle\overset{(g)}{\leq}|\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U})\mathbf{a}_{A,12}|\|\mathbf{a}^{H}_{D,12}\|_{2}\overset{(h)}{\leq}\sqrt{T_{1}}T_{2},$
(83)

where $(a)$, $(c)$, $(e)$, and $(g)$ are due to Cauchy-Schwartz inequality,
and $(b)$, $(d)$, $(f)$, and $(h)$ are due to the triangle inequality. Based
on (80), (81), (82), and (83), we prove Theorem 4 below.

First, consider Case 0. As
$\lim_{T_{1},T_{2}\rightarrow\infty}\frac{\gamma^{(0)}}{(L_{\widetilde{S1},\overline{12},\widetilde{2U}}+L_{\overline{S1},\widetilde{12},\overline{2U}}+L_{\overline{S1},\widetilde{12},\widetilde{2U}}+L_{\widetilde{S1},\widetilde{12},\overline{2U}}+L_{\widetilde{S1},\widetilde{12},\widetilde{2U}})T_{S}T_{1}T_{2}}=1,$
we can show $\gamma^{(0)}\stackrel{{\scriptstyle
T_{1},T_{2}\rightarrow\infty}}{{\sim}}(L_{\widetilde{S1},\overline{12},\widetilde{2U}}+L_{\overline{S1},\widetilde{12},\overline{2U}}+L_{\overline{S1},\widetilde{12},\widetilde{2U}}+L_{\widetilde{S1},\widetilde{12},\overline{2U}}+L_{\widetilde{S1},\widetilde{12},\widetilde{2U}})T_{S}T_{1}T_{2}$.

Next, consider Case 1. For all $\bm{\phi}_{1}$, we have:

$\displaystyle\gamma^{(1)}(\bm{\phi}_{1})\overset{(i)}{=}L_{\overline{S1},\overline{1U}}\|\mathbf{v}_{1}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{1U}^{H})\bar{\mathbf{H}}_{S1}\|_{2}^{2}+L_{\overline{S1},\overline{12},\widetilde{2U}}T_{2}\left\|\mathbf{v}_{1}^{H}\mathrm{diag}(\mathbf{a}^{H}_{D,12})\bar{\mathbf{H}}_{S1}\right\|_{2}^{2}+\gamma^{(0)}$
$\displaystyle\qquad\qquad\quad+2\mathfrak{R}\left\\{\mathbf{v}_{1}^{H}\left(\sqrt{L_{\overline{SU}}L_{\overline{S1},\overline{1U}}}\mathrm{diag}(\bar{\mathbf{h}}_{1U}^{H})\bar{\mathbf{H}}_{S1}\bar{\mathbf{h}}_{SU}+\sqrt{L_{\overline{S2},\widetilde{2U}}L_{\overline{S1},\overline{12},\widetilde{2U}}}\mathrm{diag}\left(\mathbf{a}_{D,12}^{H}\right)\bar{\mathbf{H}}_{S1}\bar{\mathbf{H}}_{S2}^{H}\mathbf{a}_{A,12}\right)\right\\}$
$\displaystyle\stackrel{{\scriptstyle(j)}}{{\leq}}L_{\overline{S1},\overline{1U}}\|\mathbf{v}_{1}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{1U}^{H})\bar{\mathbf{H}}_{S1}\|_{2}^{2}+L_{\overline{S1},\overline{12},\widetilde{2U}}T_{2}\left\|\mathbf{v}_{1}^{H}\mathrm{diag}(\mathbf{a}^{H}_{D,12})\bar{\mathbf{H}}_{S1}\right\|_{2}^{2}+\gamma^{(0)}$
$\displaystyle\quad+2\sqrt{L_{\overline{SU}}L_{\overline{S1},\overline{1U}}}\|\mathbf{v}_{1}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{1U}^{H})\bar{\mathbf{H}}_{S1}\|_{2}\|\mathbf{a}_{D,SU}\|_{2}+2\sqrt{L_{\overline{S2},\widetilde{2U}}L_{\overline{S1},\overline{12},\widetilde{2U}}}\|\mathbf{v}_{1}^{H}\mathrm{diag}\left(\mathbf{a}_{D,12}^{H}\right)\bar{\mathbf{H}}_{S1}\|_{2}\|\mathbf{a}_{D,S2}\|_{2}|r_{S2,12}|$
$\displaystyle\stackrel{{\scriptstyle(k)}}{{\leq}}L_{\overline{S1},\overline{1U}}T_{S}T_{1}^{2}+L_{\overline{S1},\overline{12},\widetilde{2U}}T_{S}T_{1}^{2}T_{2}+2\sqrt{L_{\overline{SU}}L_{\overline{S1},\overline{1U}}}T_{S}T_{1}+2\sqrt{L_{\overline{S2},\widetilde{2U}}L_{\overline{S1},\overline{12},\widetilde{2U}}}T_{S}T_{1}T_{2}+\gamma^{(0)}\triangleq\gamma^{(1)\star}_{ub},$
(84)

where $(i)$ is due to (14), $(j)$ is due to
$\mathfrak{R}\\{x\\}\leq|x|,x\in\mathbb{C}$ and Cauchy-Swartz inequality, and
$(k)$ is due to (80) and (82). On the other hand, we have:

$\displaystyle\gamma^{(1)}\left(\bm{\hat{}}{\phi}_{1}\right)=L_{\overline{S1},\overline{1U}}\|\hat{\mathbf{v}}_{1}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{1U}^{H})\bar{\mathbf{H}}_{S1}\|_{2}^{2}+L_{\overline{S1},\overline{12},\widetilde{2U}}T_{2}\left\|\hat{\mathbf{v}}_{1}^{H}\mathrm{diag}(\mathbf{a}^{H}_{D,12})\bar{\mathbf{H}}_{S1}\right\|_{2}^{2}+\gamma^{(0)}$
$\displaystyle\qquad\qquad\quad+2\mathfrak{R}\left\\{\hat{\mathbf{v}}_{1}^{H}\left(\sqrt{L_{\overline{SU}}L_{\overline{S1},\overline{1U}}}\mathrm{diag}(\bar{\mathbf{h}}_{1U}^{H})\bar{\mathbf{H}}_{S1}\bar{\mathbf{h}}_{SU}+\sqrt{L_{\overline{S2},\widetilde{2U}}}\mathrm{diag}\left(\mathbf{a}_{D,12}^{H}\right)\bar{\mathbf{H}}_{S1}\bar{\mathbf{H}}_{S2}^{H}\mathbf{a}_{A,12}\right)\right\\}$
$\displaystyle\stackrel{{\scriptstyle(l)}}{{\geq}}L_{\overline{S1},\overline{12},\widetilde{2U}}T_{2}\left\|\hat{\mathbf{v}}_{1}^{H}\mathrm{diag}(\mathbf{a}_{D,12}^{H})\bar{\mathbf{H}}_{S1}\right\|_{2}^{2}+\gamma^{(0)}-2\sqrt{L_{\overline{SU}}L_{\overline{S1},\overline{1U}}}\|\hat{\mathbf{v}}_{1}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{1U}^{H})\bar{\mathbf{H}}_{S1}\|_{2}\|\mathbf{a}_{D,SU}\|_{2}$
$\displaystyle\quad-2\sqrt{L_{\overline{S2},\widetilde{2U}}L_{\overline{S1},\overline{12},\widetilde{2U}}}\|\hat{\mathbf{v}}_{1}^{H}\mathrm{diag}\left(\mathbf{a}_{D,12}^{H}\right)\bar{\mathbf{H}}_{S1}\|_{2}\|\mathbf{a}_{D,S2}\|_{2}|r_{S2,12}|$
$\displaystyle\stackrel{{\scriptstyle(m)}}{{\geq}}L_{\overline{S1},\overline{12},\widetilde{2U}}T_{S}T_{1}^{2}T_{2}-2\sqrt{L_{\overline{SU}}L_{\overline{S1},\overline{1U}}}T_{1}T_{S}-2\sqrt{L_{\overline{S2},\widetilde{2U}}L_{\overline{S1},\overline{12},\widetilde{2U}}}T_{1}T_{S}T_{2}+\gamma^{(0)}\triangleq\gamma^{(1)\star}_{lb},$
(85)

where
$\bm{\hat{}}{\phi}_{1}\triangleq\Lambda\left(-\bm{\Delta}_{\overline{S1},\overline{12}}-\angle(r_{\overline{S1},\overline{S2}}r_{\overline{S2},\overline{12}})\mathbf{1}_{T_{1}}\right)$,
$\hat{\mathbf{v}}_{1}\triangleq\left(e^{-j\hat{\phi}_{1,t}}\right)_{t\in\mathcal{T}_{1}}$,
$(l)$ is due to $\mathfrak{R}\\{x\\}\geq-|x|,x\in\mathbb{C}$ and Cauchy-Swartz
inequality, and $(m)$ is due to (80), (82), and
$\left\|\hat{\mathbf{v}}_{1}^{H}\mathrm{diag}(\mathbf{a}_{D,12}^{H})\bar{\mathbf{H}}_{S1}\right\|_{2}^{2}=T_{S}\left|\sum_{t=1}^{T_{1}}e^{j(-\Delta_{\overline{S1},\overline{12},t}-\angle(r_{\overline{S1},\overline{S2}}r_{\overline{S2},\overline{12}})+\Delta_{\overline{S1},\overline{12},t})}\right|_{2}^{2}=T_{S}T_{1}^{2}.$
Clearly, we have
$\gamma^{(1)\star}_{lb}\leq\gamma^{(1)\star}\leq\gamma^{(1)\star}_{ub}$. As
$\lim_{T_{1},T_{2}\rightarrow\infty}\frac{\gamma^{(1)\star}_{ub}}{L_{\overline{S1},\overline{12},\widetilde{2U}}T_{S}T_{1}^{2}T_{2}}=1$
and
$\lim_{T_{1},T_{2}\rightarrow\infty}\frac{\gamma^{(1)\star}_{lb}}{L_{\overline{S1},\overline{12},\widetilde{2U}}T_{S}T_{1}^{2}T_{2}}=1$,
by the squeeze theorem, we have
$\lim_{T_{1},T_{2}\rightarrow\infty}\frac{\gamma^{(1)\star}}{L_{\overline{S1},\overline{12},\widetilde{2U}}T_{S}T_{1}^{2}T_{2}}=1.$
Thus, we can show $\gamma^{(1)\star}\stackrel{{\scriptstyle
T_{1},T_{2}\rightarrow\infty}}{{\sim}}L_{\overline{S1},\overline{12},\widetilde{2U}}T_{S}T_{1}^{2}T_{2}$.

Then, consider Case 2. For all $\bm{\phi}_{2}$, we have:

$\displaystyle\gamma^{(2)}(\bm{\phi}_{2})\stackrel{{\scriptstyle(o)}}{{=}}L_{\overline{S2},\overline{2U}}\|\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{S2}\|_{2}^{2}+L_{\widetilde{S1},\overline{12},\overline{2U}}T_{S}\left\|\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}^{H}_{2U})\bar{\mathbf{H}}_{12}\right\|_{2}^{2}+\gamma^{(0)}$
$\displaystyle\quad+2\mathfrak{R}\left\\{\mathbf{v}_{2}^{H}\left(\sqrt{L_{\overline{SU}}L_{\overline{S2},\overline{2U}}}\mathrm{diag}\left(\bar{\mathbf{h}}_{2U}^{H}\right)\bar{\mathbf{H}}_{S2}\bar{\mathbf{h}}_{SU}+\sqrt{L_{\widetilde{S1},\overline{1U}}L_{\widetilde{S1},\overline{12},\overline{2U}}}T_{S}\mathrm{diag}\left(\bar{\mathbf{h}}_{2U}^{H}\right)\bar{\mathbf{H}}_{12}\bar{\mathbf{h}}_{1U}\right)\right\\}$
$\displaystyle\stackrel{{\scriptstyle(p)}}{{\leq}}L_{\overline{S2},\overline{2U}}\|\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{S2}\|_{2}^{2}+L_{\widetilde{S1},\overline{12},\overline{2U}}T_{S}\left\|\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}^{H}_{2U})\bar{\mathbf{H}}_{12}\right\|_{2}^{2}+\gamma^{(0)}$
$\displaystyle\quad+2\sqrt{L_{\overline{SU}}L_{\overline{S2},\overline{2U}}}\|\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{S2}\|_{2}\|\mathbf{a}_{D,SU}\|_{2}+2\sqrt{L_{\widetilde{S1},\overline{1U}}L_{\widetilde{S1},\overline{12},\overline{2U}}}T_{S}\|\mathbf{v}_{2}^{H}\mathrm{diag}\left(\bar{\mathbf{h}}_{2U}^{H}\right)\bar{\mathbf{H}}_{12}\|_{2}\|\mathbf{a}_{D,1U}\|_{2}$
$\displaystyle\stackrel{{\scriptstyle(q)}}{{\leq}}L_{\overline{S2},\overline{2U}}T_{S}T_{2}^{2}+L_{\widetilde{S1},\overline{12},\overline{2U}}T_{S}T_{1}T_{2}^{2}+2\sqrt{L_{\overline{SU}}L_{\overline{S2},\overline{2U}}}T_{S}T_{2}+2\sqrt{L_{\widetilde{S1},\overline{1U}}L_{\widetilde{S1},\overline{12},\overline{2U}}}T_{S}T_{1}T_{2}+\gamma^{(0)}\triangleq\gamma^{(2)\star}_{ub},$
(86)

where $(o)$ is due to (15), $(p)$ is due to
$\mathfrak{R}\\{x\\}\leq|x|,x\in\mathbb{C}$ and Cauchy-Schwartz inequality,
and $(q)$ is due to (80) and (83). On the other hand, we have:

$\displaystyle\gamma^{(2)}(\hat{\bm{\phi}}_{2})=L_{\overline{S2},\overline{2U}}\|\hat{\mathbf{v}}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{S2}\|_{2}^{2}+L_{\widetilde{S1},\overline{12},\overline{2U}}T_{S}\left\|\hat{\mathbf{v}}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}^{H}_{2U})\bar{\mathbf{H}}_{12}\right\|_{2}^{2}+\gamma^{(0)}$
$\displaystyle\quad+2\mathfrak{R}\left\\{\hat{\mathbf{v}}_{2}^{H}\left(\sqrt{L_{\overline{SU}}L_{\overline{S2},\overline{2U}}}\mathrm{diag}\left(\bar{\mathbf{h}}_{2U}^{H}\right)\bar{\mathbf{H}}_{S2}\bar{\mathbf{h}}_{SU}+\sqrt{L_{\widetilde{S1},\overline{1U}}L_{\widetilde{S1},\overline{12},\overline{2U}}}T_{S}\mathrm{diag}\left(\bar{\mathbf{h}}_{2U}^{H}\right)\bar{\mathbf{H}}_{12}\bar{\mathbf{h}}_{1U}\right)\right\\}$
$\displaystyle\stackrel{{\scriptstyle(r)}}{{\geq}}L_{\widetilde{S1},\overline{12},\overline{2U}}T_{S}\left\|\hat{\mathbf{v}}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}^{H}_{2U})\bar{\mathbf{H}}_{12}\right\|_{2}^{2}+\gamma^{(0)}$
$\displaystyle\quad-2\sqrt{L_{\overline{SU}}L_{\overline{S2},\overline{2U}}}\|\hat{\mathbf{v}}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{S2}\|_{2}\|\mathbf{a}_{D,SU}\|_{2}-2\sqrt{L_{\widetilde{S1},\overline{1U}}L_{\widetilde{S1},\overline{12},\overline{2U}}}T_{S}\|\hat{\mathbf{v}}_{2}^{H}\mathrm{diag}\left(\bar{\mathbf{h}}_{2U}^{H}\right)\bar{\mathbf{H}}_{12}\|_{2}\|\mathbf{a}_{D,1U}\|_{2}$
$\displaystyle\stackrel{{\scriptstyle(s)}}{{\geq}}L_{\overline{S1},\overline{12},\widetilde{2U}}T_{S}T_{1}T_{2}^{2}-2\sqrt{L_{\overline{SU}}L_{\overline{S2},\overline{2U}}}T_{S}T_{2}-2\sqrt{L_{\widetilde{S1},\overline{1U}}L_{\widetilde{S1},\overline{12},\overline{2U}}}T_{S}T_{1}T_{2}+\gamma^{(0)}\triangleq\gamma^{(2)\star}_{lb},$
(87)

where
$\hat{\bm{\phi}_{2}}\triangleq\Lambda(-\bm{\Delta}_{\overline{12},\overline{2U}}-\angle(r_{\overline{12},\overline{1U}})\mathbf{1}_{T_{2}})$,
$\hat{\mathbf{v}}_{2}\triangleq\left(e^{-j\hat{\phi}_{2,t}}\right)_{t\in\mathcal{T}_{2}}$,
$(r)$ is due to $\mathfrak{R}\\{x\\}\geq-|x|,x\in\mathbb{C}$ and Cauchy-
Schwartz inequality, $(s)$ is due to (80), (83), and
$\left\|\hat{\mathbf{v}}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{12}\right\|_{2}^{2}=T_{1}|\sum_{t=1}^{T_{2}}$
$e^{j(-\Delta_{\overline{12},\overline{2U},t}-\angle(r_{\overline{12},\overline{1U}})+\Delta_{\overline{12},\overline{2U},t})}|_{2}^{2}=T_{1}T_{2}^{2}.$
Clearly, we have
$\gamma^{(2)\star}_{lb}\leq\gamma^{(2)\star}\leq\gamma^{(2)\star}_{ub}$. As
$\lim_{T_{1},T_{2}\rightarrow\infty}$
$\frac{\gamma^{(2)\star}_{ub}}{L_{\overline{S1},\overline{12},\widetilde{2U}}T_{S}T_{1}T_{2}^{2}}=1$
and
$\lim_{T_{1},T_{2}\rightarrow\infty}\frac{\gamma^{(2)\star}_{lb}}{L_{\overline{S1},\overline{12},\widetilde{2U}}T_{S}T_{1}T_{2}^{2}}=1$,
by the squeeze theorem, we have $\lim_{T_{1},T_{2}\rightarrow\infty}$
$\frac{\gamma^{(2)\star}}{L_{\overline{S1},\overline{12},\widetilde{2U}}T_{S}T_{1}T_{2}^{2}}=1,$
Thus, we can show $\gamma^{(2)\star}\stackrel{{\scriptstyle
T_{1},T_{2}\rightarrow\infty}}{{\sim}}L_{\overline{S1},\overline{12},\widetilde{2U}}T_{S}T_{1}T_{2}^{2}$.

Finally, consider Case 3. For all $\bm{\phi}_{1}$ and $\bm{\phi}_{2}$, we
have:

$\displaystyle\gamma^{(3)}(\bm{\phi}_{1},\bm{\phi}_{2})\stackrel{{\scriptstyle(t)}}{{\leq}}\gamma^{(1)\star}_{ub}+\gamma^{(2)\star}_{ub}-\gamma^{(0)}+L_{\overline{S1},\overline{12},\overline{2U}}|\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{12}\mathrm{diag}(\mathbf{a}_{A,S1})\mathbf{v}_{1}^{*}|^{2}\|\mathbf{a}^{H}_{D,S1}\|_{2}^{2}$
$\displaystyle\quad+2\sqrt{L_{\overline{S1},\overline{1U}}L_{\overline{S1},\overline{12},\overline{2U}}}|\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{12}\mathrm{diag}(\mathbf{a}_{A,S1})\mathbf{v}_{1}^{*}|\|\mathbf{a}^{H}_{D,S1}\|_{2}\|\mathbf{v}_{1}^{H}\mathrm{diag}(\bar{\mathbf{h}}^{H}_{1U})\bar{\mathbf{H}}_{S1}\|_{2}$
$\displaystyle\quad+2\sqrt{L_{\overline{S2},\overline{2U}}L_{\overline{S1},\overline{12},\overline{2U}}}|\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{12}\mathrm{diag}(\mathbf{a}_{A,S1})\mathbf{v}_{1}^{*}|\|\mathbf{a}^{H}_{D,S1}\|_{2}\|\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}^{H}_{2U})\bar{\mathbf{H}}_{S2}\|_{2}$
$\displaystyle\quad+2\sqrt{L_{\overline{SU}}L_{\overline{S1},\overline{12},\overline{2U}}}|\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{12}\mathrm{diag}(\mathbf{a}_{A,S1})\mathbf{v}_{1}^{*}||r_{S1,SU}|$
$\displaystyle\quad+2\sqrt{L_{\overline{S1},\overline{1U}}L_{\overline{S2},\overline{2U}}}\|\mathbf{v}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}^{H}_{2U})\bar{\mathbf{H}}_{S2}\|_{2}\|\mathbf{v}_{1}^{H}\mathrm{diag}(\bar{\mathbf{h}}^{H}_{1U})\bar{\mathbf{H}}_{S1}\|_{2}$
$\displaystyle\stackrel{{\scriptstyle(u)}}{{\leq}}\gamma^{(1)\star}_{ub}+\gamma^{(2)\star}_{ub}-\gamma^{(0)}+L_{\overline{S1},\overline{12},\overline{2U}}T_{S}T_{1}^{2}T_{2}^{2}+2\sqrt{L_{\overline{S1},\overline{1U}}L_{\overline{S1},\overline{12},\overline{2U}}T_{S}}T_{1}^{2}T_{2}$
$\displaystyle+2\sqrt{L_{\overline{S2},\overline{2U}}L_{\overline{S1},\overline{12},\overline{2U}}T_{S}}T_{1}T_{2}^{2}+2(\sqrt{L_{\overline{SU}}L_{\overline{S1},\overline{12},\overline{2U}}}+\sqrt{L_{\overline{S1},\overline{1U}}L_{\overline{S2},\overline{2U}}})T_{S}T_{1}T_{2}\triangleq\gamma^{(3)\star}_{ub},$

where $(t)$ is due to (84), (86), $\mathfrak{R}\\{x\\}\leq|x|,x\in\mathbb{C}$,
and Cauchy-Schwartz inequality, and $(u)$ is due to (80)-(83). On the other
hand, we have:

$\displaystyle\gamma^{(3)}(\hat{\bm{\phi}}_{1},\hat{\bm{\phi}}_{2})\stackrel{{\scriptstyle(v)}}{{\geq}}\gamma^{(1)\star}_{lb}+\gamma^{(2)\star}_{lb}-\gamma^{(0)}+L_{\overline{S1},\overline{12},\overline{2U}}|\hat{\mathbf{v}}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{12}\mathrm{diag}(\mathbf{a}_{A,S1})\hat{\mathbf{v}}_{1}^{*}|^{2}\|\mathbf{a}^{H}_{D,S1}\|_{2}^{2}$
$\displaystyle\quad-2\sqrt{L_{\overline{S1},\overline{1U}}L_{\overline{S1},\overline{12},\overline{2U}}}|\hat{\mathbf{v}}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{12}\mathrm{diag}(\mathbf{a}_{A,S1})\hat{\mathbf{v}}_{1}^{*}|\|\mathbf{a}^{H}_{D,S1}\|_{2}\|\hat{\mathbf{v}}_{1}^{H}\mathrm{diag}(\bar{\mathbf{h}}^{H}_{1U})\bar{\mathbf{H}}_{S1}\|_{2}$
$\displaystyle\quad-2\sqrt{L_{\overline{S2},\overline{2U}}L_{\overline{S1},\overline{12},\overline{2U}}}|\hat{\mathbf{v}}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{12}\mathrm{diag}(\mathbf{a}_{A,S1})\hat{\mathbf{v}}^{*}|\|\mathbf{a}^{H}_{D,S1}\|_{2}\|\hat{\mathbf{v}}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}^{H}_{2U})\bar{\mathbf{H}}_{S2}\|_{2}$
$\displaystyle\quad-2\sqrt{L_{\overline{SU}}L_{\overline{S1},\overline{12},\overline{2U}}}|\hat{\mathbf{v}}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{12}\mathrm{diag}(\mathbf{a}_{A,S1})\hat{\mathbf{v}}_{1}^{*}||r_{S1,SU}|$
$\displaystyle\quad-2\sqrt{L_{\overline{S1},\overline{1U}}L_{\overline{S2},\overline{2U}}}\|\hat{\mathbf{v}}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}^{H}_{2U})\bar{\mathbf{H}}_{S2}\|_{2}\hat{\mathbf{v}}_{1}^{H}\mathrm{diag}(\bar{\mathbf{h}}^{H}_{1U})\bar{\mathbf{H}}_{S1}\|_{2}$
$\displaystyle\stackrel{{\scriptstyle(w)}}{{\geq}}\gamma^{(1)\star}_{lb}+\gamma^{(2)\star}_{lb}-\gamma^{(0)}+L_{\overline{S1},\overline{12},\overline{2U}}T_{S}T_{1}^{2}T_{2}^{2}-2\sqrt{L_{\overline{S1},\overline{1U}}L_{\overline{S1},\overline{12},\overline{2U}}T_{S}}T_{1}^{2}T_{2}$
$\displaystyle-2\sqrt{L_{\overline{S2},\overline{2U}}L_{\overline{S1},\overline{12},\overline{2U}}T_{S}}T_{1}T_{2}^{2}-2(\sqrt{L_{\overline{SU}}L_{\overline{S1},\overline{12},\overline{2U}}}+\sqrt{L_{\overline{S1},\overline{1U}}L_{\overline{S2},\overline{2U}}})T_{S}T_{1}T_{2}\triangleq\gamma^{(3)\star}_{lb},$

where $(v)$ is due to (85), (87),
$\mathfrak{R}\\{x\\}\geq-|x|,x\in\mathbb{C}$, and Cauchy-Schwartz inequality,
and $(w)$ is due to (80)- (83), and
$\left|\hat{\mathbf{v}}_{2}^{H}\mathrm{diag}(\bar{\mathbf{h}}_{2U}^{H})\bar{\mathbf{H}}_{12}\mathrm{diag}(\mathbf{a}_{A,S1})\hat{\mathbf{v}}_{1}^{*}\right|=\left|\sum_{t=1}^{T_{1}}e^{j(-\Delta_{\overline{S1},\overline{12},t}-\angle(r_{\overline{S1},\overline{S2}}r_{\overline{S2},\overline{12}})+\Delta_{\overline{S1},\overline{12},t})}\right|^{2}$
$\left|\sum_{t=1}^{T_{2}}e^{j(-\Delta_{\overline{12},\overline{2U},t}-\angle(r_{\overline{12},\overline{1U}})+\Delta_{\overline{12},\overline{2U},t})}\right|^{2}=T_{1}^{2}T_{2}^{2}.$
Clearly, we have
$\gamma^{(3)\star}_{lb}\leq\gamma^{(3)\star}\leq\gamma^{(3)\star}_{ub}$. As
$\lim_{T_{1},T_{2}\rightarrow\infty}$
$\frac{\gamma^{(3)\star}_{ub}}{L_{\overline{S1},\overline{12},\overline{2U}}T_{S}T_{1}^{2}T_{2}^{2}}=1$
and
$\lim_{T_{1},T_{2}\rightarrow\infty}\frac{\gamma^{(3)\star}_{lb}}{L_{\overline{S1},\overline{12},\overline{2U}}T_{S}T_{1}^{2}T_{2}^{2}}=1$,
by the squeeze theorem, we have
$\lim_{T_{1},T_{2}\rightarrow\infty}\frac{\gamma^{(3)\star}}{L_{\overline{S1},\overline{12},\overline{2U}}T_{S}T_{1}^{2}T_{2}^{2}}=1,$
Thus, we can show $\gamma^{(3)\star}\stackrel{{\scriptstyle
T_{1},T_{2}\rightarrow\infty}}{{\sim}}L_{\overline{S1},\overline{12},\overline{2U}}T_{S}T_{1}^{2}T_{2}^{2}$.

Therefore, we complete the proof of Statement (i). Furthermore, by
substituting $T_{1}=cT$ and $T_{2}=(1-c)T$ into Statement (i), we can easily
show Statement (ii).

## Appendix H: Proof of Lemma 4

First, we prove Statement (i). By (17), we have
$\bar{\gamma}^{(3)}(\bm{\phi}_{1},\bm{\phi}_{2})=\mathbf{v}_{1}^{H}\bar{\mathbf{G}}_{1}^{H}\bar{\mathbf{G}}_{1}\mathbf{v}_{1}+2\mathfrak{R}\\{\bar{\mathbf{G}}_{1}^{H}\bar{\mathbf{g}}_{1}\\}+\bar{\mathbf{g}}_{1}^{H}\bar{\mathbf{g}}_{1}$
where
$\bar{\mathbf{G}}_{1}\triangleq\bar{\mathbf{H}}_{S1}^{H}\mathrm{diag}(\sqrt{\alpha_{S1}\alpha_{1U}}\bar{\mathbf{h}}_{S1}+\sqrt{\alpha_{S1}\alpha_{12}\alpha_{2U}}\mathbf{B}_{1}^{H}\mathbf{v}_{2})$
and
$\bar{\mathbf{g}}_{1}^{H}\triangleq\sqrt{\alpha_{SU}}\bar{\mathbf{h}}_{SU}^{H}+\sqrt{\alpha_{S2}\alpha_{2U}}\bar{\mathbf{h}}_{2U}^{H}$
$\mathrm{diag}(\mathbf{v}_{2}^{H})\bar{\mathbf{H}}_{S2}$. By Lemma 13, we can
show that the unique optimal solution of Problem 6 is

$\displaystyle\bar{\bm{\phi}}_{1}^{\dagger}$
$\displaystyle=\Lambda\Big{(}-\angle\left(\mathrm{diag}\left(\sqrt{\alpha_{S1}\alpha_{1U}}\bar{\mathbf{h}}_{1U}^{H}+\sqrt{\alpha_{S1}\alpha_{1U}\alpha_{S1}\alpha_{12}\alpha_{2U}}\mathbf{v}_{2}^{H}\mathbf{B}_{1}\right)\mathbf{a}_{A,S1}\right)$
$\displaystyle\quad\quad\ \
-\angle\big{(}\mathbf{a}^{H}_{D,S1}(\sqrt{\alpha_{SU}}\bar{\mathbf{h}}_{SU}+\sqrt{\alpha_{S2}\alpha_{2U}}\bar{\mathbf{H}}_{S2}^{H}\mathrm{diag}(\mathbf{v}_{2})\bar{\mathbf{h}}_{2U})\big{)}\mathbf{1}_{T_{1}}\Big{)}.$

Next, we prove Statement (ii). By (16), we have:

$\displaystyle\bar{\gamma}^{(3)}(\phi_{2,t},\bm{\phi}_{2,-t},\bm{\phi}_{1})\overset{(a)}{=}$
$\displaystyle
2\mathfrak{R}\left\\{e^{j\phi_{2,t}}\left(\sum_{k\in\mathcal{T}_{2},k\neq
t}\bar{C}_{2,t,k}e^{-j\phi_{2,k}}+\bar{d}_{2,t}\right)\right\\}+\sum_{k\in\mathcal{T}_{2},k\neq
t}\sum_{i\in\mathcal{T}_{2},i\neq
t}\bar{C}_{2,k,i}e^{j(\phi_{2,k}-\phi_{2,i})}$
$\displaystyle+\bar{C}_{2,t,t}+2\mathfrak{R}\left\\{\sum_{k\in\mathcal{T}_{2},k\neq
t}e^{j\phi_{2,k}}\bar{d}_{2,k}\right\\}+\gamma^{(0)},$

where $(a)$ is due to $\bar{\mathbf{C}}_{l}=\bar{\mathbf{C}}_{l}^{H}$. As
$\left(\bm{\phi}_{2,-t},\bm{\phi}_{1}\right)$ is given, to solve Problem 7 is
equivalent to solve:

$\displaystyle\max_{\phi_{2,t}\in[0,2\pi)}\mathfrak{R}\left\\{e^{j\phi_{l,t}}\left(\sum_{k\in\mathcal{T}_{2},k\neq
t}\bar{C}_{2,t,k}e^{-j\phi_{2,k}}+\bar{d}_{2,t}\right)\right\\}.$

Besides, we have:

$\mathfrak{R}\left\\{e^{j\phi_{2,t}}\left(\sum_{k\in\mathcal{T}_{2},k\neq
t}\bar{C}_{2,t,k}e^{-j\phi_{2,k}}+\bar{d}_{2,t}\right)\right\\}\overset{(b)}{\leq}\left|\sum_{k\in\mathcal{T}_{2},k\neq
t}\bar{C}_{2,t,k}e^{-j\phi_{2,k}}+\bar{d}_{2,t}\right|,$

where $(b)$ is due to $\mathfrak{R}\\{x\\}\leq|x|,x\in\mathbb{C}$ and holds
with equality if and only if
$\phi_{2,t}=-\angle\big{(}\sum_{k\in\mathcal{T}_{2},k\neq t}$
$\bar{C}_{2,t,k}e^{-j\phi_{2,k}}+\bar{d}_{2,t}\big{)}$. By noting that
$\phi_{2,t}\in[0,2\pi)$, the unique optimal solution of the equivalent problem
is
$\bar{\phi}_{2,t}^{\dagger}=\Lambda\left(-\angle\left(\sum_{k\in\mathcal{T}_{2},k\neq
t}\bar{C}_{2,t,k}e^{-j\phi_{2,k}}+\bar{d}_{2,t}\right)\right)$. Therefore, we
complete the proof of Lemma 4.

## References

  * [1] G. Ding, Y. Cui, L. Hu, F. Yang, L. Ding, and X. X. Chen, “Quasi-Static phase shift design for a Double-IRS cooperatively assisted system,” in _Proc. of IEEE SPAWC 2022_ , Jul., 2022.
  * [2] Q. Wu and R. Zhang, “Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network,” _IEEE Commun. Mag._ , vol. 58, no. 1, pp. 106–112, Nov., 2020.
  * [3] 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, Nov., 2019.
  * [4] H. Guo, Y. Liang, J. Chen, and E. G. Larsson, “Weighted sum-rate maximization for reconfigurable intelligent surface aided wireless networks,” _IEEE Trans. Commun._ , vol. 19, no. 5, pp. 3064–3076, May., 2020.
  * [5] X. Yu, D. Xu, and R. Schober, “Enabling secure wireless communications via intelligent reflecting surfaces,” in _Proc. of IEEE GLOBECOM 2019_ , Feb., 2019.
  * [6] C. Huang, A. Zappone, G. C. Alexandropoulos, M. Debbah, and C. Yuen, “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 8, pp. 4157–4170, Aug., 2019.
  * [7] X. Yu, D. Xu, Y. Sun, D. W. K. Ng, and R. Schober, “Robust and secure wireless communications via intelligent reflecting surfaces,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 11, pp. 2637–2652, Nov., 2020\.
  * [8] P. Wang, J. Fang, X. Yuan, Z. Chen, and H. Li, “Intelligent reflecting surface-assisted millimeter wave communications: Joint active and passive precoding design,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 12, pp. 14 960–14 973, Dec., 2020.
  * [9] Y. Han, W. Tang, S. Jin, C. K. Wen, and X. Ma, “Large intelligent surface-assisted wireless communication exploiting statistical CSI,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 8, pp. 8238–8242, Aug., 2019\.
  * [10] Y. Jia, C. Ye, and Y. Cui, “Analysis and optimization of an intelligent reflecting surface-assisted system with interference,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 12, pp. 8068–8082, Dec., 2020.
  * [11] Y. Jia, Y. Cui, and W. Jiang, “Robust optimization of instantaneous beamforming and quasi-static phase shifts in an IRS-assisted multi-cell network,” _to appear in IEEE Trans. Wireless Commun._ , Dec. 2021.
  * [12] C. Guo, Y. Cui, F. Yang, and L. Ding, “Outage probability analysis and minimization in intelligent reflecting surface-assisted MISO systems,” _IEEE Commun. Lett._ , vol. 24, no. 7, pp. 1563–1567, Feb., 2020.
  * [13] M. Zhao, A. Liu, Y. Wan, and R. Zhang, “Two-timescale beamforming optimization for intelligent reflecting surface aided multiuser communication with QoS constraints,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 9, pp. 6179–6194, Apr., 2021.
  * [14] Z. Zhang, Y. Cui, F. Yang, and L. Ding, “Analysis and optimization of outage probability in multi-intelligent reflecting surface-assisted systems,” Sept., 2019. [Online]. Available: https://arxiv.org/abs/1909.02193
  * [15] Y. Gao, J. Xu, W. Xu, D. W. K. Ng, and M. Alouini, “Distributed IRS with statistical passive beamforming for MISO communications,” _IEEE Commun. Lett._ , vol. 10, no. 2, pp. 221–225, Feb., 2021.
  * [16] C. You, B. Zheng, and R. Zhang, “Wireless communication via double IRS: Channel estimation and passive beamforming designs,” _IEEE Commun. Lett._ , vol. 10, no. 2, pp. 431–435, Feb., 2021.
  * [17] B. Zheng, C. You, and R. Zhang, “Efficient channel estimation for Double-IRS aided multi-user MIMO system,” _IEEE Trans. Commun._ , vol. 69, no. 6, pp. 3818–3832, Jun., 2021.
  * [18] Y. Han, S. Zhang, L. Duan, and R. Zhang, “Cooperative Double-IRS aided communication: Beamforming design and power scaling,” _IEEE Commun. Lett._ , vol. 9, no. 8, pp. 1206–1210, Aug., 2020.
  * [19] W. Mei, B. Zheng, C. You, and R. Zhang, “Intelligent reflecting surface aided wireless networks: From single-reflection to multi-reflection design and optimization,” Sept., 2021. [Online]. Available: https://arxiv.org/abs/2109.13641
  * [20] B. Zheng, C. You, and R. Zhang, “Double-IRS assisted multi-user MIMO: Cooperative passive beamforming design,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 7, pp. 4513–4526, Feb., 2021.
  * [21] Y. Han, S. Zhang, L. Duan, and R. Zhang, “Double-IRS aided MIMO communication under LoS channels: Capacity maximization and scaling,” _IEEE Trans. Commun._ , vol. 70, no. 4, pp. 2820–2837, Feb., 2022.
  * [22] L. Dong, H. Wang, J. Bai, and H. Xiao, “Double intelligent reflecting surface for secure transmission with inter-surface signal reflection,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 3, pp. 2912–2916, Mar., 2021\.
  * [23] W. Mei and R. Zhang, “Multi-beam multi-hop routing for intelligent reflecting surfaces aided massive MIMO,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 3, pp. 1897–1912, Mar., 2022.
  * [24] S. He and X. Dong, “High-accuracy localization platform using asynchronous time difference of arrival technology,” _IEEE Trans. Instrum. Meas._ , vol. 66, no. 7, pp. 1728–1742, Jul., 2017.
  * [25] S. Zeng, H. Zhang, B. Di, Z. Han, and L. Song, “Reconfigurable intelligent surface (RIS) assisted wireless coverage extension: RIS orientation and location optimization,” _IEEE Commun. Lett._ , vol. 25, no. 1, pp. 269–273, Jan., 2021.
  * [26] M. Jiang, M. Hosseinian, M.-i. Lee, and J. Stern-Berkowitz, “3D channel model extensions and characteristics study for future wireless systems,” in _Proc. of IEEE PIMRC 2013_ , Sept., 2013.
  * [27] T. S. Rappaport, S. Sun, R. Mayzus, H. Zhao, Y. Azar, K. Wang, G. N. Wong, J. K. Schulz, M. Samimi, and F. Gutierrez, “Millimeter wave mobile communications for 5G cellular: It will work!” _IEEE Access_ , vol. 1, pp. 335–349, May., 2013.
  * [28] X. Hu, C. Zhong, Y. Zhang, X. Chen, and Z. Zhang, “Location information aided multiple intelligent reflecting surface systems,” _IEEE Trans. Commun._ , vol. 68, no. 12, pp. 7948–7962, Dec., 2020.
  * [29] M. Biguesh and A. Gershman, “Training-based mimo channel estimation: a study of estimator tradeoffs and optimal training signals,” _IEEE Trans. Signal Process._ , vol. 54, no. 3, pp. 884–893, Mar., 2006.
  * [30] ETSI, “ETSI TR 138 901 v16.1.0: Study on channel model for frequencies from 0.5 to 100 GHz,” Nov., 2020. [Online]. Available: https://www.etsi.org/deliver/etsitr/138900138999/138901/16.01.0060/tr138901v160100p.pdf
  * [31] D. P. Bertsekas, _Nonlinear progranmming_. Athena scientific Belmont, MA, 1998.
  * [32] K. Zhi, C. Pan, H. Ren, and K. Wang, “Statistical CSI-Based design for reconfigurable intelligent surface-aided massive MIMO systems with direct links,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 5, pp. 1128–1132, May., 2021.
  * [33] J. R. Schott, _Matrix Analysis for Statistics_. Wiley, New York, 1997.
  * [34] D. Tse and P. Viswanath, _Fundamentals of Wireless Communication_. Cambridge University Press, 2005.