Characteristics of Channel Eigenvalues and Mutual Coupling Effects for Holographic Reconfigurable Intelligent Surfaces

Abstract

As a prospective key technology for the next-generation wireless communications, reconfigurable intelligent surfaces (RISs) have gained tremendous research interest in both the academia and industry in recent years. Only limited knowledge, however, has been obtained about the channel eigenvalue characteristics and spatial degrees of freedom (DoF) of systems containing RISs, especially when mutual coupling (MC) is present between the array elements. In this paper, we focus on the small-scale spatial correlation and eigenvalue properties excluding and including MC effects, for RISs with a quasi-continuous aperture (i.e., holographic RISs). Specifically, asymptotic behaviors of far-field and near-field eigenvalues of the spatial correlation matrix of holographic RISs without MC are first investigated, where the counter-intuitive observation of a lower DoF with more elements is explained by leveraging the power spectrum of the spatial correlation function. Second, a novel metric is proposed to quantify the inter-element correlation or coupling strength in RISs and ordinary antenna arrays. Furthermore, in-depth analysis is performed regarding the MC effects on array gain, effective spatial correlation, and eigenvalue architectures for a variety of element intervals when a holographic RIS works in the radiation and reception mode, respectively. The analysis and numerical results demonstrate that a considerable amount of the eigenvalues of the spatial correlation matrix correspond to evanescent waves that are promising for near-field communication and sensing. More importantly, holographic RISs can potentially reach an array gain conspicuously larger than conventional arrays by exploiting MC, and MC has discrepant impacts on the effective spatial correlation and eigenvalue structures at the transmitter and receiver.

1. Introduction

The sixth-generation (6G) communication networks is envisioned to embrace numerous new use cases and challenging requirements [1]. Among the emerging candidate physical-layer technologies for 6G, reconfigurable intelligent surfaces (RISs), sometimes also named large intelligent surfaces [2] and holographic multiple-input multiple-output (MIMO) [3], shows promising foreground in capacity and coverage enhancement, reconfigurable environment construction, intelligent sensing and control, and holographic communications [2,3,4,5,6,7,8,9]. We utilize holographic RIS [7,10] herein as an umbrella term for the two-dimensional (2D) architectures with an element spacing equal to or smaller than half a wavelength of the carrier frequency, which can be perceived as an extension of massive MIMO [11] with the ultimate form of (approximately) spatially-continuous electromagnetic (EM) aperture [12]. Holographic RISs can be employed at the base station (BS), user equipment (UE), and/or interacting objects in the propagation medium, and is likely to bring immense advantages in not only spectral efficiency, energy efficiency, and system scalability inherited from massive MIMO [13,14,15], but also the manipulation of EM waves via anomalous reflection, refraction, polarization transformation, and so on [3], reminiscent of metasurfaces in the optical regime [16,17]. In order to unleash the full potentials of holographic RISs, it is necessary to understand its fundamental properties, such as channel eigenvalues and spatial degrees of freedom (DoF).

1.1. Related Work

The use of dense antenna arrays for wireless communications was in fact explored over two decades ago (e.g., [18,19,20]), and has resurged recently as a promising 6G technology [2,3,12]. Due to the limited element spacing and 2D (as opposed to one-dimensional) structure of a dense array, the spatial correlation between array elements is not always non-zero even under isotropic scattering [12,21,22]. Among the early work involving antenna spatial correlation, the authors in [23] have considered spatial correlation among multi-element antennas, and derived upper and lower capacity bounds taking into account antenna correlation. Fading correlation has also been studied in [24] to examine the capacity growth with respect to the number of antenna elements. The impact of antenna correlation on capacity has been examined for diverse signal inputs and correlation architectures in [25]. Nevertheless, the antenna element spacing is of half-wavelength or larger in [23,24,25]. The capacity of spatially dense multiple antenna systems has been pioneered in [18], which showed that the capacity of such a system approaches a finite limit. Spatially dense MIMO arrays have also been studied in [19], where the array-gain normalized capacity has been analyzed and the performance of small wavelength-like MIMO arrays has been shown to be similar to that of arrays with larger apertures. The asymptotic capacity associated with antenna arrays of fixed length has been analyzed in [20] for uniform linear antenna arrays, revealing that the asymptotic mutual information converges almost surely as the number of antenna elements approaches infinity due to the convergence of the eigenvalues. The aforementioned work, however, focused mainly on the capacity and did not explicitly consider the spatial DoF that the arrays can offer.

The spatial DoF for holographic RISs in line-of-sight (LoS) environments have been investigated in [26], which has revealed that the DoF can be larger than one even in strong LoS channel conditions, favorable for spatial multiplexing. In [2,12], the asymptotic spatial DoF for sufficiently dense and large holographic RISs have been derived from the perspectives of channel capacity and Fourier plane-wave series expansion, respectively. The achievable DoF for more common cases with finite element spacing and aperture areas has been investigated in [22], where it has been discovered that the spatial DoF decreases as the number of elements grows which seems counter-intuitive, but the underlying causes have not been identified.

Due to the close proximity of neighboring elements in a holographic RIS, mutual coupling (MC) naturally arises. Broadly speaking, MC refers to EM interaction among the array elements and can occur because of three mechanisms: direct space coupling between array elements, indirect coupling caused by near-by scatterers, and coupling through feed network [27,28]. In this paper, we mainly refer to MC stemming from the first mechanism. MC between array elements can be characterized by a conventional multi-port circuit model, such as an impedance matrix, an admittance matrix, or a scattering matrix. Theoretically, the type of matrix used to represent the array network is not important since matrix transforms can be applied to change the type of matrix representation. There is abundant early research work on MC in antenna arrays (e.g., [28,29,30,31,32,33,34,35,36] and references therein), which considered linear arrays only or did not fix the array aperture while varying the element spacing. Aiming at exploring the physical characteristics of the emerging holographic RISs or similar structures for 6G communications, the authors in [37] have proposed an LoS communication model incorporating MC in the form of mutual impedance, and studied the associated optimal beamforming strategies. The model has then been extended in [38] to account for superdirectivity and MC effects as well as near-field propagation. An end-to-end MC-aware communication model based on mutual impedances has been propounded in [39], which is also EM-compliant and unit cell aware. In [40], a circuit-based MIMO channel model considering antenna MC, size-related antenna equivalent circuit, and channel small-scale fading has been proposed, and statistical properties including temporal autocorrelation function and spatial cross-correlation function, along with the effects of MC and size-related antenna equivalent circuit on them have been revealed. The authors in [41] have derived a beamforming vector of superdirective arrays based on a coupling matrix-enabled method, and proposed an approach to obtain the coupling matrix via spherical wave expansion.

1.2. Contributions

Despite the aforementioned extensive research work, study on small-scale spatial characteristics, the associated MC effects, and the mechanism behind some unique phenomena for holographic RISs are still in the infancy. In this paper, therefore, we carry out thorough investigation on the aspects above. Specifically, the small-scale spatial correlation, eigenvalue behavior, and spatial DoF for holographic RISs excluding and including MC are explored. The major novelty and contributions of this article lie in the following aspects:

• First, leveraging the block-Toeplitz with Toeplitz block (BTTB) matrix theory, we relate the eigenvalues of the spatial correlation matrix of the holographic RIS to the power spectrum of the spatial correlation function, and explain the counter-intuitive phenomenon of seemingly lower spatial DoF with growing numbers of elements in a holographic RIS observed in our prior work [22], which has not been addressed in the literature to our best knowledge. This analysis also helps with distinguishing the spatial DoF corresponding to the far field and near field of a holographic RIS.
• Second, we incorporate MC into the array response and spatial correlation matrix of the holographic RIS considering realistic element sizes, and demonstrate the potential of holographic RISs to reach an extraordinary array gain that is significantly higher than conventional antenna arrays with concrete examples. The results indicate that, different from the common belief that MC is always deleterious and should be avoided or compensated for, MC can be beneficial in boosting the array gain of holographic RISs even without sophisticated manipulation of excitation coefficients for the array elements.
• Furthermore, in-depth analysis and comparisons are performed regarding the MC effects on spatial correlation and the corresponding eigenvalue distributions for holographic RISs working in the transmitting (Tx) and receiving (Rx) modes, respectively, and with various element intervals as well as source and load impedance values. A metric named inter-element correlation/coupling strength indicator (ICSI) is proposed to measure the amount of inter-element correlation/coupling within an array. Results show that the effects of MC are quite discrepant for Tx and Rx arrays, and are also dependent upon element spacing, source and load impedance, among other factors, necessitating comprehensive design and implementation considerations.

