1 Introduction↩︎

Standing as a groundbreaking technique in modern wireless communication systems, riss represent a revolutionary paradigm for the smart em wave control [1]–[4]. These surfaces, composed of an array of programmable reflective or refractive unit cells, empower real-time adjustments to the radio environment. Leveraging advanced metamaterials and signal processing techniques, RIS enables precise modulation of the phase, amplitude, and polarization of em waves. Such adaptability promises to significantly benefit wireless applications including communication, sensing, localization, and mapping [5]–[9]. With the advent of densely integrated RISs such as holographic RIS [10], the phenomenon of mutual coupling emerges as a critical consideration with profound implications for the performance of RIS-assisted communications [1]. Mutual coupling refers to the interaction between adjacent RIS unit cells, wherein the em fields generated by one unit cell influence the behavior of its neighboring unit cells. This inter-element coupling effect can significantly impact the overall functionality and effectiveness of a RIS deployment, especially in scenarios involving high amplitude gains or varying operational frequencies. Therefore, understanding and mitigating the impact of mutual coupling is essential for optimizing the performance of RIS systems.

1.1 Prior Work and Motivation↩︎

RIS-aided communication has been extensively studied in recent years owing to its substantial potential in extending communication coverage [11], enhancing signal reception [12], and improving energy efficiency [13]. To date, numerous novel signal processing techniques have been developed to support RIS-aided wireless applications [2], [7], [14]. Besides those efforts in signal processing techniques, significant advancements have been made in ris hardware. On the one hand, some designs focus on novel RIS concepts and implementations. For instance, a wideband RIS that fully covers the n257 and n258 5G mmWave bands was reported in [15]. The authors of [16] designed an active RIS with consistent reflection polarization using reflection amplifiers. Furthermore, the authors of [17] proposed a RIS design featuring absorptive filtering capabilities to suppress interference. On the other hand, some RIS developments incorporate advanced fabrication and control techniques. For example, the authors of [18] designed an RIS with robust remote infrared control. Additionally, the authors of [19] and [20] introduced the concept of fully printed RIS, including both metallic patterns and switches, which substantially reduces fabrication costs, especially for large-scale RISs.

Regarding the RIS mutual coupling, numerous recent investigations have led to significant advancements in modeling and performance optimization. A physics-based end-to-end model, named PhysFad, has been proposed in [21] based on a first-principles coupled-dipole formalism. This model characterizes RIS-parametrized wireless channels with adjustable fading. Additionally, a linear series expansion-based analysis of the channel’s non-linearity mechanism has been presented in [22], and an experimental validation using a 1-bit RIS prototype operating in the cmWave band has been reported in [23]. Stemming from fundamental em principles, the PhysFad model inherently accounts for proximity-induced mutual coupling, reverberation-induced multi-reflection, and environmental fading, providing essential insights into accurate ris channel modeling, particularly in rich-scattering environments. With its EM-compliant characteristics, this channel model serves as a reliable benchmark for generating physically precise wireless channels and verifying advanced RIS configuration designs [24].

Another suitable theory to model mutual coupling in RIS is the microwave network theory. Notably, existing microwave network-based communication modeling can be categorized into two types: impedance matrix-based models and scattering matrix-based models. In the impedance matrix-based approach, the first effort can be found in [25], where an em-compliant and mutual coupling-aware communication model is proposed by adopting the mutual impedance analysis among RIS unit cells. This model has proven effective in guiding RIS configuration, facilitating applications such as end-to-end received power maximization in siso systems [26] and sum-rate optimization in multi-user interference mimo channels [27]. On the other hand, the scattering matrix-based model was initially explored by the authors of [28] through a scattering parameter network analysis. Subsequently, the authors of [29] derived a scattering matrix-based mutual coupling-aware communication model. The scattering parameter is regarded as a model that is more directly related to the radiation pattern and easier to measure [30]. Importantly, the derivation in [29] and [31] establish the equivalence between the impedance matrix-based and scattering matrix-based communication models. Recently, microwave network-based models, either in the impedance matrix form or the scattering matrix form, have been successfully applied, for example, in em modeling for arbitrarily shaped ris[32] and in ris-parametrized wireless networks-on-chip [33].

This paper aims to provide a practical method for evaluating mutual coupling in RIS, thereby facilitating channel modeling and downstream applications. Due to the passive reflective nature of RIS, the accurate coupling parameters between its unit cells cannot be directly acquired through simulation or measurement. Notice that in far-field and sparse-scattering environments, the wireless channel can be accurately characterized by the conventional cascaded channel model, as evidenced by various studies [34]–[36]. In pursuit of both accurate mutual coupling description and model tractability, this paper adopts microwave network theory-based modeling to propose a novel scattering parameter estimation method. This method features a limited number of unknowns, independent of the number of ris elements, as long as it constitutes a uniform array. Tailored for sparse-scattering propagation in mmWave/THz frequency, the adopted microwave network theory enables the separation of RIS mutual coupling and environmental scattering in a cascaded manner, with negligible loss of accuracy, by disregarding mutual coupling with the transmitter and receiver under the far-field assumption [28], [31]. Leveraging this model, the proposed approach can determine these limited scattering parameters based on the sample observations from a single full-wave simulation. The accuracy of the estimated results is validated through both simulation and experiment, encompassing assessment of the cross-angle radiation pattern and end-to-end S-parameter predictions.

1.2 Main Contributions↩︎

The main contributions of this paper are summarized as follows.

Model Training: Utilizing the microwave network-based channel model,a realistic evaluation method of the coupling parameters among RIS unit cells is proposed. By only considering the coupling between adjacent RIS unit cells and neglecting very weak coupling effects beyond the distance of two unit cells, we parameterize the scattering matrix of the RIS into a finite set of scattering parameters. Based on this parameterized approximation, along with the far-field cascaded wireless channel model, we design a training approach to estimate these scattering parameters based on a single 3D full-wave em simulation of the radiation pattern. The simulated RIS radiation pattern accurately captures the 3D em behavior including mutual coupling, thus we can estimate these scattering parameters by optimizing their fitness to the simulated pattern. We show that the estimation can be formulated into a convex optimization problem via the Neumann series expansion, accompanied by a rigorous proof of convexity. Simulation Validation: The mutual coupling-aware communication model with the trained scattering parameters is validated through the full-wave simulation. First, we perform the full-wave simulations of 3 different far-field radiation patterns of the RIS by varying the position of the Tx antenna. In the meantime, the theoretical predictions of these radiation patterns are obtained using the mutual coupling-aware model (with the trained scattering parameters) and the conventional mutual coupling-unaware model, respectively. Benchmarked against the simulated radiation patterns, the trained model exhibits superior accuracy compared to the conventional model without accounting for mutual coupling, validating the precision of the trained model.

Experimental Validation: To further confirm the trained model, experimental validation is also carried out in this work. The experimental measurements are conducted on the 1-bit mmWave RIS prototype reported in [15]. We use two horn antennas to respectively transmit and receive em signals, and measure the S-parameter between the two horn antennas through a vector network analyzer. Afterward, this S-parameter is predicted by the trained model and the conventional model, respectively. The comparisons over different setups demonstrate that the trained model can predict the S-parameter closer to the real measurement, revealing both the accuracy of the adopted mutual coupling-aware communication model and the effectiveness of the proposed scattering parameters training approach. Moreover, our experimental results demonstrate that while these S-parameters are derived from far-field modeling, they can also accurately predict near-field radiation.

The paper is organized as follows. Section II introduces the fundamentals of RIS hardware. Section III recaps the RIS-aided communication models, including both the conventional mutual coupling-unaware model and the recently proposed mutual coupling-aware model based on scattering matrices. Leveraging a 3D full-wave simulation of the RIS radiation pattern, a training approach for the scattering parameters among RIS unit cells is proposed in Section IV. Based on the trained model, the simulation validation and the experimental validation are conducted in Section V and Section VI, respectively. Finally, the conclusions of the paper are drawn in Section VII.

2 Fundamentals of RIS Hardware↩︎