Isotropic scattering is considered in this paper, since it is a typical type of environment involved in an enormous amount of theoretical research work, and also encountered by low-frequency bands (e.g., sub-1 GHz) which are still crucial even in 6G to guarantee wide coverage and high reliability [42], and the corresponding results can serve as theoretical upper bounds for non-isotropic scattering scenarios. The contributions in this paper unveil some fundamental channel eigenvalue features and practical spatial DoF of holographic RISs, and provide valuable hints on channel estimation and beamforming strategies for holographic RISs [38,41,43].

1.3. Article Outline and Notation

The remainder of this paper is organized as follows: in Section 2, we describe the system model, and formulate array responses and spatial correlation excluding and including MC for holographic RISs. Analyses of asymptotic eigenvalue distributions of the spatial correlation matrix without and with MC are provided in Section 3 and Section 4, respectively. Conclusions are drawn in Section 5.

The following notations will be utilized throughout the paper: $\mathbf{A}$ for matrix, $\mathbf{a}$ for column vector, ${\left[\mathbf{A}\right]}_{i,j}$ for the $(i,j)$th entry of $\mathbf{A}$; ${\mathbf{A}}^{\mathrm{T}}$, ${\mathbf{A}}^{*}$, and ${\mathbf{A}}^{\mathrm{H}}$ for the transpose, conjugate, and Hermitian of $\mathbf{A}$, respectively; $det\left(\mathbf{A}\right)$ for determinant of the square matrix $\mathbf{A}$, $\mathrm{tr}\left(\mathbf{A}\right)$ for the trace of $\mathbf{A}$, ${\mathbf{I}}_N$ for the $N\times N$ identity matrix, while $\lceil a\rceil$ and $\lfloor a\rfloor$ for the ceiling and floor of the scalar a, respectively.

2. System Model

We consider a holographic RIS in a wireless communication system where the holographic RIS can be employed at the BS, UE, and/or interacting objects that can radiate, receive, reflect, or refract wireless signals, as illustrated in Figure 1. The horizontal and vertical lengths of the holographic RIS are $L_x$ and $L_z$, with element spacing of $d_x$ and $d_z$, respectively. Each element in the holographic RIS is modeled as a cylindrical thin wire of perfectly conducting material, and is connected to a tunable load, where the load can be a positive-intrinsic-negative diode whose inductance and capacitance are adaptable to reconfigure the response of each element [39].

The signal sent from or impinging on the holographic RIS is generally composed of a superposition of multipath components which can be regarded as a continuum of plane waves, hence the channel capturing small-scale fading can be expressed as [44]

$$\mathbf{h}={\int }_0^{\pi }{\int }_0^{\pi }\tilde{s}(\varphi ,\theta )\mathbf{a}(\varphi ,\theta )d\varphi d\theta$$

(1)

where $\tilde{s}(\varphi ,\theta )$ denotes the angular distribution function that contains the channel gain and phase shift corresponding to the direction $(\varphi ,\theta )$ with $\varphi$ and $\theta$ representing the azimuth and zenith angles, respectively, while $\mathbf{a}(\varphi ,\theta )$ is the array response vector. The correlation function is given by

$$\mathbf{R}=\mathbb{E}\left\{\mathbf{h}{\mathbf{h}}^{\mathrm{H}}\right\}={\int }_0^{\pi }{\int }_0^{\pi }s(\varphi ,\theta )\mathbf{a}(\varphi ,\theta ){\mathbf{a}}^{\mathrm{H}}(\varphi ,\theta )d\varphi d\theta$$

(2)

where $s(\varphi ,\theta )$ denotes the normalized spatial scattering function satisfying

$$\mathbb{E}\left\{\tilde{s}(\varphi ,\theta ){\tilde{s}}^{*}({\varphi }^{\prime },{\theta }^{\prime })\right\}=s(\varphi ,\theta )\delta (\varphi -{\varphi }^{\prime })\delta (\theta -{\theta }^{\prime })$$

(3)

and ${\int }_0^{\pi }{\int }_0^{\pi }s(\varphi ,\theta )d\varphi d\theta =1$. Physically, the presence of spatial correlation means that the signal strengths at different elements do not vary independently, but may rise or fade simultaneously.

2.1. Array Response and Spatial Correlation Excluding MC

For an array with N elements, where each element has the same pattern function of $p(\varphi ,\theta )$ (Strictly speaking, if MC exists, the central elements and the ones near the array edges may not maintain the same element pattern when embedded in an array, even if their isolated element patterns are identical [27]. Nevertheless, it is possible to compensate for the pattern distortion via predetermined illumination, and here we assume the same embedded element pattern for all the elements in an array. The element pattern variation owing to MC is another topic and is deferred to future work), the far-field radiation pattern of the array is expressed as

$$f(\varphi ,\theta )=\sum _{n=1}^N{w}_{n}p(\varphi ,\theta )e^{j\kappa \widehat{\mathbf{d}}·{\mathbf{d}}_n}$$

(4)

where $w_n$ denotes the complex excitation coefficient proportional to the current on the n-th element, $\kappa =2\pi /\lambda$ is the wavenumber with $\lambda$ being the carrier wavelength, $\widehat{\mathbf{d}}$ represents the unit vector of the far-field direction $(\varphi ,\theta )$ in the spherical coordinate system, and ${\mathbf{d}}_n$ is the position vector of the n-th element. (4) holds if there is no MC among the array elements, which is usually the case when the spacing between adjacent elements is sufficiently large (e.g., more than a couple of wavelengths). Accordingly, the conventional MC-unaware array response vector is defined as

$${\mathbf{a}}_0={\left[e^{j\kappa \widehat{\mathbf{d}}·{\mathbf{d}}_1},\dots ,e^{j\kappa \widehat{\mathbf{d}}·{\mathbf{d}}_n},\dots ,e^{j\kappa \widehat{\mathbf{d}}·{\mathbf{d}}_N}\right]}^{\mathrm{T}}.$$

(5)

Then, the correlation matrix ${\mathbf{R}}_0$ excluding the MC effects is given by

$${\mathbf{R}}_0={\int }_0^{\pi }{\int }_0^{\pi }s(\varphi ,\theta ){\mathbf{a}}_0(\varphi ,\theta ){\mathbf{a}}_0^{\mathrm{H}}(\varphi ,\theta )d\varphi d\theta .$$

(6)

The concrete expression and properties of ${\mathbf{R}}_0$ will be provided in Section 3 to investigate the characteristics of its eigenvalues.

2.2. Array Response and Spatial Correlation Including MC

In a holographic RIS, the elements are densely arranged so that MC is usually non-negligible. The element pattern considering MC can be formulated as [32]

$$\tilde{p}(\varphi ,\theta )=p(\varphi ,\theta )\sum _{m=1}^N{c}_{mn}{e}^{j\kappa \widehat{\mathbf{d}}·{\mathbf{d}}_m}$$

(7)

where $c_{mn}$ represents the coupling coefficient between the m-th and the n-th elements. Namely, the pattern of the n-th element can be modeled as the original element pattern $p(\varphi ,\theta )$ times a weighted sum of impact from all the other elements. Thus, the radiation pattern of the array incorporating MC follows

$$\tilde{f}(\varphi ,\theta )=\sum _{n=1}^N\sum _{m=1}^N{w}_{n}p(\varphi ,\theta )c_{mn}{e}^{j\kappa \widehat{\mathbf{d}}·{\mathbf{d}}_m}.$$

(8)

The coupling matrix collecting the MC coefficients is written as

$$\mathbf{C}=\left[\begin{array}{cccc}{c}_{11}& c_{12}& \cdots & c_{1N}\\ c_{21}& c_{22}& \cdots & c_{2N}\\ ⋮& ⋮& \ddots & ⋮\\ c_{N1}& c_{N2}& \cdots & c_{NN}\end{array}\right].$$

(9)

When there is no MC among the elements, $\mathbf{C}$ reverts to an identity matrix. It can be derived from (5), (7) and (9) that the effective array response vector $\mathbf{a}$ involving MC is given by

$$\mathbf{a}={\mathbf{C}}^{\mathrm{T}}{\mathbf{a}}_0$$

(10)

in which ${\mathbf{a}}_0$ stands for the original MC-unaware array response vector in (5).

Now, let us look at a crucial metric relevant to the array response—array gain, which is defined herein as the increase in radiation power of an array compared with that of a single element under the same total excitation power. It is well known in the antenna literature [45,46,47] that the array gain of an array with closely-spaced elements can grow with the square of the number of elements, and rigorous proof for a linear array was provided in [46] with optimal beamforming in the end-fire direction. For an arbitrary pointing direction $(\varphi ,\theta )$, the array gain $G_{\mathrm{array}}$ incorporating MC can be formulated based on (8) as

$$\begin{array}{cc}\hfill G_{\mathrm{array}}(\varphi ,\theta )& \hfill =\frac{|{\sum }_{n=1}^N{\sum }_{m=1}^N{w}_{n}p(\varphi ,\theta )c_{mn}{e}^{j\kappa \widehat{\mathbf{d}}·{\mathbf{d}}_m}{|}^2}{|w_{0}p{(\varphi ,\theta )|}^2}\hfill \\ \hfill & \hfill =\frac{|{\sum }_{n=1}^N{\sum }_{m=1}^N{w}_n{c}_{mn}{e}^{j\kappa \widehat{\mathbf{d}}·{\mathbf{d}}_m}{|}^2}{|w_0{|}^2}\hfill \\ \hfill & \hfill =\frac{|{\mathbf{a}}^{\mathrm{T}}\mathbf{w}{|}^2}{|w_0{|}^2}{,\mathrm{s}.\mathrm{t}.\left|\right|\mathbf{w}\left|\right|}_2=\left|w_0\right|\hfill \end{array}$$

(11)

where $w_0$ stands for the complex excitation coefficient for a single element such that $|w_0{|}^2$ represents the total excitation power, $\mathbf{a}$ is given in (10), and $\mathbf{w}$ denotes the beamforming column vector consisting of the complex excitation coefficients $w_n$ for an array. Note that $\mathbf{w}$ is inherently a function of the pointing direction $(\varphi ,\theta )$ due to $\mathbf{a}$. Consequently, the optimal beamforming vector maximizing $G_{\mathrm{array}}(\varphi ,\theta )$ is

$$\begin{array}{c}\hfill {\mathbf{w}}_{\mathrm{opt}}=\zeta {\mathbf{a}}^{*}=\zeta {\mathbf{C}}^{\mathrm{H}}{\mathbf{a}}_0^{*}\end{array}$$

(12)

in which $\zeta$ is a normalization factor equal to $\frac{|w_0|}{\left|\right|{\mathbf{C}}^{\mathrm{H}}{\mathbf{a}}_0^{*}{\left|\right|}_2}$ to satisfy the power constraint. Plugging (12) into (11) produces the maximum array gain at the direction $(\varphi ,\theta )$

$$\begin{array}{c}\hfill G_{\mathrm{array}}{(\varphi ,\theta )}_{\max }=\left|{\mathbf{a}}_0^{\mathrm{T}}\mathbf{C}{\mathbf{C}}^{\mathrm{H}}{\mathbf{a}}_0^{*}\right|.\end{array}$$

(13)

It is straightforward to observe from (13) that, when excluding MC, i.e., when $\mathbf{C}$ is an identity matrix, the maximum array gain is $\left|{\mathbf{a}}_0^{\mathrm{T}}{\mathbf{a}}_0^{*}\right|=N$ for an arbitrary direction. When taking MC into account, the maximum array gain becomes $\left|{\mathbf{a}}_0^{\mathrm{T}}\mathbf{C}{\mathbf{C}}^{\mathrm{H}}{\mathbf{a}}_0^{*}\right|$ which is usually larger than N as will be shown later by simulations, and varies with pointing angles. The array gain in (13) of a holographic RIS is highly promising, since it indicates that, compared with a traditional array with N elements, the received signal-to-noise ratio (SNR) can be enhanced by up to $\frac{\left|{\mathbf{a}}_0^{\mathrm{T}}\mathbf{C}{\mathbf{C}}^{\mathrm{H}}{\mathbf{a}}_0^{*}\right|}{N}$ fold for a fixed transmit power using a holographic RIS with the same number of elements; or, equivalently, the transmit power can be scaled down by up to a factor of $\frac{N}{\left|{\mathbf{a}}_0^{\mathrm{T}}\mathbf{C}{\mathbf{C}}^{\mathrm{H}}{\mathbf{a}}_0^{*}\right|}$ without compromising the received SNR. The disadvantages of the proposed method in (12) are that the coupling matrix $\mathbf{C}$ needs to be known (e.g., through rigorous theoretical analysis or measurements) before conducting the beamforming, and that the achievable array gain in (13) might be smaller than N in some corner cases (as will be shown later via simulations) depending on the properties of the coupling matrix $\mathbf{C}$ and target pointing angles.

Based upon the array response vector incorporating MC in (10), the effective spatial correlation matrix $\mathbf{R}$ in (2) is expanded as

$$\begin{array}{cc}\hfill \mathbf{R}& \hfill =\mathbb{E}\left\{\mathbf{h}{\mathbf{h}}^{\mathrm{H}}\right\}\hfill \\ \hfill & \hfill ={\int }_0^{\pi }{\int }_0^{\pi }s(\varphi ,\theta ){\mathbf{C}}^{\mathrm{T}}{\mathbf{a}}_0{\mathbf{a}}_0^{\mathrm{H}}(\varphi ,\theta ){\mathbf{C}}^{*}d\varphi d\theta \hfill \\ \hfill & \hfill ={\mathbf{C}}^{\mathrm{T}}\left({\int }_0^{\pi }{\int }_0^{\pi }s(\varphi ,\theta ){\mathbf{a}}_0{\mathbf{a}}_0^{\mathrm{H}}(\varphi ,\theta )d\varphi d\theta \right){\mathbf{C}}^{*}\hfill \\ \hfill & \hfill ={\mathbf{C}}^{\mathrm{T}}{\mathbf{R}}_0{\mathbf{C}}^{*}\hfill \end{array}$$

(14)

which implies that each entry in the effective spatial correlation matrix $\mathbf{R}$ is determined by the collective effects of all the entries in the original spatial correlation matrix ${\mathbf{R}}_0$ weighted by the relevant entries in the coupling matrix $\mathbf{C}$. The theoretical analysis in (12)–(14) will be applied in Section 4 to study the performance of the proposed beamforming approach and the influence of MC on the effective spatial correlation of holographic RISs.

3. Eigenvalue Distributions Without MC

Denote the number of elements in a holographic RIS along the x and z directions as $N_x$ and $N_z$, respectively, i.e., the total number of elements $N=N_x{N}_z$. For the isotropic scattering environment, the spatial scattering function in (2) is [22]