2.1 Functional Principle of RIS↩︎

A ris is a planar surface that consists of an array of scattering unit cells, each of which can impose an independent phase shift, and possibly an amplitude gain, on the incident em waves. Fig. 1 depicts a conceptual structure of the RIS. As illustrated, each RIS unit cell is equivalently comprised of a reflecting element and a functional circuit with tunable load \(Z_{\text{RIS},m}\) [28], [37], [38], where \(m\) denotes the index of the unit cell. By carefully adjusting the tunable loads of all the unit cells, the desired phase shifts (and amplitude gains) can be achieved thus the reradiated em waves can be beamformed toward specified directions. As shown in Fig. 1, different ris unit cells, e.g., \(Z_{\text{RIS},m}\) and \(Z_{\text{RIS},n}\), can be further connected through a reconfigurable impedance \(Z_{\text{RIS},m,n}\), which is referred to as the beyond diagonal ris or group/fully-connected ris as recently proposed in [29], [39]. For clarification, this work uses the classical ris with independent unit cells, which is also called the single-connected ris.

Each unit cell of the ris can be treated as a 1-port network. The voltage waves impinging on the port (denoted as \(V_\mathrm{in}\)) and reflected from the port (denoted as \(V_\mathrm{re}\)) are related through the reflection coefficient \(\Theta\) as \[V_\mathrm{re} = \Theta V_\mathrm{in}.\] That is, a phase shift \(\angle\Theta\) and an amplitude gain \(|\Theta|\) are imposed on the incident em wave. The reflection coefficient \(\Theta\) is determined by the load impedance \(Z_\mathrm{RIS}\) of the ris unit cell through the relationship [28], [40] \[\label{eq95ThetaSin} \Theta = \frac{Z_\mathrm{RIS}-Z_0}{Z_\mathrm{RIS}+Z_0},\tag{1}\] where \(Z_0\) is the reference impedance (or the characteristic impedance) and usually \(Z_0=\unit[50]{\Omega}\). Depending on the real part of the load impedance \(Z_\mathrm{RIS}\), various types of RIS can be distinguished as follows.

• \(\mathrm{Re}(Z_\mathrm{RIS})>0\): This yields an amplitude gain \(|\Theta|<1\) and is the case of the most real hardware of the quasi-passive RIS [15].

• \(\mathrm{Re}(Z_\mathrm{RIS})=0\): This yields a unit-modulus constraint \(|\Theta|=1\) which is widely adopted in the existing theoretical work for simplicity. The unit-modulus reflection coefficient can be realized through a purely reactive RIS load [2], [8], [13], [28], [41].

• \(\mathrm{Re}(Z_\mathrm{RIS})<0\): This yields an amplitude gain \(|\Theta|>1\) which is known as the active ris [16]. Several technologies can realize a negative resistance. For example, by imposing proper bias voltage, the tunnel diode can work in the negative differential resistance regions thus leading to an active power gain [42].

2.2 A Fabricated Quasi-Passive RIS Prototype↩︎

Based on the functional principle in Section 2.1, a practical wideband 1-bit quantized quasi-passive RIS operating at the 5G mmWave band has been designed and fabricated for the verification of the mutual coupling effect. This RIS prototype utilizes a single-connected network architecture [39]. The RIS unit cell design, array synthesis, simulation, and fabricated prototype are demonstrated as follows.

2.2.1 RIS Unit Cell Design↩︎

The detailed design methodology for wideband RIS unit element has been investigated in our previous work [15]. The reflection phase and reflection magnitude of the designed unit cell (with \(^\circ\) incidence ) when the PIN switch at ON/OFF states are shown in Fig 2 Such a reflection performance is realized by employing patches with circular cutouts [43], [44], two long vias, one short via, and PIN diodes to excite the unit cell with two adjacent high-order harmonic resonances. Basically, the 1-bit RIS unit cell design has two states when the PIN diodes are switched ON and OFF. Ideally, these two states should have full reflection (magnitude equal to 1 for quasi-passive RIS) and \(\unit[180]{^\circ}\) phase difference for 1-bit phase quantization. However, considering the imperfections of the practical realization, a phase difference of \(\unit[180]{^\circ}\pm\unit[20]{^\circ}\) is acceptable, which maintains a decent RIS performance. We can observe a reflection amplitude larger than \(\unit[-2.8]{dB}\) and a phase difference within \(\unit[180]{^\circ}\pm\unit[20]{^\circ}\) for the bandwidth from \(\unit[22.7]{}\) to \(\unit[30.5]{GHz}\). For example, at the frequency of \(\unit[28]{GHz}\) (as highlighted in Fig. 2), each RIS unit cell has the phase and amplitude responses as follows. \[\angle\Theta = \begin{cases} \unit[-193.3]{^\circ} & \text{if } \text{PIN state} = \text{ON},\\ \unit[-388.1]{^\circ} & \text{if } \text{PIN state} = \text{OFF}, \end{cases}\] \[|\Theta| = \begin{cases} \unit[-1.57]{dB} & \text{if } \text{PIN state} = \text{ON},\\ \unit[-0.47]{dB} & \text{if } \text{PIN state} = \text{OFF}. \end{cases}\] Here, the reduced reflection magnitude in the ON state compared to the OFF state is attributed to the PIN diode losses, whereas there are no such losses in the OFF state.

The designed RIS has a small unit cell size of \(0.35\lambda\times 0.35\lambda\) at the center operation frequency of 27.5 GHz. Based on the existing results in the literature, the mutual coupling effect among RIS unit cells becomes severe when the inter-element distance is shorter than half-wavelength [26]. Therefore, the designed RIS prototype should have an observable strong mutual coupling between adjacent RIS cells, suggesting its suitability for the experimental validation of mutual coupling.

2.2.2 RIS Array Simulation↩︎

Before the fabrication and measurement, the radiation performance of the designed RIS can be evaluated through the full-wave simulation based on the finite element method (FEM), which captures precise em characteristics including mutual coupling. For this design, the RIS is composed of 400 unit elements in a 20 × 20 configuration, as demonstrated in Fig. 3. The simulation of the RIS with 400 elements and an incident horn antenna has almost reached the computation limit of the workstation with GB RAM. Fig. 4 illustrates a sample of the simulated 3D radiation pattern along with the RIS and in a polar plot. It should be noted that different radiation patterns can be achieved by varying the signal incident angle, RIS beamforming, and the distance between RIS and the incident horn antenna. However, the mutual coupling parameters, as will be shown in Section 3.2, are independent of these system setups and only depend on the physical layout of RIS reflecting elements. Therefore, the simulated radiation patterns under different system setups can be utilized to train and test the model parameters of the mutual coupling, since it is unable to measure them directly.

2.2.3 Fabricated Prototype↩︎

Typically, a completely functional RIS design is comprised of a periodic structure with magnitude/phase modulation and a digital circuit for the array pattern generation and control. In this work, the designed control circuits and the packaged RIS prototype are shown in Fig. 5. Specifically, the designed circuit board links with a microcontroller unit (MCU) WROOM-32 module that can compute the array patterns, as illustrated in Fig. 5-(a). The circuit is designed with 8-bit phase registers (74HC595D) and PMOS to bias the PIN switches, as shown in Fig. 5-(b). When the output of the phase register is at the high voltage level, the PIN switch is turned on, thus the RIS unit cell, as demonstrated in Fig. 5-(c), is in the ON state. Otherwise, the corresponding unit cell is in the OFF state. The phase distribution map of the RIS can be visualized through an LED array in the back of the control circuit. Fig. 5-(d) demonstrates the overall fabricated RIS prototype, including the RIS array, control circuit board, connectors, and flexible printed circuit (FPC), which are packaged together through 3D-printed mechanical supports using a Raise3D Pro2 printer with polylactic acid filaments. Note that the function of FPC here is to connect the RIS array with the whole designed circuits inside the package.

3 RIS-Aided Communication Models↩︎