$$s(\varphi ,\theta )=\frac{\sin \theta }{2\pi },\varphi \in \left[0,\pi \right],\theta \in \left[0,\pi \right].$$

(15)

Substituting (15) into (6) in Section 2.1 yields the spatial correlation matrix ${\mathbf{R}}_0$of the holographic RIS under isotropic scattering. As proven in [21,22], ${\mathbf{R}}_0$ can be characterized by a sinc function as follows:

$$\begin{array}{c}\hfill {\left[{\mathbf{R}}_0\right]}_{n_1,n_2}=\sin \mathrm{c}\left(\frac{2\left|\right|{\mathbf{d}}_{n_1}-{\mathbf{d}}_{n_2}{\left|\right|}_2}{\lambda }\right),n_1,n_2=1,\dots ,N\end{array}$$

(16)

where $\sin \mathrm{c}\left(x\right)\triangleq \frac{\sin \left(\pi x\right)}{\pi x}$ is the sinc function, ${\mathbf{d}}_{n_1}$ and ${\mathbf{d}}_{n_2}$ denote the coordinates of the $n_1$-th and $n_2$-th elements in the holographic RIS, respectively. The behavior of small-scale spatial correlation is depicted in Figure 2 for element spacing up to three times the wavelength $\lambda$. It is observed from (16) and Figure 2 that the spatial correlation is minimal only for some element spacing, instead of between any two elements, thus the classical independent and identically distributed (i.i.d.) Rayleigh fading model is not applicable in such a system [12,21,22].

Figure 3 portrays the eigenvalues of ${\mathbf{R}}_0/N$ in non-increasing order for various N, or equivalently various element spacing, with $L_x=L_z=12\lambda$. In addition, the dotted vertical line represents the asymptotic spatial DoF $\lceil \frac{\pi L_x{L}_z}{{\lambda }^2}\rceil$ for $\min (L_x,L_z)/\lambda \to \infty$[48,49]. A few key remarks can be drawn from Figure 3: First, while the popular i.i.d. Rayleigh fading channel has almost identical eigenvalues whose amount equals the number of antenna elements deployed, the correlated channel herein has uneven and fewer dominant eigenvalues (those larger than the ones highlighted by the dark circles in the inset of Figure 3) and smaller rank. Second, the accuracy of the spatial DoF $\lceil \frac{\pi L_x{L}_z}{{\lambda }^2}\rceil$ increases with the element density. More importantly, the number of dominant eigenvalues, highlighted by the dark circles in the inset of Figure 3, declines as the total number of elements N increases, which seems counter-intuitive. Therefore, in the following subsections, we will dig into the underlying causes of the aforementioned uncommon phenomenon of seemingly reduced spatial DoF with an increased number of elements in a holographic RIS.

3.1. Relationship between Eigenvalues and Power Spectrum

Note that the spatial correlation matrix ${\mathbf{R}}_0$ of the 2D holographic RIS in (16) can be formulated as

$${\mathbf{R}}_0=\left[\begin{array}{ccccc}{\mathbf{B}}_0& {\mathbf{B}}_1& {\mathbf{B}}_2& \cdots & {\mathbf{B}}_{N_z-1}\\ {\mathbf{B}}_{-1}& {\mathbf{B}}_0& {\mathbf{B}}_1& \cdots & {\mathbf{B}}_{N_z-2}\\ {\mathbf{B}}_{-2}& {\mathbf{B}}_{-1}& {\mathbf{B}}_0& \cdots & {\mathbf{B}}_{N_z-3}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\mathbf{B}}_{1-N_z}& {\mathbf{B}}_{2-N_z}& {\mathbf{B}}_{3-N_z}& \cdots & {\mathbf{B}}_0\end{array}\right]$$

(17)

where each block ${\mathbf{B}}_m\in {\mathbb{C}}^{N_x\times N_x}$, $\left|m\right|<N_z$, is a symmetric Toeplitz matrix by itself, and ${\mathbf{B}}_{-m}={\mathbf{B}}_m^T$, i.e., the entire matrix ${\mathbf{R}}_0$ is symmetric as well. Moreover, the matrix blocks ${\mathbf{B}}_m$’s along each diagonal of ${\mathbf{R}}_0$ are identical. Therefore, the spatial correlation matrix ${\mathbf{R}}_0$ has a symmetric BTTB structure under isotropic scattering, thus we resort to the relevant theory of the asymptotic distribution of eigenvalues of BTTB forms to investigate the properties of the eigenvalues of ${\mathbf{R}}_0/N$.

Denote the $(u,v)$th entry of ${\mathbf{B}}_m$ by $b_{l,m}=b_{u-v,m}$, $u,v=1,\dots ,N_x$, and define the function

$$g({\mathrm{\Delta }}_x,{\mathrm{\Delta }}_z)=\sin \mathrm{c}\left(\frac{2\sqrt{{\mathrm{\Delta }}_x^2+{\mathrm{\Delta }}_z^2}}{\lambda }\right)$$

(18)

where ${\mathrm{\Delta }}_x$ and ${\mathrm{\Delta }}_z$ are the spacing along the x- and z-axes, respectively, between a pair of spatial points of interest in the $xoz$ plane. $\left\{b_{l,m}\right\}$ can thus be regarded as a finite truncated bi-sequence generated from $g({\mathrm{\Delta }}_x,{\mathrm{\Delta }}_z)$, and more precisely,

$$b_{l,m}=\sin \mathrm{c}\left(\frac{2\sqrt{{\left(\frac{l{L}_x}{N_x}\right)}^2+{\left(\frac{m{L}_z}{N_z}\right)}^2}}{\lambda }\right).$$

(19)

Let $G({\omega }_x,{\omega }_z)$ be the 2D Fourier transform of $\left\{b_{l,m}\right\}$ given by

$$\begin{array}{cc}\hfill & \hfill G({\omega }_x,{\omega }_z)=\frac{1}{N}\sum _{l=-(N_x-1)}^{N_x-1}\sum _{m=-(N_z-1)}^{N_z-1}{b}_{l,m}{e}^{-j(l{\omega }_x+m{\omega }_z)},\hfill \\ \hfill & \hfill {\omega }_x=-\frac{(N_x-1)\pi }{N_x},-\frac{(N_x-3)\pi }{N_x},\dots ,\frac{(N_x-1)\pi }{N_x},\hfill \\ \hfill & \hfill {\omega }_z=-\frac{(N_z-1)\pi }{N_z},-\frac{(N_z-3)\pi }{N_z},\dots ,\frac{(N_z-1)\pi }{N_z}.\hfill \end{array}$$

(20)

Since the double-index sequence $\left\{b_{l,m}\right\}$ consists of spatial samples of the continuous sinc function $g({\mathrm{\Delta }}_x,{\mathrm{\Delta }}_z)$ in (18), $G({\omega }_x,{\omega }_z)$ in (20) can be looked upon as the power spectrum of the discretized and truncated spatial correlation function. According to the properties of BTTB matrices, the eigenvalues of ${\mathbf{R}}_0/N$ behave asymptotically the same as the spectral sampling points of $G({\omega }_x,{\omega }_z)$ as $N\to \infty$, if the double-index sequence formed by $b_{l,m}$ is absolutely summable [50,51], i.e.,

$$\underset{N_x,N_z\to \infty }{lim}\sum _{l=-(N_x-1)}^{N_x-1}\sum _{m=-(N_z-1)}^{N_z-1}\left|b_{l,m}\right|\le \mathrm{Constant}<\infty .$$

(21)

The condition in (21) is indeed satisfied as shown by (22),