As a foundation of the subsequent analyses, this section reviews the RIS-aided communication models, including both the conventional model that does not account for RIS mutual coupling [45] and the recently proposed mutual coupling-aware communication model [29]. Consider a RIS-aided mimo communication with an \(M_\mathrm{T}\)-antenna transmitter, an \(M_\mathrm{R}\)-antenna receiver, and a RIS with \(N\) unit cells. Let \({\boldsymbol{v}}_\mathrm{R}\in\mathbb{C}^{M_\mathrm{R}\times 1}\) and \({\boldsymbol{v}}_\mathrm{T}\in\mathbb{C}^{M_\mathrm{T}\times 1}\) denote the total voltages at the receiver and the transmitter, respectively. The wireless channel is represented by the channel matrix \({\boldsymbol{H}}\in\mathbb{C}^{M_\mathrm{R}\times M_\mathrm{T}}\), which builds the connection between \({\boldsymbol{v}}_\mathrm{R}\) and \({\boldsymbol{v}}_\mathrm{T}\) as \[{\boldsymbol{v}}_\mathrm{R} = {\boldsymbol{H}}{\boldsymbol{v}}_\mathrm{T}.\]

3.1 Conventional Channel Model without Mutual Coupling↩︎

A commonly adopted RIS-aided communication channel \({\boldsymbol{H}}\) can be expressed as [2], [8], [13], [38], [45], [46] \[\label{eq95Hconven} {\boldsymbol{H}}= {\boldsymbol{H}}_\mathrm{RT} + {\boldsymbol{H}}_\mathrm{RI}\boldsymbol{\Theta}{\boldsymbol{H}}_\mathrm{IT},\tag{2}\] where \({\boldsymbol{H}}_\mathrm{RT}\in\mathbb{C}^{M_\mathrm{R}\times M_\mathrm{T}}\) denotes the Tx-Rx channel, \({\boldsymbol{H}}_\mathrm{IT}\in\mathbb{C}^{N\times M_\mathrm{T}}\) denotes the Tx-RIS channel, \({\boldsymbol{H}}_\mathrm{RI}\in\mathbb{C}^{M_\mathrm{R}\times N}\) denotes the RIS-Rx channel, and \(\boldsymbol{\Theta}\in\mathbb{C}^{N\times N}\) is the RIS reflection matrix. Here, \(\boldsymbol{\Theta}\) is a diagonal matrix whose diagonal entries are given by 1 . The expressions of these channel matrices \({\boldsymbol{H}}_\mathrm{RT}\), \({\boldsymbol{H}}_\mathrm{RI}\), and \({\boldsymbol{H}}_\mathrm{IT}\) will be specified in Section 6.

Although channel model 2 is widely utilized, it disregards the interactions between adjacent RIS cells, i.e., the mutual coupling effect. In general, the mutual coupling in an antenna array/RIS can be reasonably ignored if the inter-element distance is large enough (e.g., larger than the half wavelength [26], [47]). Nonetheless, the significance of mutual coupling becomes increasingly pronounced with the advent of techniques like holographic mimo with extremely compact array integration  [48], [49], yielding an inter-element distance much shorter than the half wavelength. Therefore, the mutual coupling-aware communication models are required.

3.2 Mutual Coupling-Aware Channel Model↩︎

To account for the mutual coupling effect in RIS-aided communication, several novel communication models have been recently proposed. These models can be categorized into two types: the impedance matrix (Z-parameters)-based model [25] and the scattering matrix (S-parameters)-based model [28]. It should be noted that based on microwave network theory, these two types of models are essentially equivalent and can be seamlessly converted to each other. Such a relationship has been thoroughly analyzed in [29] and [31].

In this work, we adopt the scattering matrix-based communication model, which was first reported in [29]² as \[\label{eq95simpH2} {\boldsymbol{H}}= {\boldsymbol{H}}_\mathrm{RT}+{\boldsymbol{H}}_\mathrm{RI}\big(\boldsymbol{\Theta}^{-1}-{\boldsymbol{S}}\big)^{-1}{\boldsymbol{H}}_\mathrm{IT},\tag{3}\] where \({\boldsymbol{S}}\in\mathbb{C}^{N\times N}\) denotes the scattering matrix of the RIS reflecting element network and thus characterizes the mutual coupling. Specifically, the entry \(S_{ij}\) in \({\boldsymbol{S}}\) is the S-parameter between the \(i\)-th element and the \(j\)-th element on the RIS, which indicates the voltage wave measured at the \(i\)-th element when a unit voltage wave is driven at the \(j\)-th element [40]. By reciprocity, we have \(S_{ij}=S_{ji}\) thus \({\boldsymbol{S}}\) is a symmetric matrix. Additionally, the diagonal entry \(S_{ii}\) denotes the self-scattering of the \(i\)-th reflecting element. This model 3 can also be found in, e.g., [31] or [50].

4 Scattering Model Training↩︎

Now, we have two communication models 2 and 3 . However, as can be observed in 3 , to utilize this mutual coupling-aware model, we still need to determine the scattering matrix \({\boldsymbol{S}}\). Since our RIS prototype only reflects em waves passively, this scattering matrix \({\boldsymbol{S}}\) cannot be measured or simulated directly. Hence, in this paper, we first propose a model training method that can estimate matrix \({\boldsymbol{S}}\) based on a single 3D full-wave simulation of the RIS radiation pattern. Then, based on the estimated scattering matrix \(\hat{{\boldsymbol{S}}}\), we can evaluate and compare the accuracy of the two communication models. In this section, we first introduce the theoretical calculation of the RIS radiation pattern based on the two communication models. Next, we propose a model training method leveraging such theoretical models and a 3D full-wave simulation of the RIS radiation pattern.

4.1 Theoretical Radiation Pattern Calculation↩︎

Since our measurement facility with horn Tx and Rx antennas form a siso setup, we now present the calculation of the RIS radiation pattern based on the siso version of 2 and 3 . Aiming to evaluate the RIS mutual coupling, we focus on the Tx-RIS-Rx channel only and block the Tx-Rx path that is unrelated to RIS. The radiation pattern refers to a function of the radiation properties of an array as a function of the directional coordinates [51].³ By fixing the position of the Tx horn antenna, we can calculate the radiation pattern of the Tx\(+\)RIS system. We denote the directional coordinates by the aod \(\boldsymbol{\theta}=[\theta_\mathrm{az},\theta_\mathrm{el}]^\mathsf{T}\) departing from the RIS, which consists of an azimuth angle \(\theta_\mathrm{az}\) and an elevation angle \(\theta_\mathrm{el}\). In this paper, the elevation angle is defined as the angle between the departure direction and the normal direction of the array plane, as utilized in [51].

4.1.1 Radiation Pattern Calculation Based on 2↩︎

Based on the conventional communication model 2 , different RIS unit cells reflect signals independently without interaction. Hence, the (amplitude) radiation pattern in linear scale can be calculated by [34] \[\begin{align} \label{eq95E} &E(\boldsymbol{\theta})\\ =&\Big|\!\sum_{n=1}^N \Theta_{n,n}\frac{\cos^{q_\mathrm{e}}(\theta_\mathrm{el})\cos^{q_\mathrm{f}}(\theta_{\mathrm{T},n})}{\|{\boldsymbol{p}}_n-{\boldsymbol{p}}_\mathrm{T}\|_2}e^{jk_0({\boldsymbol{p}}_n^\mathsf{T}{\boldsymbol{u}}(\boldsymbol{\theta})-\|{\boldsymbol{p}}_n-{\boldsymbol{p}}_\mathrm{T}\|_2)}\Big|,\notag \end{align}\tag{4}\] where \(k_0={2\pi}/{\lambda}\) is the wavenumber, \({\boldsymbol{p}}_n\in\mathbb{R}^3\) denotes the position of the \(n\)-th RIS element, \({\boldsymbol{p}}_\mathrm{T}\in\mathbb{R}^3\) denotes the position of the Tx horn, and \({\boldsymbol{u}}(\boldsymbol{\theta})\) denotes the unit directional vector of the aod \(\boldsymbol{\theta}\). All these positions and vectors are expressed in the body coordinate system of the RIS. Here, \(\lambda=c/f\) denotes the em wavelength, where \(f\) is the signal frequency and \(c\) is the speed of light. Besides, \(\theta_{\mathrm{T},n}\) denotes the elevation angle of the aod from Tx to the \(n\)-th RIS unit cell expressed in the coordinate system of the Tx horn, and \(q_\mathrm{e}\) and \(q_\mathrm{f}\) are the feed power patterns which determine the radiation directivities of the RIS unit cell and Tx horn antenna, respectively. In addition, \(\Theta_{n,n}\) denotes the reflection coefficient of the \(n\)-th RIS unit cell, i.e., the \(n\)-th diagonal entry of \(\boldsymbol{\Theta}\).

We can reform 4 into a more compact form as \[\label{eq95Econv} E(\boldsymbol{\theta}) = \big|{\boldsymbol{h}}_{\mathrm{out}}^\mathsf{T}(\boldsymbol{\theta})\boldsymbol{\Theta}{\boldsymbol{h}}_\mathrm{in}\big|,\tag{5}\] where the entries in \({\boldsymbol{h}}_\mathrm{in}\in\mathbb{C}^{N\times 1}\) and \({\boldsymbol{h}}_{\mathrm{out}}(\boldsymbol{\theta})\in\mathbb{C}^{N\times 1}\) are given by \[\begin{align} \tag{6} h_{\mathrm{in},n} &= \frac{\cos^{q_\mathrm{f}}(\theta_{\mathrm{T},n})}{\|{\boldsymbol{p}}_n-{\boldsymbol{p}}_\mathrm{T}\|_2}e^{-jk_0\|{\boldsymbol{p}}_n-{\boldsymbol{p}}_\mathrm{T}\|_2},\\ \tag{7} h_{\mathrm{out},n}(\boldsymbol{\theta}) &= \cos^{q_\mathrm{e}}(\theta_\mathrm{el})e^{jk_0{\boldsymbol{p}}_n^\mathsf{T}{\boldsymbol{u}}(\boldsymbol{\theta})}. \end{align}\]

4.1.2 Radiation Pattern Calculation Based on 3↩︎

By replacing the RIS response \(\boldsymbol{\Theta}\) with \((\boldsymbol{\Theta}^{-1}-{\boldsymbol{S}})^{-1}\) according to 3 , the mutual coupling-aware RIS radiation pattern can be computed as \[\label{eq95rMC} {E}(\boldsymbol{\theta}) = \big|{\boldsymbol{h}}_{\mathrm{out}}^\mathsf{T}(\boldsymbol{\theta})\big(\boldsymbol{\Theta}^{-1}-{\boldsymbol{S}}\big)^{-1}{\boldsymbol{h}}_\mathrm{in}\big|.\tag{8}\]

For clarification, we denote the radiation pattern calculated based on 8 as \(E(\boldsymbol{\theta};{\boldsymbol{S}})\), showing it is a function of scattering matrix \({\boldsymbol{S}}\). In contrast, we denote the radiation pattern calculated through 5 as \(E(\boldsymbol{\theta};\mathbf{0})\) for distinction. Besides, radiation pattern of an array is usually normalized with respect to their maximum value, yielding normalized radiation patterns. To do so, we define the normalized version of 5 and 8 as \[\begin{align} E_\mathrm{n}(\boldsymbol{\theta};\mathbf{0}) &= E(\boldsymbol{\theta};\mathbf{0})/\big(\max_{\boldsymbol{\theta}} E(\boldsymbol{\theta};\mathbf{0})\big),\\ E_\mathrm{n}(\boldsymbol{\theta};{\boldsymbol{S}}) &= E(\boldsymbol{\theta};{\boldsymbol{S}})/\big(\max_{\boldsymbol{\theta}} E(\boldsymbol{\theta};{\boldsymbol{S}})\big).\label{eq95Entilde} \end{align}\tag{9}\]

4.2 Scattering Parameters Training↩︎

Based on the derivation in Section 4.1, we now design a training procedure to estimate the unknown scattering matrix \({\boldsymbol{S}}\) based on 8 . As mentioned in Section 3.2, matrix \({\boldsymbol{S}}\) is symmetric, indicating that the dof of \({\boldsymbol{S}}\) is  \(N(N+1)/2\). Considering the geometry of the upa of RIS elements, the dof of \({\boldsymbol{S}}\) can be further reduced.

Typically, introducing a third element between two existing element will result in an em interaction characterized by indirect coupling, which is generally exceedingly weak. Moreover, this attenuation intensifies with an increasing number of elements in between, implying that the mutual coupling between two elements separated by multiple intervening elements is negligible. Therefore, in this study, we only consider the mutual coupling between element pairs with 0 or 1 intervening element, while ignoring the rest, as shown in Fig. 7. Furthermore, we note that the unit cells in our RIS prototype exhibit fold symmetry along the \(y\)-axis rather than rotational symmetry. In consideration of such geometric configuration, we can summarize that there are only 8 different values of S-parameters in \({\boldsymbol{S}}\), which we denote as \(S_1,S_2,\dots,S_8\), as depicted in Fig. 7. Additionally, we denote the self-scattering parameter as \(S_0\). Hence, we can approximately express the scattering matrix \({\boldsymbol{S}}\) as \[{\boldsymbol{S}}\approx \sum_{i=0}^8S_i{\boldsymbol{A}}_i,\] where \({\boldsymbol{A}}_i\) is the known support matrix of the value \(S_i\) in \({\boldsymbol{S}}\) according to array configuration. Namely, the entry at the \(j\)-th row and \(k\)-th column is \(A_{i,jk}=1\) for all \(S_{jk}=S_i\); otherwise, \(A_{i,jk}=0\). For example, the support matrix for the self-coupling \(S_0\) is \({\boldsymbol{A}}_0=\mathbf{I}_{N^2}\). This expression suggests the dof of \({\boldsymbol{S}}\) is reduced to 9, which is independent of the RIS size \(N\). In other words, the scattering matrix \({\boldsymbol{S}}\) can be determined by these 9 parameters \(\{S_0,\dots,S_8\}\). In addition, given that the adopted RIS prototype exhibits robust and stable unit cell magnitude and phase performance within \(0-50^\circ\) incidences at [ 27 5 ] GHz [15] , we assume for simplicity in our analysis that these 9 parameters are independent of the incident angle of em waves.

According to 9 , the mutual coupling-aware normalized radiation pattern \(E_\mathrm{n}(\boldsymbol{\theta};{\boldsymbol{S}})\) is a function of \({\boldsymbol{S}}\) and thus is a function of parameters \(\{S_0,\dots,S_8\}\). Suppose there is a normalized radiation pattern obtained through a 3D full-wave simulation based on the actual RIS prototype, denoted as \(\bar{E}_\mathrm{n}(\boldsymbol{\theta})\). We sample this simulated pattern over aod \(\boldsymbol{\theta}\), thus acquiring the discrete version of the simulated pattern as \(\bar{E}_\mathrm{n}(\boldsymbol{\theta}_\ell),\ell=1,\dots,L_{\rm S}\). Here, \(L_{\rm S}\) denotes the number of samples. Based on the RIS radiation prediction formula 9 and the obtained samples of the actual radiation pattern, we can design an optimization problem to estimate the scattering parameters \(\{S_0,\dots,S_8\}\). The objective is to optimize these scattering parameters to minimize the discrepancy between the predicted radiation pattern and the actual radiation pattern. Consequently, a least-squares problem is formulated as follows: \[\begin{align} &\{\hat{S}_0,\!\dots,\!\hat{S}_8\} \notag \\ &=\!\arg\min_{S_0,\dots,S_8}\sum_{\ell=1}^{L_\mathrm{S}}\big(\frac{E(\boldsymbol{\theta}_\ell;\!S_0,\!\dots,\!S_8\!)}{\max\limits_\ell E(\boldsymbol{\theta}_\ell;\!S_0,\!\dots,\!S_8\!)}-\bar{E}_\mathrm{n}(\boldsymbol{\theta}_\ell)\big)^2. \label{eq95min} \end{align}\tag{10}\] Subsequently, the scattering matrix can be recovered as \[\hat{{\boldsymbol{S}}} = \sum_{i=0}^8\hat{S}_i{\boldsymbol{A}}_i.\]