$$\begin{array}{cc}\hfill & \hfill \underset{N_x,N_z\to \infty }{lim}\sum _{l=-(N_x-1)}^{N_x-1}\sum _{m=-(N_z-1)}^{N_z-1}\left|b_{l,m}\right|\hfill \\ \hfill =& \hfill \underset{N_x,N_z\to \infty }{lim}\sum _{l=-(N_x-1)}^{N_x-1}\sum _{m=-(N_z-1)}^{N_z-1}\left|\sin \mathrm{c}\left(\frac{2\sqrt{{\left(\frac{l{L}_x}{N_x}\right)}^2+{\left(\frac{m{L}_z}{N_z}\right)}^2}}{\lambda }\right)\right|\hfill \\ \hfill =& \hfill {\int }_{-1}^1{\int }_{-1}^1\left|\sin \mathrm{c}\left(\frac{2\sqrt{{\left(\varpi L_x\right)}^2+{\left(\varsigma L_z\right)}^2}}{\lambda }\right)\right|d\varpi d\varsigma \hfill \\ \hfill \le & \hfill {\int }_{-1}^1{\int }_{-1}^{1}1d\varpi d\varsigma =4<\infty \hfill \end{array}$$

(22)

hence the eigenvalues of ${\mathbf{R}}_0/N$ and the spectral sampling points of $G({\omega }_x,{\omega }_z)$ in (20) are asymptotically equally distributed. Consequently, insights on the eigenvalues of ${\mathbf{R}}_0/N$ can be drawn via the investigation of $G({\omega }_x,{\omega }_z)$. In what follows, we explore how the power spectrum $G({\omega }_x,{\omega }_z)$ changes with the element spacing of the holographic RIS, or equivalently, the spatial sampling frequency, in order to explain the unconventional observation of seemingly lower spatial DoF with growing numbers of elements in a holographic RIS as manifest in Figure 3.

3.2. Analysis on Eigenvalues via Power Spectrum

When the element spacing in a holographic RIS is ${\eta }_x\lambda$ and ${\eta }_z\lambda$ (${\eta }_x,{\eta }_z>0$) along the x and z directions, respectively, and $L_x={\beta }_x\lambda$, $L_z={\beta }_z\lambda$ (${\beta }_x,{\beta }_z>0$), we obtain $N_x={\beta }_x/{\eta }_x+1$, $N_z={\beta }_z/{\eta }_z+1$. The power spectrum $G({\omega }_x,{\omega }_z)$ in (20) can be recast as (23),

$$\begin{array}{cc}\hfill & \hfill G({\kappa }_x,{\kappa }_z)=\frac{1}{N}\sum _{l=-\frac{{\beta }_x}{{\eta }_x}}^{\frac{{\beta }_x}{{\eta }_x}}\sum _{m=-\frac{{\beta }_z}{{\eta }_z}}^{\frac{{\beta }_z}{{\eta }_z}}\sin \mathrm{c}\left(2\sqrt{{\left(\frac{l{\eta }_x(N_x-1)}{N_x}\right)}^2+{\left(\frac{m{\eta }_z(N_z-1)}{N_z}\right)}^2}\right)e^{-j\lambda (l{\eta }_x{\kappa }_x+m{\eta }_z{\kappa }_z)},\hfill \\ \hfill & \hfill {\kappa }_x=-\frac{(N_x-1)\kappa }{2N_x{\eta }_x},-\frac{(N_x-3)\kappa }{2N_x{\eta }_x},\dots ,\frac{(N_x-1)\kappa }{2N_x{\eta }_x},\hfill \\ \hfill & \hfill {\kappa }_z=-\frac{(N_z-1)\kappa }{2N_z{\eta }_z},-\frac{(N_z-3)\kappa }{2N_z{\eta }_z},\dots ,\frac{(N_z-1)\kappa }{2N_z{\eta }_z}\hfill \end{array}$$

(23)

in which $\kappa =2\pi /\lambda$ denotes the wavenumber; ${\kappa }_x$ and ${\kappa }_z$ represent the wavenumber along the x and z directions, respectively. The wavenumber along the positive y direction in Figure 1 is defined as

$$\varkappa ({\kappa }_x,{\kappa }_z)=\sqrt{{\kappa }^2-\left({\kappa }_x^2+{\kappa }_z^2\right)}$$

(24)

and the wave vector is

$$\mathsf{\kappa }={\kappa }_x\widehat{\mathbf{x}}+\varkappa ({\kappa }_x,{\kappa }_z)\widehat{\mathbf{y}}+{\kappa }_z\widehat{\mathbf{z}}$$

(25)

where $\widehat{\mathbf{x}}$, $\widehat{\mathbf{y}}$, and $\widehat{\mathbf{z}}$ denote the unit vector along the x-, y-, and z-axis, respectively. It is noteworthy that a real-valued $\varkappa ({\kappa }_x,{\kappa }_z)$ corresponds to a wave propagating along the y-direction, while an imaginary-valued $\varkappa ({\kappa }_x,{\kappa }_z)$ indicates an evanescent wave that usually exists around the surface of an object and decays exponentially in space [12,17].

The formulation in (23) allows us to study the power spectrum $G({\kappa }_x,{\kappa }_z)$ as a function of the wavenumbers ${\kappa }_x$ and ${\kappa }_z$, so as to gain more insight on the influence of the spatial sampling frequency on $G({\kappa }_x,{\kappa }_z)$, and, equivalently, on the eigenvalues of the normalized spatial correlation matrix ${\mathbf{R}}_0/N$, thanks to the time-frequency and space-wavenumber duality [48]. For a continuous and infinitely large holographic RIS aperture, $G({\kappa }_x,{\kappa }_z)$ in (23) approaches the following distribution [52]:

$$\begin{array}{c}\hfill G^{(\infty )}({\kappa }_x,{\kappa }_z)=\frac{\mathrm{\Pi }\left(\frac{\sqrt{{\kappa }_x^2+{\kappa }_z^2}}{2\kappa }\right)}{\frac{\kappa }{2\pi }\sqrt{{\kappa }^2-({\kappa }_x^2+{\kappa }_z^2)}}\end{array}$$

(26)

where $\mathrm{\Pi }(·)$ is the rectangle function. As implied by (26), $G^{(\infty )}({\kappa }_x,{\kappa }_z)$ assumes a bowl-like shape for a given $\kappa$, which increases monotonically with $\sqrt{{\kappa }_x^2+{\kappa }_z^2}$ inside the region

$$\mathcal{D}\left(\kappa \right)=\left\{({\kappa }_x,{\kappa }_z)\in {\mathbb{R}}^2:{\kappa }_x^2+{\kappa }_z^2\le {\kappa }^2\right\}$$

(27)

achieving the maximum value when $\sqrt{{\kappa }_x^2+{\kappa }_z^2}=\kappa$, and then transits abruptly to zero outside $\mathcal{D}\left(\kappa \right)$. In practical implementations, however, a holographic RIS often has a finite aperture size and is composed of discrete elements, which is equivalent to spatiallytruncating and sampling an originally infinite and continuous holographic RIS aperture. Truncating can be thought of as applying a windowing function to an originally infinitely-long signal in the spatial domain, which will cause the signal to appear outside $\mathcal{D}\left(\kappa \right)$ in the wavenumber domain after performing the Fourier transform, as displayed in Figure 4. Regarding spatial sampling, it is known from time-frequency-domain signal processing that a finer sampling granularity in the time domain yields a higher resolution in the (traditional) frequency domain. Analogously, for a holographic RIS, a smaller spatial sampling interval, which entails a smaller element spacing, gives rise to a higher resolution in the spatial frequency (i.e., wavenumber) domain. Accordingly, taking the x direction as an example, the highest absolute value of the resolvable wavenumber ${\kappa }_x$, namely $\frac{(N_x-1)\kappa }{2N_x{\eta }_x}$, is inversely proportional to the element spacing ${\eta }_x$ in (23). Specifically, if $N_x\gg 1$, for the most common element spacing of half a wavelength (i.e., ${\eta }_x=1/2$), the highest resolvable wavenumber is $\kappa$ as expected, and it increases to $\frac{\kappa }{2{\eta }_x}$ as ${\eta }_x$ decreases.