In the following subsection, we present a solution to 10 by reformulating the optimization problem based on the Neumann series expansion.

4.3 Solving 10 via Neumann Series Approximation↩︎

Due to the matrix inversion in 8 and the normalization operation in 10 , the objective function in 10 is non-convex and non-differentiable. Although gradient-based iterative algorithms can be applied using the numerical gradient (e.g., through finite difference approximation), the local-minimum issue is unavoidable, and there is no guarantee that such algorithms can converge to the global minimum. To this end, here we propose a more tractable method to estimate the scattering parameters by reformulating 10 based on the Neumann series expansion.

The Neumann series expansion [52] reads \[\label{eq:Neumann} (\boldsymbol{\Theta}^{-1}-{\boldsymbol{S}})^{-1}\!=\!\Big(\sum_{n=0}^\infty \big(\boldsymbol{\Theta}{\boldsymbol{S}}\big)^n\Big)\boldsymbol{\Theta} \!=\! \boldsymbol{\Theta}\! +\! \boldsymbol{\Theta}{\boldsymbol{S}}\boldsymbol{\Theta}\! +\! \boldsymbol{\Theta}{\boldsymbol{S}}\boldsymbol{\Theta}{\boldsymbol{S}}\boldsymbol{\Theta}\!+\! \cdots,\tag{11}\] if all the eigenvalues of \(\boldsymbol{\Theta}{\boldsymbol{S}}\) are within the unit circle, i.e., \(|\lambda_i(\boldsymbol{\Theta}{\boldsymbol{S}})|<1\), \(\forall\;i=1,\dots,N\).⁴ Therefore, based on 8 and 11 , we have \[\begin{align} {E}(\boldsymbol{\theta}) &= \big|{\boldsymbol{h}}_{\mathrm{out}}^\mathsf{T}(\boldsymbol{\theta})\big(\boldsymbol{\Theta}^{-1}-{\boldsymbol{S}}\big)^{-1}{\boldsymbol{h}}_\mathrm{in}\big|,\\ &\overset{(a)}{\approx} \big|{\boldsymbol{h}}_{\mathrm{out}}^\mathsf{T}(\boldsymbol{\theta})\big(\boldsymbol{\Theta}+\boldsymbol{\Theta}{\boldsymbol{S}}\boldsymbol{\Theta}\big){\boldsymbol{h}}_\mathrm{in}\big|, \label{eq:appro} \end{align}\tag{12}\] where step (a) is obtained by only preserving the first two terms of the Neumann series expansion.

Then, we can reformulate the objective function in 10 as \[\begin{gather} f(S_0,\dots,S_8,\varepsilon) = \Big\|\Big({\boldsymbol{H}}_{\mathrm{out}}^\mathsf{T}\big(\boldsymbol{\Theta}\!+\!\boldsymbol{\Theta}\big(\sum_{i=0}^8S_i{\boldsymbol{A}}_i\big)\boldsymbol{\Theta}\big){\boldsymbol{h}}_\mathrm{in}\Big)\!\\ \odot\!\Big({\boldsymbol{H}}_{\mathrm{out}}^\mathsf{H}\big(\boldsymbol{\Theta}^*\!+\!\boldsymbol{\Theta}^*\big(\sum_{i=0}^8S_i^*{\boldsymbol{A}}_i\big)\boldsymbol{\Theta}^*\big){\boldsymbol{h}}_\mathrm{in}^*\Big)\!-\!\varepsilon\bar{{\boldsymbol{e}}}\Big\|_2,\label{eq:f} \end{gather}\tag{13}\] where \(\varepsilon\in\mathbb{R}^+\) is an auxiliary variable counteracting the unknown normalization scale, \(\odot\) denotes Hadamard (elementwise) product, \({\boldsymbol{H}}_\mathrm{out}\overset{\Delta}{=}[{\boldsymbol{h}}_\mathrm{out}(\boldsymbol{\theta}_1),\dots,{\boldsymbol{h}}_\mathrm{out}(\boldsymbol{\theta}_{L_{\rm S}})]\), and \(\bar{{\boldsymbol{e}}}\overset{\Delta}{=}[\bar{E}^2_\mathrm{n}(\boldsymbol{\theta}_1),\dots,\bar{E}^2_\mathrm{n}(\boldsymbol{\theta}_{L_{\rm S}})]^\mathsf{T}\). Subsequently, we can estimate the scattering parameters by solving the following optimization problem: \[\label{eq:opt2} \begin{align} \min_{S_0,\dots,S_8,\varepsilon}\;&f(S_0,\dots,S_8,\varepsilon),\\ \mathrm{s. t.}\quad &S_0,\dots,S_8\in\mathbb{C},\;\varepsilon\in\mathbb{R}^+. \end{align}\tag{14}\] This is a convex optimization problem, the convexity of which is presented in Appendix. By applying a gradient descent algorithm, the scattering parameters \(S_0,\dots,S_8\) can be estimated. Note that the direction where the function \(f\) decreases fastest is given by \(-\partial f/\partial [S_0^*,\dots,S_8^*,\varepsilon]^\mathsf{T}\) [53].

5 Simulation Results↩︎

This section first conducts the model training procedure as described in Section 4.2 and Section 4.3, and then assesses the accuracy of the trained model via full-wave simulation.

5.1 Model training↩︎

Before conducting the estimation process proposed in Section 4.3, we first perform a full-wave simulation based on a similar transmitting mimo setup to obtain a set of rough estimates of the unknown scattering parameters, which can serve as (i) the initialization of the iterative gradient descent method and (ii) an additional regularization term to the objective function \(f(S_0,\dots,S_8,\varepsilon)\) in 13 to mitigate over-fitting.

5.1.1 Simulation of Transmitting MIMO Antennas↩︎

As mentioned in Section 4, the scattering matrix \({\boldsymbol{S}}\) for the RIS cannot be measured or simulated directly since RIS is excited through incident EM waves instead of direct feeding ports. Nevertheless, from the em coupling perspective, the mutual coupling of the reflective RIS elements and the transmitting mimo antennas should possess a similar tendency. Hence, a similar antenna model with mimo configuration can be simulated using commercial full-wave simulation software such as Ansys HFSS, to obtain an approximated scattering matrix to initialize 14 . We simulate a \(3\times 3\) mimo antenna array with a similar physical layout as RIS using Ansys HFSS, as presented in Fig. 8.⁵

The simulated results of the scattering parameters among such a mimo antenna array are listed in Table 1. To be clear, we denote these simulated mimo scattering parameters as \(\tilde{S}_0,\tilde{S}_1,\dots,\tilde{S}_8\) and the corresponding scattering matrix as \(\tilde{{\boldsymbol{S}}}=\sum_{i=0}^8\tilde{S}_i{\boldsymbol{A}}_i\). By observing the data in Table 1, insights are obtained that the adjacent coupling in the \(x\)-axis direction is stronger than that in the \(y\)-axis direction, since \(|S_1|>|S_3|\). This phenomenon results from the non-rotational symmetry of the antenna structure. To further validate and better understand this difference. The surface current distribution of the mimo antenna configuration is also demonstrated in Fig. 8. Here, the unit in the right-bottom corner is excited, leading to the induction of currents in neighboring units through the coupling effect. Therefore, the EM coupling along the \(x\) and \(y\)-axes can be intuitively observed through the surface current distribution on adjacent unit cells. It is clear that the coupling along the \(x\)-axis results in a stronger coupled surface current than that along the \(y\)-axis, coinciding with the simulated results in Table 1. These observations imply that the mutual coupling in the RIS can also exhibit such a trend.

5.1.2 RIS Model Training↩︎

Now we conduct the estimation of the S-parameters based on the method presented in Section 4.2 and 4.3. Using the simulated \(\tilde{{\boldsymbol{S}}}\) in Table 1, we can easily verify that the condition \(|\lambda_i(\boldsymbol{\Theta}\tilde{{\boldsymbol{S}}})|<1\), \(\forall\;i=1,\dots,N\) holds. Thus, the Neumann series expansion in 11 is valid based on the fact that the true \({\boldsymbol{S}}\) is close to \(\tilde{{\boldsymbol{S}}}\). Besides, Fig. 10 plots the truncation error of this Neumann series expansion. It is shown that by considering only the first two terms, we can obtain an approximation with a truncation error of approximately \(0.018\%\) only, which justifies the accuracy of the approximation in 12 .

Besides justifying the accuracy of the Neumann series-based approximation, the simulated \(\tilde{{\boldsymbol{S}}}\) can further benefit our estimation by providing a regularization penalty. Specifically, we rewrite 14 as \[\label{eq:opt3} \begin{align} \min_{S_0,\dots,S_8,\varepsilon}\;&f(S_0,\dots,S_8,\varepsilon) + \sum_{i=0}^8\zeta_i|S_i-\tilde{S}_i|^2,\\ \mathrm{s. t.}\quad &S_0,\dots,S_8\in\mathbb{C},\;\varepsilon\in\mathbb{R}^+, \end{align}\tag{15}\] which is still a convex problem and \(\zeta_i\) controls the importance/trust of the each regularization term.

Utilizing the simulated data in Table 1 as an initialization, this convex optimization problem can be solved using the gradient descent method. The observations \(\bar{E}_\mathrm{n}(\boldsymbol{\theta}_\ell),\ell=1,2,\dots,L_{\rm S}\), are obtained through a single 3D full-wave simulation of RIS radiation pattern using Ansys HFSS. Here, we sample the 3D pattern uniformly such that \(L_{\rm S}=91\times 91\). This simulation is conducted at frequency of \(\unit[27.5]{GHz}\), where the Tx horn is placed at a distance of \(\unit[0.18]{m}\) to the RIS with incident angle \(\phi_\mathrm{az}=90^\circ,\;\phi_\mathrm{el}=20^\circ\).⁶ The RIS phase shifts and amplitude gains are assigned following the beamforming method in [15]. The overall model training procedure and results are visualized in Fig. 9, and the values of the estimated scattering parameters \(\hat{S}_0,\hat{S}_1,\dots,\hat{S}_8\) are listed in Table 2.

5.2 Simulation Validation↩︎

Since we have learned a scattering matrix \(\hat{{\boldsymbol{S}}}\) based on the simulated RIS radiation pattern in Fig. 9, we now evaluate the accuracy of the estimated \(\hat{{\boldsymbol{S}}}\) as well as model 3 . We conduct another 3 different full-wave simulations of RIS radiation pattern by varying the system setup. Specifically, we change the incident elevation angle of the Tx horn to \(\phi_\mathrm{el}=30^\circ\), \(40^\circ\), and \(50^\circ\), respectively.⁷ The corresponding RIS beamforming is also changing according to the method in [15]. Fig. 11 compares the different theoretical predictions and the full-wave simulation (serving as ground truth) of RIS radiation patterns over these 3 different setups. The evaluated theoretical results include the prediction based on the conventional mutual coupling-unaware model (i.e., \({E}_\mathrm{n}(\boldsymbol{\theta};\mathbf{0})\)), the mutual coupling-aware model using simulated \(\tilde{{\boldsymbol{S}}}\) in the MIMO setup (i.e., \({E}_\mathrm{n}(\boldsymbol{\theta};\tilde{{\boldsymbol{S}}})\)), and the mutual coupling-aware model using trained \(\hat{{\boldsymbol{S}}}\) (i.e., \({E}_\mathrm{n}(\boldsymbol{\theta};\hat{{\boldsymbol{S}}})\)).

Typically, for the normalized radiation pattern of an array, the predicting accuracy of its sidelobe (especially the strongest sidelobe) is a major concern. As demonstrated in Fig. 11, the prediction based on the mutual coupling-aware model using trained \(\hat{{\boldsymbol{S}}}\) (red dashed curve) outperforms the other two theoretical predictions, showing the highest accuracy aligned with the full-wave simulations across the three tested setups. The corresponding prediction errors are listed in Table 3, which are computed as \[\label{eq:PE1} {\rm P}{\rm E}= \frac{\|{\boldsymbol{r}}({\boldsymbol{Q}})-\bar{{\boldsymbol{r}}}\|_2}{\|\bar{{\boldsymbol{r}}}\|_2},\tag{16}\] where \(\bar{{\boldsymbol{r}}}\overset{\Delta}{=}[\bar{E}_\mathrm{n}(\boldsymbol{\theta}_1),\dots,\bar{E}_\mathrm{n}(\boldsymbol{\theta}_{L_{\rm R}})]^\mathsf{T}\in\mathbb{R}^{L_{\rm R}}\), \({\boldsymbol{r}}({\boldsymbol{Q}})\overset{\Delta}{=}[{E}_\mathrm{n}(\boldsymbol{\theta}_1;{\boldsymbol{Q}}),\dots,{E}_\mathrm{n}(\boldsymbol{\theta}_{L_{\rm R}};{\boldsymbol{Q}})]^\mathsf{T}\in\mathbb{R}^{L_{\rm R}}\), with \(L_{\rm R}\) denoting the sample size of these testing radiation patterns and \({\boldsymbol{Q}}\in\{\mathbf{0},\tilde{{\boldsymbol{S}}},\hat{{\boldsymbol{S}}}\}\). Here, we set \(L_{\rm R}=90\) by uniformly sampling over \(-90^\circ<\theta_\mathrm{el}<90^\circ\) with a \(2^\circ\) step size. It is observed that the trained model \(\hat{{\boldsymbol{S}}}\) surpasses the other two models in predicting the RIS radiation pattern with higher accuracy. Particularly, at \(\phi_\mathrm{el}=40^\circ\), there is a reduction of  [ 10 7 ] % in prediction error observed These results reveal that : ( i ) the mutual coupling-aware communication model 3 can correctly describe the mutual coupling effect in RIS , thus presenting a more precise performance in predicting RIS radiation pattern compared to the mutual coupling - unaware model \({\boldsymbol{Q}}=\mathbf{0}\). (ii) the proposed model training method can accurately estimate the scattering matrix \({\boldsymbol{S}}\), as the prediction result using trained \(\hat{{\boldsymbol{S}}}\) significantly and consistently outperforms that using simulated \(\tilde{{\boldsymbol{S}}}\) in MIMO setup. An additional observation is that at certain AoAs, such as \(\phi_\mathrm{el}=40^\circ\), the deviation between the conventional model and the ground truth is much more significant than in the other two cases. This indicates an angle-selective property of the impact of mutual coupling.

5.3 Complexity Analysis↩︎

5.3.1 Resources Consumption of the Full-Wave Simulation↩︎

Our proposed S-parameters estimation method is based on the observations obtained through the full-wave simulation. To evaluate the complexity of this full-wave simulation process, Table 4 summarizes the computational resources consumptions when simulating over various RIS sizes. Here, the processor of the workstation is Intel(R) Xeon(R) Gold 6230 CPU @[ 2 10 ] GHz ( 2 processors ) with installed RAM of GB , and the utilized simulation software is ANSYS HFSS Version 2021 1 0 It can be observed that as the RIS array becomes larger , the consumed computational resources significantly increase , making full - wave simulation on large - scale RISs impossible However , since our scattering model is parameterized into limited amount of S - parameters , it is always feasible to perform full - wave simulations on a smaller RIS as long as all the defined S - parameters ( see Fig 7 ) are involved Nevertheless , this reduction can diminish model robustness , thus indicating a trade - off between the computational complexity and the simulation precision.

5.3.2 Computational Complexity of the Model Training↩︎