To examine the impact of the spatial sampling interval on the power spectrum $G({\kappa }_x,{\kappa }_z)$, we perform numerical simulations with a series of ${\eta }_x$ and ${\eta }_z$ values, and the corresponding results are shown in Figure 5. Figure 5a,c,e illustrate the power spectrum $G({\kappa }_x,{\kappa }_z)$ evaluated at ${\eta }_x={\eta }_z=1/2$, ${\eta }_x={\eta }_z=1/4$, and ${\eta }_x={\eta }_z=1/8$, respectively. For fair comparison, ${\kappa }_x$ and ${\kappa }_z$ range from $-4\kappa$ to $4\kappa$ in all of the three subfigures, while the actual wavenumber regime that the corresponding holographic RIS can resolve is within the white dashed frame in each of the three subfigures. Figure 5b,d,f are the zoom-in views of Figure 5a,c,e, respectively, where $|{\kappa }_x|\le \kappa$ and $|{\kappa }_z|\le \kappa$. Since $b_{l,m}$ in (19) can be treated as an aperiodic discrete-space signal, its Fourier transform (i.e., spectrum in the wavenumber domain) $G({\kappa }_x,{\kappa }_z)$ in (23) is periodic with periodicities of ${\kappa }_x/{\eta }_x$ and ${\kappa }_z/{\eta }_z$ along the x- and z-directions, respectively, as demonstrated in Figure 5a,c,e. Therefore, the spectrum in the wavenumber domain is actually superpositions of the original spectrum and its replicas, and the smaller the element spacing, the more serious the aliasing is. Due to the presence of spectrum leakage outside the region $\mathcal{D}\left(\kappa \right)$ in (27) as shown by Figure 4, the spectrum aliasing magnifies some of the spectrum values around the periphery of $\mathcal{D}\left(\kappa \right)$, and this magnification is more evident for smaller element spacing, as evidenced in Figure 5b,d,f. Recall that the eigenvalues of ${\mathbf{R}}_0/N$ are asymptotically equally distributed with the corresponding power spectrum $G({\kappa }_x,{\kappa }_z)$, hence the number of significant eigenvalues appears larger for smaller element spacing, which explains the observation highlighted by the dark circles in Figure 3. In other words, the number of eigenvalues lying in the far-field propagating wave region $\mathcal{D}\left(\kappa \right)$ in (27) actually does not change with the element spacing, the specious more dominant eigenvalues and hence higher spatial DoF for a smaller number of elements are contributed by evanescent waves outside the region $\mathcal{D}\left(\kappa \right)$.

Evanescent waves are usually localized in the reactive near field of an array. They contain the high-spatial-frequency components of an object, and normally do not contribute to the far-field channel capacity, but have a huge potential in near-field communication, sensing, and power transfer [53,54]. Moreover, diverse approaches have been put forward to convert evanescent waves to propagating waves that can be radiated to the far field, e.g., via simple metal strip gratings, metasurfaces, or randomly distributed metal wires [55,56,57]. Overall, holographic RISs, along with the associated propagating and evanescent waves, can find expansive potential applications in future communications, sensing, and related realms.

4. Eigenvalue Distributions with MC

In this section, we investigate eigenvalue distributions taking into account MC in holographic RISs.

4.1. Coupling Matrix

The impedance matrix $\mathbf{Z}$ of an antenna array can be expressed as

$$\mathbf{Z}=\left[\begin{array}{cccc}{z}_{\mathrm{A}}& z_{12}& \cdots & z_{1N}\\ z_{21}& z_{\mathrm{A}}& \cdots & z_{2N}\\ ⋮& ⋮& \ddots & ⋮\\ z_{N1}& z_{N2}& \cdots & z_{\mathrm{A}}\end{array}\right]$$

(28)

with $z_{\mathrm{A}}$ denoting the antenna impedance, and $z_{mn}$ the mutual impedance between the m-th and n-th elements. Given $\mathbf{Z}$, the (unnormalized) coupling matrix of the array at the Tx side can be calculated based on circuit theory as [30,34]

$${\mathbf{C}}_{\mathrm{T}}^{\prime }=\mathbf{Z}{\left(\mathbf{Z}+z_{\mathrm{S}}{\mathbf{I}}_N\right)}^{-1}$$

(29)

where $z_{\mathrm{S}}$ is the source impedance. If there is no MC between the Tx elements, then $\mathbf{Z}$ is diagonal with entries $z_{\mathrm{A}}$, hence ${\mathbf{C}}_{\mathrm{T}}^{\prime }$ is diagonal as well with entries ${\left[{\mathbf{C}}_{\mathrm{T}}^{\prime }\right]}_{nn}=z_{\mathrm{A}}/(z_{\mathrm{A}}+z_{\mathrm{S}})$. Consequently, we normalize the coupling matrix by $z_{\mathrm{A}}/(z_{\mathrm{A}}+z_{\mathrm{S}})$ to obtain

$${\mathbf{C}}_{\mathrm{T}}=\left(1+z_{\mathrm{S}}/z_{\mathrm{A}}\right)\mathbf{Z}{\left(\mathbf{Z}+z_{\mathrm{S}}{\mathbf{I}}_N\right)}^{-1}.$$

(30)

On the other hand, the (unnormalized) coupling matrix of the array at the Rx side is given by [30,34]

$${\mathbf{C}}_{\mathrm{R}}^{\prime }=z_{\mathrm{L}}{\mathbf{I}}_N{\left(\mathbf{Z}+z_{\mathrm{L}}{\mathbf{I}}_N\right)}^{-1}$$

(31)

where $z_{\mathrm{L}}$ is the load or termination impedance. Analogous to the Tx side, if no MC exists between the Rx elements, then ${\mathbf{C}}_{\mathrm{R}}^{\prime }$ is diagonal with entries ${\left[{\mathbf{C}}_{\mathrm{R}}^{\prime }\right]}_{nn}=z_{\mathrm{L}}/(z_{\mathrm{A}}+z_{\mathrm{L}})$, such that the normalized coupling matrix is

$${\mathbf{C}}_{\mathrm{R}}=\left(z_{\mathrm{A}}+z_{\mathrm{L}}\right){\left(\mathbf{Z}+z_{\mathrm{L}}{\mathbf{I}}_N\right)}^{-1}.$$

(32)

Given the coupling matrix, the effective spatial correlation matrix $\mathbf{R}$ for a holographic RIS can now be computed by applying the theoretical analysis in Section 2. The steps for deriving $\mathbf{R}$ are summarized below:

• Step 1: Calculate the conventional (i.e., MC-unaware) array response vector ${\mathbf{a}}_0$ in (5) based on the geometry and element topology of the RIS;
• Step 2: Obtain the spatial scattering function $s(\varphi ,\theta )$ in (3) based on theoretical analysis or measurements;
• Step 3: Calculate the spatial correlation matrix ${\mathbf{R}}_0$ that excludes MC according to (6);
• Step 4: Obtain the impedance matrix $\mathbf{Z}$ in (28) for the elements in a holographic RIS based on theoretical analysis or measurements;
• Step 5: Calculate the coupling matrix at the Tx and Rx, ${\mathbf{C}}_{\mathrm{T}}$ and ${\mathbf{C}}_{\mathrm{R}}$, according to (30) and (32), respectively;
• Step 6: Calculate the effective spatial correlation matrix $\mathbf{R}$ according to (14) using ${\mathbf{R}}_0$ from Step 3 and ${\mathbf{C}}_{\mathrm{T}}$ (${\mathbf{C}}_{\mathrm{R}}$) from Step 5 for the Tx (Rx).

4.2. Inter-Element Correlation/Coupling Strength Indicator

To facilitate the quantitative description of the physical MC or correlation strength represented by a coupling or correlation matrix $\mathbf{Q}\in {\mathbb{C}}^{N\times N}$, we propound a new metric called ICSI defined as