Now we evaluate the computational complexity of the proposed model training process. According to the gradient descent principle and the derivative calculation presented in Appendix, we can determine the computational complexity of the proposed scattering parameters estimation method as \(\mathcal{O}(L_\mathrm{S} N^2)\). The complexity of the estimation process increases quadratically with the RIS size and linearly with the number of radiation pattern observation samples. While the size of the RIS is fixed following the previous full-wave simulation step, the number of observation samples is adjustable, posing a potential impact on estimation accuracy. Fig. 12 shows the prediction errors evaluated according to 16 with \({\boldsymbol{Q}}=\hat{{\boldsymbol{S}}}\) versus sample sizes \(L_\mathrm{S}\). It is observed that a greater \(L_\mathrm{S}\) can generally result in a lower prediction error, which also reveals a trade-off between the computational complexity and the estimation accuracy.

6 Experimental Results↩︎

Although the accuracy of the mutual coupling-aware model and the trained scattering matrix have been verified by full-wave simulation in Section 5.2, further validations are performed in this section through experimental measurement on the fabricated RIS prototype described in Section 2.2. The measurement setup is demonstrated in Fig. 13. Here, we use two PE9851B/SF-20 horn antennas to transmit and receive wireless signals, while the fabricated 1-bit quasi-passive RIS is utilized to reflect the em waves. As mentioned, an MCU WROOM-32 module is used to control the beamforming at the RIS, which can be visualized through the LED array at the back surface of the RIS. Based on this setup, the S-parameter between the Tx and Rx horns can be measured using an Anritsu ME7828A vector network analyzer. This real measured Tx-Rx S-parameter is denoted as \(\bar{S}_\mathrm{RT}\). Both the distance between Tx and RIS and the distance between Rx and RIS are set as \(\unit[0.18]{m}\), which are within Rayleigh distance (approximately [ 2 14 ] m ) of the fabricated RIS and form a near - field scenario.

On the other hand, both communication model 2 and 3 (with trained \(\hat{{\boldsymbol{S}}}\)) can be adopted to predict this S-parameter \(S_\mathrm{RT}\). By comparing both theoretical predictions with the measured results, we can assess the precision of model 2 and 3 . Specifically, the S-parameter between Tx and Rx horns can be theoretically calculated as \[\begin{align} \text{\boldsymbol{MC-unaware:} } & S_\mathrm{RT}^\mathrm{conv} = {\boldsymbol{h}}_\mathrm{RI}\boldsymbol{\Theta}{\boldsymbol{h}}_\mathrm{IT}, \\ \text{\boldsymbol{MC-aware:} } & S_\mathrm{RT}^\mathrm{new} = {\boldsymbol{h}}_\mathrm{RI}\big(\boldsymbol{\Theta}^{-1}-\hat{{\boldsymbol{S}}}\big)^{-1}{\boldsymbol{h}}_\mathrm{IT}. \end{align}\] Here, each entry of the wireless channels \({\boldsymbol{h}}_\mathrm{IT}\in\mathbb{C}^{N\times 1}\) and \({\boldsymbol{h}}_\mathrm{RI}\in\mathbb{C}^{N\times 1}\) are computed as \[\begin{align} h_{\mathrm{IT},n} &= \frac{\cos^{q_\mathrm{f}}(\theta_{\mathrm{T},n})}{\|{\boldsymbol{p}}_n-{\boldsymbol{p}}_\mathrm{T}\|_2}e^{-jk_0\|{\boldsymbol{p}}_n-{\boldsymbol{p}}_\mathrm{T}\|_2},\\ h_{\mathrm{RI},n} &= \frac{\cos^{q_\mathrm{e}}(\theta_{\mathrm{e},n}) \cos^{q_\mathrm{f}}(\theta_{\mathrm{R},n})}{\|{\boldsymbol{p}}_n-{\boldsymbol{p}}_\mathrm{R}\|_2}e^{-jk_0\|{\boldsymbol{p}}_n-{\boldsymbol{p}}_\mathrm{R}\|_2}, \end{align}\] where \({\boldsymbol{p}}_\mathrm{T}\) and \({\boldsymbol{p}}_\mathrm{R}\) stand for the positions of the Tx and Rx horns, respectively. Note that these channels \({\boldsymbol{h}}_\mathrm{IT}\) and \({\boldsymbol{h}}_\mathrm{RI}\) have been adjusted to the near-field channel model to fit the measurement setup. Before the final comparison, we process the following minimization to remove the global phase shift and amplitude loss: \[\label{eq95Glo} \{\hat{\alpha},\hat{\beta}\} = \arg\min_{\alpha,\beta}\;\sum_{i=1}^{D}|\alpha e^{j\beta} S_{\mathrm{RT},i}^\mathrm{X}-\bar{S}_{\mathrm{RT},i}|^2,\tag{17}\] where \(i\) indices the prediction or measurement over different system setups and \(\mathrm{X}\in\{\mathrm{conv},\mathrm{new}\}\). Finally, we adopt \(\hat{S}_{\mathrm{RT},i}^\mathrm{X}=\hat{\alpha} e^{j\hat{\beta}} S_{\mathrm{RT},i}^\mathrm{X},\;i=1,2,\dots,\) as the theoretical prediction for both models.

The measured and predicted S-parameters over different setups are plotted in Fig. 14 and the corresponding numerical errors between the measurements and predictions are computed and presented in Table 5. Here, we show the results of 3 different setups of the Tx horn, which correspond to the setups in Fig. 11. For each Tx setup, we test the cases where the Rx horn is placed at the 3 chosen aods \(\theta_\mathrm{el}\) corresponding to the peak values in Fig. 11, aiming to attain higher SNRs. As shown in Table 5, the mutual coupling-aware model consistently offers more accurate predictions than the conventional model when benchmarked against real measured data. Interestingly, for \(\phi_\mathrm{el}=30^\circ\), the mutual coupling-aware model yields only minor improvements over the conventional model. This can be attributed to the minimal impact of mutual coupling in this configuration, where the conventional model already provides accurate predictions of RIS radiation, as evidenced in Table 3. However, for \(\phi_\mathrm{el}=40^\circ\) and \(\phi_\mathrm{el}=50^\circ\), the mutual coupling-aware model exhibits significantly lower prediction errors. This observation further validates the correctness of the scattering matrix-based model 3 and the accuracy of our estimation on \({\boldsymbol{S}}\).

7 Conclusion↩︎

This paper conducts a realistic evaluation of the mutual coupling in RIS-aided communication based on an authentic RIS prototype. Adopting the scattering matrix-based model, we parameterize the RIS scattering matrix by analyzing the geometric configuration of RIS reflecting elements. A practical training approach for these scattering parameters is proposed leveraging a single 3D full-wave simulation of the RIS radiation pattern. Both full-wave simulations and experimental measurements are carried out to verify the accuracy of the trained model. The results of this work validate the correctness of the RIS-aided wireless communication model based on scattering matrices. Furthermore, it offers an effective method to estimate RIS scattering parameters, foreseeing substantial potential in applications such as RIS-aided channel estimation, beamforming, radio localization, etc. Nonetheless, the proposed model training method requires full-wave simulation, which is time-consuming and thus not suitable for online mutual coupling calibration. Future work can explore more complex RIS designs, such as multi-polarization, nonuniform array configurations, and more realistic modeling that considers the frequency and configuration dependencies of mutual coupling.

This appendix demonstrates the convexity of 14 . The following derivations are based on the complex-valued matrix differentiation theory [54]. Considering a complex-value function \(f({\boldsymbol{x}})\): \(\mathbb{C}^{N_x}\rightarrow \mathbb{R}\), the adopted differentiation rule derives complex derivatives by treating \({\boldsymbol{x}}\) and \({\boldsymbol{x}}^*\) as two independent variables. In addition, for \(y\in\mathbb{C}\), \({\boldsymbol{y}}\in\mathbb{C}^{N_y}\), and \({\boldsymbol{x}}\in\mathbb{C}^{N_x}\), we use notations \(\frac{\partial y}{\partial {\boldsymbol{x}}}\in\mathbb{C}^{N_x\times 1}\), \(\frac{\partial {\boldsymbol{y}}}{\partial {\boldsymbol{x}}}\in\mathbb{C}^{N_y\times N_x}\) to denote the first-order partial derivatives.

We first define the unknown variables into the vector form as \({\boldsymbol{s}}\overset{\Delta}{=}[S_0,\dots,S_8]^\mathsf{T}\) and \({\boldsymbol{v}}\overset{\Delta}{=}[{\boldsymbol{s}}^\mathsf{T},\varepsilon]^\mathsf{T}\). The domain of \({\boldsymbol{v}}\) is the manifold \(\mathbb{C}^9\times\mathbb{R}^+\), which is a convex set. Based on 13 , we define \[\begin{gather} g_\ell({\boldsymbol{v}}) = \Big({\boldsymbol{h}}_{\mathrm{out}}^\mathsf{T}(\boldsymbol{\theta}_\ell)\big(\boldsymbol{\Theta}\!+\!\boldsymbol{\Theta}\big(\sum_{i=0}^8S_i{\boldsymbol{A}}_i\big)\boldsymbol{\Theta}\big){\boldsymbol{h}}_\mathrm{in}\Big)\!\\ \Big({\boldsymbol{h}}_{\mathrm{out}}^\mathsf{H}(\boldsymbol{\theta}_\ell)\big(\boldsymbol{\Theta}^*\!+\!\boldsymbol{\Theta}^*\big(\sum_{i=0}^8S_i^*{\boldsymbol{A}}_i\big)\boldsymbol{\Theta}^*\big){\boldsymbol{h}}_\mathrm{in}^*\Big)\!-\!\varepsilon\bar{E}^2_\mathrm{n}(\boldsymbol{\theta}_\ell), \\ \ell=1,\dots,L_{\rm S}. \end{gather}\] Then we have \[\begin{align} \frac{\partial g_\ell({\boldsymbol{v}})}{\partial {\boldsymbol{s}}} &\!=\! {\boldsymbol{a}}_\ell\Big({\boldsymbol{h}}_{\mathrm{out}}^\mathsf{H}(\boldsymbol{\theta}_\ell)\big(\boldsymbol{\Theta}^*\!+\!\boldsymbol{\Theta}^*\big(\sum_{i=0}^8S_i^*{\boldsymbol{A}}_i\big)\boldsymbol{\Theta}^*\big){\boldsymbol{h}}_\mathrm{in}^*\Big),\\ \frac{\partial g_\ell({\boldsymbol{v}})}{\partial {\boldsymbol{s}}^*} &= {\boldsymbol{a}}_\ell^*\Big({\boldsymbol{h}}_{\mathrm{out}}^\mathsf{T}(\boldsymbol{\theta}_\ell)\big(\boldsymbol{\Theta}\!+\!\boldsymbol{\Theta}\big(\sum_{i=0}^8S_i{\boldsymbol{A}}_i\big)\boldsymbol{\Theta}\big){\boldsymbol{h}}_\mathrm{in}\Big),\\ \frac{\partial g_\ell({\boldsymbol{v}})}{\partial \varepsilon} &= \frac{\partial g_\ell({\boldsymbol{v}})}{\partial \varepsilon^*} = \bar{E}^2_\mathrm{n}(\boldsymbol{\theta}_\ell), \end{align}\] where \([{\boldsymbol{a}}_\ell]_i = {\boldsymbol{h}}_{\mathrm{out}}^\mathsf{T}(\boldsymbol{\theta}_\ell)(\boldsymbol{\Theta}\!+\!\boldsymbol{\Theta}{\boldsymbol{A}}_i\boldsymbol{\Theta}){\boldsymbol{h}}_\mathrm{in}\). Therefore, we can derive \[\begin{align} \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial {\boldsymbol{v}}\partial {\boldsymbol{v}}} &\!=\! \begin{bmatrix} \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial {\boldsymbol{s}}\partial{\boldsymbol{s}}} & \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial \varepsilon\partial{\boldsymbol{s}}}\\[2pt] \big(\frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial {\boldsymbol{s}}\partial\varepsilon}\big)^\mathsf{T}& \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial \varepsilon\partial\varepsilon} \end{bmatrix} \!= \mathbf{0}_{10\times 10},\\ \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial {\boldsymbol{v}}^* \partial {\boldsymbol{v}}} &\!=\! \begin{bmatrix} \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial {\boldsymbol{s}}^*\partial{\boldsymbol{s}}} & \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial \varepsilon^*\partial{\boldsymbol{s}}}\\[2pt] \big(\frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial {\boldsymbol{s}}^*\partial\varepsilon}\big)^\mathsf{T}& \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial \varepsilon^*\partial\varepsilon} \end{bmatrix} \!=\! \begin{bmatrix} {\boldsymbol{a}}_\ell{\boldsymbol{a}}_\ell^\mathsf{H}& \mathbf{0}_{9\times 1}\\[3pt] \mathbf{0}_{1\times 9} & 0 \end{bmatrix},\\ \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial {\boldsymbol{v}}\partial {\boldsymbol{v}}^*} &\!=\! \begin{bmatrix} \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial {\boldsymbol{s}}\partial{\boldsymbol{s}}^*} & \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial \varepsilon\partial{\boldsymbol{s}}^*}\\[2pt] \big(\frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial {\boldsymbol{s}}\partial\varepsilon^*}\big)^\mathsf{T}& \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial \varepsilon\partial\varepsilon^*} \end{bmatrix} \!=\! \begin{bmatrix} {\boldsymbol{a}}_\ell^*{\boldsymbol{a}}_\ell^\mathsf{T}& \mathbf{0}_{9\times 1}\\[3pt] \mathbf{0}_{1\times 9} & 0 \end{bmatrix},\\ \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial {\boldsymbol{v}}^* \partial {\boldsymbol{v}}^*} &\!=\! \begin{bmatrix} \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial {\boldsymbol{s}}^*\partial{\boldsymbol{s}}^*} & \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial \varepsilon^*\partial{\boldsymbol{s}}^*}\\[2pt] \big(\frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial {\boldsymbol{s}}^*\partial\varepsilon^*}\big)^\mathsf{T}& \frac{\partial^2 g_\ell({\boldsymbol{v}})}{\partial \varepsilon^*\partial\varepsilon^*} \end{bmatrix} \!= \mathbf{0}_{10\times 10}. \end{align}\] According to [54], the full complex Hessian matrix of \(g_\ell({\boldsymbol{v}})\) is given by \[\begin{align} {\boldsymbol{H}}_{g_\ell({\boldsymbol{v}})} = \mathrm{blkdiag}({\boldsymbol{a}}_\ell{\boldsymbol{a}}_\ell^\mathsf{H},0,{\boldsymbol{a}}_\ell^*{\boldsymbol{a}}_\ell^\mathsf{T},0) \succeq \mathbf{0}. \end{align}\] Thus, \(g_\ell({\boldsymbol{v}})\) is convex, \(\forall\ell=1,\dots,L_{\rm S}\), based on the second-order condition of convexity [55].

Now we consider \(f({\boldsymbol{v}})\) in 13 , which can be written as \(f({\boldsymbol{v}}) = \|{\boldsymbol{g}}({\boldsymbol{v}})\|_2\) with \({\boldsymbol{g}}({\boldsymbol{v}})\overset{\Delta}{=}[g_1({\boldsymbol{v}}),\dots,g_{L_{\rm S}}({\boldsymbol{v}})]^\mathsf{T}\). Since the \(\ell_2\)-norm \(\|\cdot\|_2\) is convex and nondecreasing in each argument \(g_\ell({\boldsymbol{v}})\), along with the fact that each \(g_\ell({\boldsymbol{v}})\) is convex, the function \(f({\boldsymbol{v}})\) is convex according to the vector composition rule [55]. As a conclusion, 14 is a convex optimization problem.

References↩︎