$$\begin{array}{cc}\hfill \mathrm{ICSI}& \hfill =\frac{1}{N}\sum _{n=1}^N\frac{\frac{1}{N-1}{\sum }_{m=1,m\ne n}^N\left|{\left[\mathbf{Q}\right]}_{n,m}\right|}{\left|{\left[\mathbf{Q}\right]}_{n,n}\right|}\hfill \\ \hfill & \hfill =\frac{1}{N(N-1)}\sum _{n=1}^N\sum _{m=1,m\ne n}^N\frac{\left|{\left[\mathbf{Q}\right]}_{n,m}\right|}{\left|{\left[\mathbf{Q}\right]}_{n,n}\right|}\hfill \end{array}$$

(33)

where $\frac{\left|{\left[\mathbf{Q}\right]}_{n,m}\right|}{\left|{\left[\mathbf{Q}\right]}_{n,n}\right|}$ represents the inter-element correlation/coupling strength between the n-th and m-th ($m\ne n$) elements normalized by the self-correlation/coupling magnitude of the n-th element to ease the comparison between different matrices, then the normalized inter-element correlation/coupling magnitude is averaged over all pairs of distinct array elements to arrive at the mean normalized inter-element correlation/coupling strength, i.e., the ICSI. The value of ICSI lies between 0 and 1, which indicate no MC/correlation and maximum MC/correlation between different array elements, and roughly correspond to the highest and lowest level of orthogonality among the columns/rows of the correlation or coupling matrix $\mathbf{Q}$, respectively. A large value of ICSI entails intense correlation/coupling between distinct array elements on average, which is likely to happen when the correlation/coupling between some pairs of elements is quite strong and/or a large amount of elements are mutually coupled or correlated, and this usually gives rise to unevenly distributed eigenvalues of the coupling/correlation matrix, as will be shown by the numerical results in the subsection below.

4.3. Numerical Simulations Including MC

Due to the complexity of the MC phenomenon, it is impossible to show generic formulas or curves for the impedance matrix or coupling matrix that apply to all types of elements or arrays. For ease of exposition, we adopt the analytical equations for MC between identical small dipoles as an example. Although the expressions are derived based on two dipoles, they may be applied to any pair of elements in a linear or planar array with or without a ground plane [58]. When the elements are small and resonant, such as dipoles, the first-order result of MC is to alter the impedance of each of the array elements [27]. Thus, the coupling can be described in terms of mutual impedance between elements (when higher order effects can be neglected), which is also a common practice in plenty of previous research work (e.g., [35,37,39]). The mutual impedance and MC models in this paper are intended to be illustrative; for arrays composed of dipoles with other layouts or a different type of elements, it is necessary to employ more suitable mutual impedance/MC models, or to measure these parameters in the actual array, which is deferred to future work. In this paper, we employ the impedance matrix of half-wave dipole elements in a parallel-in-echelon configuration as an example, as illustrated in Figure 6, where $z_{\mathrm{A}}\approx 73.1+j42.5\mathrm{\Omega }$ [27,59]. In the simulations, $L_x=4\lambda$ and the number of half-wavelength dipoles along the z-axis is set to eight. The spacing between the upper end of a dipole and the lower end of the adjacent one above it along the z-axis is set to $\lambda /50$, which is negligibly small compared with $\lambda$ so that $d_z\approx \lambda /2$ and $L_z\approx 4\lambda$.

As mentioned in Section 2.2, one of the distinctive features of a holographic RIS is that its array gain may be significantly larger than a conventional array; thus, let us first investigate the array gain of holographic RISs. As an instance, Figure 7 depicts the Tx array gain as a function of the azimuth angle $\varphi$ (see Figure 1) for both without and with MC cases. The zenith angle $\theta$ (see Figure 1) is set to ${90}^{\circ }$. $d_x$, $\mathbf{w}$, $\mathbf{C}$, and ${\mathbf{a}}_0$ denote the element spacing along the x-axis, beamforming vector, coupling matrix in (30), and MC-unaware array response vector in (5), respectively. The element spacing along the x-axis $d_x$ is set to $\lambda /2$ and $\lambda /8$, respectively, and four MC plus beamforming schemes are considered for each $d_x$. (Scheme 1): MC is included, and the proposed beamforming vector $\mathbf{w}$ in (12) is adopted. (Scheme 2): MC is included, and the ordinary conjugate beamforming ${\mathbf{a}}_0^{*}$ (followed by power normalization) is employed. (Scheme 3): MC is included, and the existing beamforming method in [37] is utilized which maximizes the directivity of the RIS and is equivalent to ${\mathbf{C}}^{-1}{\mathbf{a}}_0^{*}$ (followed by power normalization). (Scheme 4): MC is excluded, and the optimal conjugate beamforming ${\mathbf{a}}_0^{*}$ (followed by power normalization) is applied. The following remarks can be drawn from Figure 7: First, comparing the cases without and with MC, the array gain grows significantly in most cases when MC is included, even with the ordinary MC-unaware conjugate beamforming, except for some angles with the half-wavelength spacing. For example, the maximum achievable array gain using the proposed beamforming method is over $2.8$ times that without MC for $d_x=\lambda /8$, and the gap is likely to expand as the RIS becomes denser, which is promising for SNR enhancement, coverage extension, and energy-efficient transmission for green communications. Second, when excluding MC, the array gain ratio between $d_x=\lambda /8$ and $d_x=\lambda /2$ equals the ratio of the corresponding number of elements ($264/72\approx 3.7$ herein) as expected, while this ratio reaches $8.8$ when MC is included, indicating that the array gain incorporating MC increases more rapidly with the number of elements than the without MC scenario. Third, the proposed MC-aware beamforming approach outperforms its two counterparts including the one in [37]. The reason may lie in that the beamforming method in [37] is aimed to maximize the directivity of the RIS, which represents the radiation power for a certain pointing direction against that over all pointing directions, while the proposed beamforming vector maximizes the radiation power of an array against that of a single element, i.e., their optimization objects are slightly different. In general, the performance gaps among the three beamforming schemes increase with the element density and the proximity to the end-fire direction.

Now, we inspect the eigenvalue behaviors of the effective spatial correlation matrix in (14) at the Tx and Rx, respectively, for different element intervals with a fixed holographic RIS size. Figure 8 depicts the eigenvalue magnitude of the effective Tx spatial correlation matrix for both without and with MC cases. Ideally, the source impedance would be the conjugate of the element impedance, i.e., $z_{\mathrm{S}}=z_{\mathrm{A}}^{*}$, but this is difficult to realize in practice. We thereby consider both perfect and imperfect impedance match by setting $z_{\mathrm{L}}$ to $z_{\mathrm{A}}^{*}$ and a common value of $50\mathrm{\Omega }$ to examine the influence of different impedance matching.

Table 1. Inter-element correlation/coupling strength indicator (ICSI) in (33) for Tx spatial correlation matrix without and with MC under various element spacing.

ICSI in (33)	Tx Spatial Correlation Matrix without MC	Tx Spatial Correlation Matrix with MC, ${\mathit{z}}_{\mathbf{S}}={\mathit{z}}_{\mathbf{A}}^{*}$	Tx Spatial Correlation Matrix with MC, ${\mathit{z}}_{\mathbf{S}}=50\mathbf{\Omega }$	Tx Spatial Correlation Matrix with MC, ${\mathit{z}}_{\mathbf{S}}\mathbf{=}300\mathbf{\Omega }$
$d_x=\lambda /2$	0.0495	0.0927	0.0752	0.1500
$d_x=\lambda /4$	0.0646	0.0750	0.0665	0.0974
$d_x=\lambda /8$	0.0702	0.0705	0.0671	0.0851

lists the associated ICSI calculated using (33), which includes another very large source impedance of $300\mathrm{\Omega }$ for comparison purposes. The following main remarks can be drawn from Figure 8 and Table 1:

• When excluding MC, the most prominent inflection point of the eigenvalues in Figure 8 shifts to the left, i.e., the number of dominant eigenvalues decreases, as the element spacing shrinks, which is consistent with the observations for the larger holographic RIS aperture in Figure 3 and is well explained by the analysis in Section 3.
• In most cases, MC increases the effective spatial correlation at Tx, as indicated by the more rapidly-decaying eigenvalues compared with the no MC case in Figure 8 as well as the ICSI values in Table 1, and also elevates the largest eigenvalues for all element spacings studied. Furthermore, the correlation enhancement by MC diminishes as the array becomes denser, and MC may even have a decorrelation effect for sufficiently dense RISs. The reason lies in the fact that the product term $\mathbf{Z}{\left(\mathbf{Z}+z_{\mathrm{S}}{\mathbf{I}}_N\right)}^{-1}$ in (30) approaches an identity matrix when $z_{\mathrm{S}}\to 0$; meanwhile, for non-zero $z_{\mathrm{S}}$, it becomes a banded symmetric block-Toeplitz matrix that is sparse with non-zero entries confined to a diagonal band, and the off-diagonal entries are significantly smaller than the diagonal ones. In addition, the sparsity becomes more pronounced as the number of elements increases. Consequently, the behavior of ${\mathbf{C}}_{\mathrm{T}}$ in (30) gradually resembles that of a diagonal matrix as the element spacing dwindles, such that the effective spatial correlation becomes weaker for smaller element spacing when including MC.
• Different source impedance values exert a noticeable effect on the eigenvalue structure. Specifically, the effective spatial correlation is enhanced more substantially by perfect impedance match $z_{\mathrm{S}}=z_{\mathrm{A}}^{*}$ in contrast to $z_{\mathrm{S}}=50\mathrm{\Omega }$, as demonstrated by the corresponding eigenvalue trends in Figure 8 and ICSI values in Table 1. This can be explained by similar reasons stated above: the properties of ${\mathbf{C}}_{\mathrm{T}}$ in (30) are farther apart from those of a diagonal matrix with a larger $z_{\mathrm{S}}$ (i.e., $73.1-j42.5\mathrm{\Omega }$, as opposed to $50\mathrm{\Omega }$), thus higher correlation is induced. This is further verified by the ICSI for an even greater $z_{\mathrm{S}}$ of $300\mathrm{\Omega }$ in Table 1.
• In general, the eigenvalue magnitude increases with the density of the holographic RIS, which is consistent with the array gain enhancement shown by Figure 7.

The eigenvalue magnitude of the effective spatial correlation matrix at the Rx side is depicted in Figure 9 under the same simulation settings as in Figure 8, and the corresponding ICSI values are provided in

Table 2. ICSI in (33) for Rx spatial correlation matrix without and with MC under various element spacing.

ICSI in (33)	Rx Spatial Correlation Matrix without MC	Rx Spatial Correlation Matrix with MC, ${\mathit{z}}_{\mathbf{L}}={\mathit{z}}_{\mathbf{A}}^{*}$	Rx Spatial Correlation Matrix with MC, ${\mathit{z}}_{\mathbf{L}}=50\mathbf{\Omega }$	Rx Spatial Correlation Matrix with MC, ${\mathit{z}}_{\mathbf{L}}=300\mathbf{\Omega }$
$d_x=\lambda /2$	0.0495	0.0682	0.0787	0.0550
$d_x=\lambda /4$	0.0646	0.0867	0.1001	0.0698
$d_x=\lambda /8$	0.0702	0.1045	0.1213	0.0828

. Some interesting phenomena are observed and described below.

• MC at the Rx reduces the magnitude of eigenvalues in most cases (except for the first few largest eigenvalues). Compared with the Tx coupling matrix ${\mathbf{C}}_{\mathrm{T}}$ in (30), the Rx coupling matrix ${\mathbf{C}}_{\mathrm{R}}$ in (32) is mainly different by a term $z_{\mathrm{A}}{\mathbf{Z}}^{-1}$, which causes the reduction of the eigenvalue magnitude after multiplying with the original spatial correlation matrix.
• As seen in Table 2, MC increases the effective spatial correlation at the Rx, which is similar to the Tx side. However, contrary to the observation for the Tx side, MC with a larger load impedance tends to have less correlation enhancement effect at Rx in general, as shown by Figure 9 and Table 2, since a larger load impedance renders the term $z_{\mathrm{L}}{\mathbf{I}}_N$ more dominant in ${\left(\mathbf{Z}+z_{\mathrm{L}}{\mathbf{I}}_N\right)}^{-1}$ in (32) so that the coupling matrix is more like a diagonal matrix with smaller off-diagonal entries and hence weaker correlation.
• Rx MC increases the effective spatial DoF for element spacings less than half a wavelength, as implied by the most prominent inflection point of the eigenvalues which shifts to the right when considering MC. This is beneficial for spatial diversity and multistreaming.

To analyze and compare the MC effects of different element types, we also look at the element with an isotropic radiation pattern which is commonly assumed in the existing work. The real part of the impedance matrix of isotropic elements is [47]

$$\begin{array}{c}\hfill {\left[\mathfrak{R}\left\{{\mathbf{Z}}_{\mathrm{iso}}\right\}\right]}_{n_1,n_2}=r_{\mathrm{iso}}\sin \mathrm{c}\left(\frac{2\left|\right|{\mathbf{d}}_{n_1}-{\mathbf{d}}_{n_2}{\left|\right|}_2}{\lambda }\right),n_1,n_2=1,\dots ,N\end{array}$$

(34)

where $r_{\mathrm{iso}}$ denotes the radiation resistance of an isotropic element, and the definitions of $\sin \mathrm{c}\left(x\right)$, ${\mathbf{d}}_{n_1}$, and ${\mathbf{d}}_{n_2}$ are aligned with those in (16). By applying appropriate matching techniques, the imaginary part of the impedance matrix can be removed while the real part remains [47], thus we only need to consider the real part. The Rx coupling matrix for isotropic elements is obtained by plugging (34) into (32) and setting $z_{\mathrm{L}}=r_{\mathrm{iso}}$. Figure 10 displays the eigenvalue magnitude versus eigenvalue index of the effective Rx spatial correlation matrix for half-wavelength dipoles and isotropic elements, which shows that the effective spatial correlation for isotropic elements is obviously lower than that for half-wavelength dipoles, as manifest from the more evenly distributed eigenvalue magnitude for isotropic elements. This is because the MC is zero whenever the spacing between two isotropic elements is multiples of half a wavelength, as indicated by (34), while the MC is all non-zero for the dipoles herein, hence the overall MC is lower for isotropic elements.

5. Conclusions

In this paper, we have investigated the array response and small-scale spatial correlation for holographic RISs excluding and including MC. In-depth analysis is conducted on the asymptotic eigenvalue distribution and spatial DoF for the spatial correlation matrix under isotropic scattering, by linking the eigenvalues to the power spectrum of the spatial correlation function based on the BTTB matrix theory. It is demonstrated that the specious more dominant eigenvalues for fewer elements in holographic RISs are due to the spectrum aliasing in the wavenumber domain which enhances the near-field evanescent waves, thus the far-field spatial DoF actually does not increase. Furthermore, an MC-aware beamforming scheme is proposed which is aimed to maximize the array gain, and is shown to outperform existing methods. In addition, it is found that MC exerts discrepant effects on Tx and Rx modes: For the Tx, MC increases the eigenvalue magnitude and effective spatial correlation in most cases, while Rx MC often reduces the eigenvalue magnitude while increasing the spatial DoF, especially for dense RISs. The array gain and channel eigenvalue enhancement by Tx MC is beneficial to SNR elevation, coverage extension, and energy saving, and the growing spatial DoF caused by Rx MC indicates that efficient spatial multiplexing is possible even with compact arrays.