Ziwei Wan
email: ziweiwan@bit.edu.cn

Member, IEEE Zhen Gao
email: gaozhen16@bit.edu.cn

Fellow, IEEE Feifei Gao
email: feifeigao@ieee.org

Fellow, IEEE Marco Di Renzo
email: marco.di-renzo@universite-paris-saclay.fr

Fellow, IEEE Mohamed-Slim Alouini
email: slim.alouini@kaust.edu.sa

M.-S Alouini

Terahertz Massive MIMO with Holographic Reconfigurable Intelligent Surfaces

Keywords: Terahertz communications, reconfigurable intelligent surface, massive MIMO, holographic communications, compressive sensing (CS), channel estimation

We propose a holographic version of a reconfigurable intelligent surface (RIS) and investigate its application to terahertz (THz) massive multiple-input multiple-output systems. Capitalizing on the miniaturization of THz electronic components, RISs can be implemented by densely packing subwavelength unit cells, so as to realize continuous or quasicontinuous apertures and to enable holographic communications. In this paper, in particular, we derive the beam pattern of a holographic RIS. Our analysis reveals that the beam pattern of an ideal holographic RIS can be well approximated by that of an ultra-dense RIS, which has a more practical hardware architecture. In addition, we propose a closed-loop channel estimation (CE) scheme to effectively estimate the broadband channels that characterize THz massive MIMO systems aided by holographic RISs. The proposed CE scheme includes a downlink coarse CE stage and an uplink finer-grained CE stage. The uplink pilot signals are judiciously designed for obtaining good CE performance. Moreover, to reduce the pilot overhead, we introduce a compressive sensing-based CE algorithm, which exploits the dual sparsity of THz MIMO channels in both the angular domain and delay domain. Simulation results demonstrate the superiority of holographic RISs over the non-holographic ones, and the effectiveness of the proposed CE scheme.

I. INTRODUCTION

Over the past few years, the demand for wireless data traffic has increased significantly due to the explosive growth of mobile devices and multimedia applications [START_REF] Elayan | Terahertz band: The last piece of RF spectrum puzzle for communication systems[END_REF], [START_REF] Rappaport | Wireless communications and applications above 100 GHz: Opportunities and challenges for 6G and beyond[END_REF]. To accommodate these demands, the possible use of the terahertz (THz) band has attracted great interest from both industry and academia [START_REF] Elayan | Terahertz band: The last piece of RF spectrum puzzle for communication systems[END_REF]- [START_REF] Yuan | Hybrid beamforming for terahertz multi-carrier systems over frequency selective fading[END_REF]. The THz band can provide more abundant bandwidth (from 0.1 THz to 10 THz), higher data rates (from tens of Gbps to several Tbps), and lower latency (of the order of micro-seconds [START_REF] Elayan | Terahertz band: The last piece of RF spectrum puzzle for communication systems[END_REF]), as compared with the millimeter-wave (mmWave) band. The THz band is considered to be a promising candidate to enable beyond 5G and 6G communications.

In spite of the appealing advantages of the THz band, establishing a reliable transmission link at THz frequencies is a non-trivial task. This is because (i) there exist strong atmospheric attenuation and extremely high free-space losses in the THz band; and (ii) the line-of-sight (LoS) link is very sensitive to blockage effects in the THz band and thus the links are usually intermittent. These disadvantages may negatively affect the communication range and may severely degrade the service coverage of THz communication systems. The deployment of massive [START_REF] Lin | Indoor terahertz communications: How many antenna arrays are needed?[END_REF] or even ultra-massive [START_REF] Sarieddeen | Terahertz-band ultra-massive spatial modulation MIMO[END_REF] multipleinput multiple-output (MIMO) systems in the context of THz communications may provide considerable beamforming gain in order to overcome the mentioned limitations, but it may result in an unaffordable power consumption and may put an overweight burden on the overall communication system design.

Recently, the emerging technology of reconfigurable intelligent surface (RIS) has been proposed and applied to wireless communications in order to enhance the communication performance in various scenarios, such as MIMO communications, physical layer security, unmanned aerial vehicle communications, simultaneous wireless information and power transfer, cognitive radio systems [START_REF] Renzo | Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and road ahead[END_REF]- [START_REF] Qiao | Secure transmission for intelligent reflecting surface-assisted mmWave and terahertz systems[END_REF]. Made of passive and metamaterial-based reconfigurable elements, an RIS can manipulate both the phase and amplitude of the incident electromagnetic (EM) signals so as to reflect them towards the desired directions. More importantly, unlike other transmission technologies such as active relay [START_REF] Renzo | Reconfigurable intelligent surfaces vs. relaying: Differences, similarities, and performance comparison[END_REF], an RIS does not need power-hungry radio frequency chains (RFCs) and power amplifiers, which may be beneficial for developing green and cost-efficient communications. Although the application of RISs to the mmWave band has been investigated recently [START_REF] He | Cascaded channel estimation for large intelligent metasurface assisted massive MIMO[END_REF]- [START_REF] Liu | Deep denoising neural network assisted compressive channel estimation for mmWave intelligent reflecting surfaces[END_REF], the utilization of RIS for THz communications is still at its infancy. [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF] propose a joint active and passive beamforming scheme based on convex optimization to maximize the signal-to-interference-plus-noise ratio (SINR) at the receivers. A similar scenario is considered in [START_REF] Nadeem | Asymptotic max-min SINR analysis of reconfigurable intelligent surface assisted MISO systems[END_REF], where the phases of the RIS that maximize the SINR are computed via the projected gradient ascent algorithm.

A. Prior Work

RIS-aided MIMO communications have attracted lots of research interest lately. The authors of

In [START_REF] Huang | Reconfigurable intelligent surface assisted multiuser MISO systems exploiting deep reinforcement learning[END_REF], the advanced deep reinforcement learning (DRL) method is considered for the beamforming design in RIS-aided multiuser systems. The authors of [START_REF] Najafi | Physicsbased modeling and scalable optimization of large intelligent reflecting surfaces[END_REF] develop a scalable optimization framework for large RISs, based on which a two-stage method, including an offline design and an online optimization, is proposed to optimize the RIS. A practical amplitude and phase shift model for the reflecting element of an RIS is investigated in [START_REF] Abeywickrama | Intelligent reflecting surface: Practical phase shift model and beamforming optimization[END_REF], where the beamforming optimization is conducted by considering the phase-dependent amplitude variation of each reflecting element. In [START_REF] Gradoni | End-to-end mutual coupling aware communication model for reconfigurable intelligent surfaces: An electromagnetic-compliant approach based on mutual impedances[END_REF], the authors introduce a new communication model for RISs that accounts for the mutual coupling and the interplay between the amplitude and phase response of the scattering elements of the RIS.

The designs in [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF]- [START_REF] Gradoni | End-to-end mutual coupling aware communication model for reconfigurable intelligent surfaces: An electromagnetic-compliant approach based on mutual impedances[END_REF] rely on the knowledge of global channel state information (CSI). Since no active elements are used in RISs, channel estimation (CE) is a challenging task and an essential prerequisite in RIS-aided MIMO systems. In [START_REF] Mishra | Channel estimation and low-complexity beamforming design for passive intelligent surface assisted MISO wireless energy transfer[END_REF] and [START_REF] Zheng | Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization[END_REF], the authors propose a CE algorithm for RISaided systems, based on the least square (LS) estimator, for application to frequency-flat fading channels and frequencyselective fading channels, respectively. To further facilitate the CE task when large-size arrays are deployed, compressive sensing (CS)-based CE schemes are investigated in [START_REF] He | Cascaded channel estimation for large intelligent metasurface assisted massive MIMO[END_REF], [START_REF] Wan | Broadband channel estimation for intelligent reflecting surface aided mmWave massive MIMO systems[END_REF], where the inherent sparsity of mmWave channels in the angular domain is exploited in order to reduce the pilot overhead. In [START_REF] Taha | Deep learning for large intelligent surfaces in millimeter wave and massive MIMO systems[END_REF] and [START_REF] Liu | Deep denoising neural network assisted compressive channel estimation for mmWave intelligent reflecting surfaces[END_REF], a new architecture of RIS for application to mmWave band is proposed, where a few active RFCs are available at the RIS in order to learn the channel in real time. Based on this architecture, deep learning (DL) techniques are adopted for CE and beamforming.

As far as the application of RISs to THz communications is concerned, the authors of [START_REF] Ma | Joint channel estimation and data rate maximization for intelligent reflecting surface assisted terahertz MIMO communication systems[END_REF] present an RIS-aided THz MIMO communication system for indoor applications. Joint CE and data rate maximization schemes are proposed in [START_REF] Ma | Joint channel estimation and data rate maximization for intelligent reflecting surface assisted terahertz MIMO communication systems[END_REF] based on both CS and DL techniques. In [START_REF] Qiao | Secure transmission for intelligent reflecting surface-assisted mmWave and terahertz systems[END_REF], the problem of RIS-assisted secure transmission in the THz band is investigated. The phase shifts at the RIS that maximize the secrecy rate are obtained based on convex optimization. Based on the current state of research, we evince that the design of RIS-aided THz communication systems is still at an early development stage.

Recently, the concept of holographic communication has been proposed as a new paradigm shift in MIMO [START_REF] Pizzo | Spatially-stationary model for holographic MIMO small-scale fading[END_REF]- [START_REF] Hu | Beyond massive MIMO: The potential of data transmission with large intelligent surfaces[END_REF] and RIS-aided [START_REF] Huang | Holographic MIMO surfaces for 6G wireless networks: Opportunities, challenges, and trends[END_REF] communications. One of the main features of holographic communications is the integration of very large numbers of tiny and inexpensive antennas or reconfigurable elements into a compact space in order to realize a holographic array with a spatially continuous aperture [START_REF] Pizzo | Spatially-stationary model for holographic MIMO small-scale fading[END_REF]- [START_REF] Huang | Holographic MIMO surfaces for 6G wireless networks: Opportunities, challenges, and trends[END_REF]. This holographic architecture is easier to realize in the THz band thanks to the miniaturization of THz electronic components. In [START_REF] Ma | Joint channel estimation and data rate maximization for intelligent reflecting surface assisted terahertz MIMO communication systems[END_REF], for example, the size of each graphenebased reflecting element is 200 µm × 190 µm at a carrier frequency of 0.22 THz that corresponds to a wavelength λ ≈ 1360 µm. Therefore, the reflecting elements can be spaced more densely than λ/2 so as to form a spatially continuous surface [START_REF] Huang | Holographic MIMO surfaces for 6G wireless networks: Opportunities, challenges, and trends[END_REF], since the resulting surface is homogenizable [START_REF] Renzo | Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and road ahead[END_REF]. This densely spaced or continuous implementation of RISs is referred to as holographic RIS. Channel modeling and data transmission schemes based on active holographic surfaces are investigated in [START_REF] Pizzo | Spatially-stationary model for holographic MIMO small-scale fading[END_REF] and [START_REF] Huang | Holographic MIMO surfaces for 6G wireless networks: Opportunities, challenges, and trends[END_REF], respectively. However, to the best of our knowledge, no current research works have tackled the physical layer transmission design of passive holographic communications, where nearly-passive holographic RISs with spatially continuous apertures are deployed.

B. Paper Contributions

In this work, we focus our attention on the analysis of holographic RISs for application to massive MIMO systems in the THz band. In particular, the main contributions of this paper can be summarized as follows:

• We derive the beam pattern of an RIS made of discrete elements, and propose an angular-domain beamforming framework. By applying Fourier analysis to the reflection coefficients of the elements at the RIS, we prove that the beam pattern of an RIS with discrete elements can be represented as a weighted integral of Dirichlet kernel functions. On this basis, we propose an angulardomain beamforming framework. The weighting factors in the beam pattern are designed and the corresponding reflection coefficients of the RIS are reconstructed via the Fourier transform of the obtained weighting factors. • We generalize the analysis and design to (continuous) holographic RISs. Based on the proposed beamforming framework, we derive and obtain closed-form solutions for the beamforming design in two important cases, i.e., narrow beam steering (NBS) and spatial bandpass filtering (SBF), which play an important role in RISaided communication systems. We further extend these solutions to holographic RISs, in which the elements are closely spaced so as to yield a virtually continuous spatial aperture. The results reveal that the beam pattern of an ideal holographic RIS can be well approximated by an ultra-dense RIS, which has a practical hardware architecture.

• We propose a closed-loop CE scheme to effectively estimate the broadband channels of THz massive MIMO systems based on holographic RISs. The proposed approach consists of downlink and uplink transmissions.

In the downlink transmission, the holographic RIS uses SBF beamforming so that the users can coarsely estimate the range of LoS angles. In the subsequent uplink transmission, the users with similar LoS angles are scheduled into the same group, and the coarsely-estimated LoS angles are exploited to design the uplink pilot signals for the finer-grained uplink CE. To further reduce the uplink pilot overhead, a CS-based CE scheme is introduced, where the dual sparsity of THz MIMO channels in both the angular domain and delay domain is leveraged.

C. Notation

Column vectors and matrices are denoted by lower-and upper-case boldface letters, respectively. (•) * , (•) T , (•) H and (•) † denote the conjugate, transpose, conjugate transpose and the pseudo-inverse, respectively. j = √ -1 is the imaginary I.

II. SYSTEM DESCRIPTION AND CHANNEL MODEL

In this section, we present the system model and the effective baseband channel model of the considered RIS-aided THz massive MIMO system over frequency-selective fading channels. We consider an RIS-aided THz massive MIMO system that operates in time division duplex (TDD) mode, as illustrated in Fig. 1. The base station (BS) and userequipments (UEs) are equipped with half-wavelength spaced uniform planar arrays (UPAs) that consist of

M B = M B x ×M B y and M U = M U x × M U y antennas, respectively, where M B x (M U
x ) and M B y (M U y ) are the numbers of antennas along the azimuth and elevation directions, respectively. An RIS of physical size A x × A y is deployed to enhance the effective coverage of the BS. A low-power RIS controller is connected to the RIS and is operated by the BS through a control link. The BS carries out beamforming optimization and feeds the results to the RIS controller through the control link so that the reflection coefficients of each reflecting element of the RIS can be adjusted accordingly. As shown in Fig. 1(a), we assume that the LoS link between the BS and the UE is blocked by an obstacle or by a human body. Thus, the UE communicates with the BS only via the RIS, which is regarded as a virtual LoS transmission in RIS-aided systems [START_REF] Wan | Broadband channel estimation for intelligent reflecting surface aided mmWave massive MIMO systems[END_REF]. To reduce the power consumption and hardware cost, a hybrid analog-digital architecture is considered at the BS, i.e., there are only N RF M B RFCs at the BS, and each of them is connected to M B antennas through M B phase shifters. In addition, each UE employs analog beamforming, in which only one RFC is connected to M U antennas through M U phase shifters. An orthogonal frequency division multiplexing (OFDM) transmission scheme with K subcarriers and sampling period T s is adopted. The cyclic prefix (CP) of length N CP T s is added before each OFDM symbol to avoid inter-symbol interference. The center-carrier frequency is f c corresponding to a wavelength λ.

Through an appropriate deployment of the RIS, we assume that there exists a LoS path between the BS and the RIS. The LoS angles between the BS and the RIS are assumed to be known in advance based on the location of the RIS [START_REF] Wan | Broadband channel estimation for intelligent reflecting surface aided mmWave massive MIMO systems[END_REF]. In THz channels, the path loss of the non-LoS (NLoS) paths is known to be much larger than that of the LoS paths. Therefore, we neglect the NLoS paths in the channel between the BS and the RIS.

In the following text, we first introduce the physical channel model that is widely considered in previous works on RISaided MIMO systems [START_REF] He | Cascaded channel estimation for large intelligent metasurface assisted massive MIMO[END_REF]- [START_REF] Ma | Joint channel estimation and data rate maximization for intelligent reflecting surface assisted terahertz MIMO communication systems[END_REF]. Then, we introduce an effective baseband channel model by taking into account the beamforming design at the BS, the RIS and the UE.

a) Physical Channel Model with Discrete RISs: Under the assumptions that the RIS has N reflecting elements, the downlink spatial channel G ∈ C N ×M B from the BS to the RIS can be modeled as

G = α DL a R ψ R a H B ψ B , (1)
where α DL is the channel coefficient, ψ B = [ψ azi BS , ψ ele BS ] T and ψ R = [ψ azi RIS , ψ ele RIS ] T are the LoS angle of departure (AoD) and LoS angle of arrival (AoA) of the BS-RIS channel, respectively. ψ B (ψ R ) includes both the azimuth part ψ azi BS (ψ azi RIS ) and elevation part ψ ele BS (ψ ele RIS ), which are assumed to be fixed and known as detailed in previous text. a B ψ B ∈ C M B ×1 and a R ψ R ∈ C N ×1 denote the steering vectors at the BS and the RIS, respectively. a B ψ B is given by

a B ψ B = 1, ..., e -jdUPA(mxψ azi BS +myψ ele BS ) , ..., e -jdUPA[(M B x -1)ψ azi BS +(M B y -1)ψ ele BS ] T , (2)
where d UPA = λ/2 is the element spacing of the UPA at the BS, and

0 ≤ m x ≤ (M B x -1), 0 ≤ m y ≤ (M B y -1)
. a R ψ R can be written by using a similar notation and assumptions.

As for the RIS-UE spatial channel H ∈ C M U ×N , we consider a Rician fading channel model that consists of one LoS path and L NLoS paths, as shown in Fig. 1(b). In particular, we have

H =β DL,LoS a U ν LoS a H R µ LoS + 1 √ LK f L l=1 β DL l a U (ν l ) a H R (µ l ), (3)
where β DL,LoS and β DL l are the channel coefficients of the LoS component and the l-th NLoS component, respectively, K f is the Rician factor that denotes the ratio of the energy between the LoS and NLoS channels, µ LoS and ν LoS are the LoS AoD and the LoS AoA, respectively, µ l and ν l are the NLoS AoD and the NLoS AoA of the l-th NLoS path. The steering vector at the UE a U (•) ∈ C M U ×1 can be formulated similar to [START_REF] Rappaport | Wireless communications and applications above 100 GHz: Opportunities and challenges for 6G and beyond[END_REF].

b) Effective Baseband Channel Model: The effective baseband channel is defined as the inner product between the analog beamforming vector at the transceiver and the steering vector of the physical channels. Assuming that the r-th RFC (1 ≤ r ≤ N RF ) of the BS and the UE use the analog beamforming vectors f r ∈ C M B ×1 and w ∈ C M U ×1 , respectively, and that the RIS uses ϕ ∈ C N ×1 as the reflection coefficients, the effective baseband channel h DL,eff r can be written as

h DL,eff r = w H HΦGf r = α DL g B r ψ B β DL,LoS g µ LoS , ψ R g U ν LoS + 1 √ LK f L l=1 β DL l g µ l , ψ R g U (ν l ) , (4)
where

Φ = diag(ϕ) ∈ C N ×N , and g B r ψ B = a H B ψ B f r , g µ, ψ R = a H R (µ) Φa R ψ R
, and g U (ν) = w H a U (ν) are the beam patterns of the BS, the RIS, and the UE, respectively.

The effective baseband channel model in ( 4) can be generalized to the case of frequency-selective channel. In a frequencyselective channel, in particular, the BS-RIS-UE effective baseband channel in the delay domain can be formulated as

h DL r (τ ) = α DL g B r ψ B BS to RIS h DL,LoS (τ ) + h DL,NLoS (τ ) RIS to UE , (5)
where h DL,LoS (τ ) = β DL,LoS g µ LoS , ψ R g U ν LoS p τ -τ DL,LoS , (6) and

h DL,NLoS (τ ) = 1 √ LK f L l=1 β DL l g µ l , ψ R g U (ν l )p τ -τ DL l . (7)
In ( 6) and [START_REF] Yuan | Hybrid beamforming for terahertz multi-carrier systems over frequency selective fading[END_REF], τ DL,LoS and τ DL l denote the delay offsets corresponding to the LoS path and the l-th NLoS path, respectively, and p(τ ) is the pulse shaping filter function.

It is worth noting that, in contrast with RISs made of discrete elements [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF]- [START_REF] Abeywickrama | Intelligent reflecting surface: Practical phase shift model and beamforming optimization[END_REF], [START_REF] Mishra | Channel estimation and low-complexity beamforming design for passive intelligent surface assisted MISO wireless energy transfer[END_REF]- [START_REF] Ma | Joint channel estimation and data rate maximization for intelligent reflecting surface assisted terahertz MIMO communication systems[END_REF], a holographic RIS is modeled as an array with a spatially continuous aperture [START_REF] Huang | Holographic MIMO surfaces for 6G wireless networks: Opportunities, challenges, and trends[END_REF] (i.e., N → ∞). The physical channels associated with a holographic RIS cannot be represented in terms of the finitedimensional matrices in (1) and [START_REF] Jansen | Diffuse scattering from rough surfaces in THz communication channels[END_REF]. Therefore, in this paper, we utilize the effective baseband channel model in ( 5)- [START_REF] Yuan | Hybrid beamforming for terahertz multi-carrier systems over frequency selective fading[END_REF] to describe the channels associated with a holographic RIS. The beam patterns in (4) are elaborated in detail in the next section.

Since we consider that the system operates in TDD mode, we assume that the channel reciprocity between the uplink and downlink transmissions holds. In particular, the uplink channels can be modeled based on the downlink channels reported in previous text. We omit the details of the uplink channel model for brevity.

Remark: Although it has been reported that the overall MIMO channel model in the THz band is consistent with that in the mmWave band (see, e.g., [START_REF] Yan | A dynamic array-of-subarrays architecture and hybrid precoding algorithms for terahertz wireless communications[END_REF], [START_REF] Yuan | Hybrid beamforming for terahertz multi-carrier systems over frequency selective fading[END_REF], [START_REF] Ma | Joint channel estimation and data rate maximization for intelligent reflecting surface assisted terahertz MIMO communication systems[END_REF]), some distinct differences between the THz and the mmWave bands, as detailed in the following text, may cause different setting of channel parameters and thus deserve more attention. First, the atmospheric attenuation caused by molecular absorption [START_REF] Elayan | Terahertz band: The last piece of RF spectrum puzzle for communication systems[END_REF], [START_REF] Rappaport | Wireless communications and applications above 100 GHz: Opportunities and challenges for 6G and beyond[END_REF] becomes non-negligible in the THz band, thus introducing an extra path loss for link budget. Second, THz signals with an extremely short wavelength are likely to undergo diffuse scattering on the scattering surface [START_REF] Jansen | Diffuse scattering from rough surfaces in THz communication channels[END_REF]. This results in a limited number of effective NLoS paths (e.g., L = 1) and a large Rician factor (e.g. K f = 30 dB). Third, the THz are more vulnerable to the presence of blockages, which may result in weak LoS directive paths. These factors motivate the deployment of RISs to provide strong virtual LoS links for THz communications in the absence of reliable LoS/NLoS links. In addition, due to the large number of reflecting elements, which yield a better controllability of the radio signals and a high beamforming gain, holographic-based RISs provide further advantages compared with their non-holographic counterpart. This is further discussed next.

III. BEAMFORMING DESIGN FOR HOLOGRAPHIC RISS

To complete the formulation of the effective baseband channel model in ( 5)- [START_REF] Yuan | Hybrid beamforming for terahertz multi-carrier systems over frequency selective fading[END_REF], we first derive and analyze the beam pattern of RISs based on discrete planar arrays (DPAs). Then, we extend the results to RISs with spatially continuous apertures, which are referred to as continuous metasurfaces (CMSs). In addition, we derive closed-form beamforming solutions in two important cases, i.e., narrow beam steering and spatial bandpass filtering, which both play an essential role in the proposed CE scheme for RIS-aided THz massive MIMO systems.

A. RISs Based on Discrete Planar Arrays

As shown in Fig. 2(a) and (b), a DPA-based RIS placed on the x-y plan consists of numerous evenly-spaced reflecting elements. The number of elements along the xand y-directions is N x ∈ Z and N y ∈ Z, respectively. The distance between two adjacent elements is d, and we assume d ≤ λ/2 due to the Nyquist sampling theorem. We define the total physical size of a DPA-based RIS as A x × A y with A x = N x d and A y = N y d, and assume that A x and A y remain unchanged unless stated otherwise. Let (x m , y n ) be the coordinate of the (m, n)-th reflecting element, then we have

x m = (m -1)d, y n = (n -1)d, (8)
where

1 ≤ m ≤ N x , 1 ≤ n ≤ N y .
Assume that a narrowband reference signal "1" impinges on the RIS with an azimuth AoA θ azi in ∈ [0, 2π) and an elevation AoA θ ele in ∈ [0, π/2] (defined in Fig. 2(a)). The phase-difference of the incident signal at the (m, n)-th element (compared to the reference point (0, 0)) can be written as

a (x m , y n ; ψ in ) = e j 2π λ p T r = e j 2π λ (xm cos θ azi in sin θ ele in +yn sin θ azi in sin θ ele in ) = e j(xmψ azi in +ynψ ele in ) , (9)
where p = [x n , y m , 0] T is the position vector of the (m, n)th element, r = [cos θ azi in sin θ ele in , sin θ azi in sin θ ele in , cos θ ele in ] T is the vector of the incident direction, and λ sin θ azi out sin θ ele out . The definitions of ψ in and ψ out simplify the notation. The received signal along the direction of observation (i.e., reflection) ψ out after the incident signal is reflected by the RIS is denoted by g(ψ out , ψ in ), which is given by the superposition of the signals reflected by all the individual elements of the RIS

ψ in ∆ = ψ azi in , ψ ele in T with ψ azi in ∆ = 2π λ cos θ azi in sin θ ele in , ψ ele in ∆ = 2π λ sin θ azi in sin θ ele in is a 2-tuple variable
g (ψ out , ψ in ) = Nx n=1 Ny m=1 a * (x m , y n ; ψ out )Φ (m, n) a (x m , y n ; ψ in ) = Nx n=1 Ny m=1 Φ (m, n) e -j(xm∆x+yn∆y) , ( 10
)
where ∆ x = ψ azi out -ψ azi in , ∆ y = ψ ele out -ψ ele in , and Φ (m, n) ∈ C is the reflection coefficient of the (m, n)-th element of the RIS, whose amplitude and phase are software-programmable via the RIS controller.

We refer to g (ψ out , ψ in ) as the beam pattern of the DPAbased RIS, which is consistent with the definition in (4). The amplitude of g (ψ out , ψ in ) can be used to evaluate the intensity of the signal along the direction ψ out after reflection from the RIS. In general terms, Φ (m, n) can be treated as a two-dimensional discrete signal in the spatial domain 1 We assume that ψ in = ψ azi in , ψ ele in T is known and is treated as a constant in the beamforming design. This follows from the assumption that the LoS angles of the BS-RIS channel are fixed and known.

2π d ω (k, l)Ξ Nx [d (k -∆ x )] Ξ Ny [d (l -∆ y )] dkdl, (13)
with spatial sampling period equal to d in both the azimuth and elevation directions. By using the discrete-time Fourier transform (DTFT), Φ (m, n) can be represented as

Φ (m, n) = d 2π 2 2π d 2π d ω (k, l) e j(dmk+dnl) dkdl, (11)
where

ω (k, l) = Nx m=1 Ny n=1 Φ (m, n)e -j(dmk+dnl) (12)
is a 2-dimensional periodic function whose period is 2π d with respect to both k and l, and 2π d (•)d(•) denotes the integral in an arbitrary interval of length 2π d . Substituting ( 11) into (10), we obtain [START_REF] Li | Reconfigurable intelligent surface assisted UAV communication: Joint trajectory design and passive beamforming[END_REF], as shown at the top of the next page, where

ω (k, l) = ω (k, l) e j Nx +1 2 dk e j Ny +1 2 dl . (14)
Equation ( 13) reveals that the beam pattern g (ψ out , ψ in ) can be formulated as a weighted integral of Dirichlet kernel functions Ξ(•) [START_REF] Lee | Channel estimation via orthogonal matching pursuit for hybrid MIMO systems in millimeter wave communications[END_REF] whose weighting factors are ω (k, l). This is a generalization of the far-field results in [START_REF] Tang | Wireless communications with reconfigurable intelligent surface: Path loss modeling and experimental measurement[END_REF]. Equations ( 11)-( 13) shed light on the design of angular-domain beamforming for DPA-based RIS. In particular, the proposed design is based on two steps: (i) first, ω (k, l) in ( 13) is optimized in order to design g (ψ out , ψ in ) that corresponds to the desired direction of observation; and (ii) then, the corresponding reflection coefficients Φ (m, n) are reconstructed via [START_REF] Wu | Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network[END_REF]. In the next sub-section, we further explain this procedure by providing the optimal beamforming for two specific cases.

B. Beamforming Design for DPA-Based RIS

In this sub-section, we describe the following two use cases for beamforming design based on the proposed approach. a) Narrow beam steering. Given the desired beamforming direction ψ opt = ψ azi opt , ψ ele opt T , the target of NBS is to maximize |g (ψ out , ψ in ) | for ψ out = ψ opt , and to minimize (null)

|g (ψ out , ψ in ) | for ψ out = ψ opt . Since |Ξ Nx [d (k -∆ x )]| and Ξ Ny [d (l -∆ y )] attain their maximum if and only if k = 2πu d + ∆ x , u ∈ Z and l = 2πv d + ∆ y , v ∈ Z
, respectively, we can impose the following design for the angular-domain coefficients of NBS [START_REF] Xu | Resource allocation for IRS-assisted full-duplex cognitive radio systems[END_REF] where

ω NBS (k, l; ψ opt ) ∝ u∈Z δ(k -k opt,u ) v∈Z δ(l -l opt,v ),
k opt,u = 2πu d + ψ azi opt -ψ azi in and l opt,v = 2πv d + ψ ele opt -ψ ele in .
It can be readily verified that (15) ensures the required periodicity of 2π d with respect to both k and l. The corresponding beam pattern for NBS can be obtained by substituting ( 15) into ( 13) and ( 14), which yields

g NBS (ψ out , ψ in ; ψ opt ) ∝ e j d-Ax 2 ∆x e j d-Ay 2 ∆y × Ξ Nx [d (k opt -∆ x )] Ξ Ny [d (l opt -∆ y )] , (16)
where k opt ∆ = k opt,0 and l opt ∆ = l opt,0 , since in [START_REF] Li | Reconfigurable intelligent surface assisted UAV communication: Joint trajectory design and passive beamforming[END_REF] we choose the integral intervals containing k opt,0 and l opt,0 for k and l, respectively. From [START_REF] Wu | Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network[END_REF], the reflection coefficients for NBS can be formulated as follows Φ NBS (m, n; ψ opt ) ∝ e j(dmkopt+dnlopt) . ( 17) This NBS design results in a narrowest beam pattern towards a certain transmission direction ψ opt , and thus provides the optimal beamforming gain in the point-to-point communications.

b) Spatial bandpass filtering. As far as SBF is concerned, we aim to design |g (ψ out , ψ in ) | to be quasi-constant for ψ azi min ≤ ψ azi out ≤ θ azi max and ψ ele min ≤ ψ ele out ≤ θ ele max , and to be almost zero, i.e., |g (ψ out , ψ in ) | ≈ 0 otherwise, where ψ min = ψ azi min , ψ ele min T and ψ max = ψ azi max , ψ ele max T are referred to as the cut-off angles. To clearly explain the design of SBF, we first decompose the values of k in one period as follows

k min ≤ k ≤ k max , (18)
a < k < k min or k max < k ≤ a + 2π d , (19)
where k min = ψ azi min -ψ azi in , k max = ψ azi max -ψ azi in , and a ∈ R is an arbitrary number that satisfies the constraints a < k min and a + 2π d > k max . Since we assume d ≤ λ/2, the existence of a is guaranteed. Similarly, the values of l in one period can be decomposed as follows

l min ≤ l ≤ l max , ( 20
) b < l < l min or l max < l ≤ b + 2π d , (21)
where l min = ψ ele min -ψ ele in , l max = ψ ele max -ψ ele in , and b ∈ R is an arbitrary number that satisfies the constraints b < l min and b + 2π d > l max . The angular-domain coefficients for SBF (in one period of k and l) can be designed as follows

ω SBF (k, l; ψ min , ψ max ) ∝ e -j Nx+1
2 dk e -j Ny +1 2 dl , 0, for [START_REF] Qian | Beamforming through reconfigurable intelligent surfaces in single-user MIMO systems: SNR distribution and scaling laws in the presence of channel fading and phase noise[END_REF] and [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF] for [START_REF] Danufane | On the path-loss of reconfigurable intelligent surfaces: An approach based on Green's Theorem applied to vector fields[END_REF] and [START_REF] Nadeem | Asymptotic max-min SINR analysis of reconfigurable intelligent surface assisted MISO systems[END_REF] .

( 22) Substituting ( 22) into ( 13) and ( 14), and choosing the intervals of integration equal to a, a + 2π

Ξ Nx [d (k -∆ x )]Ξ Ny [d (l -∆ y )] dkdl. (23)
The corresponding reflection coefficients can be calculated from [START_REF] Wu | Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network[END_REF], which yields

where m = m -Nx+1 2 and n = n -Ny+1 2 . To validate the proposed designs, we show two realizations of |g NBS (ψ out , ψ in ; ψ opt ) | and |g SBF (ψ out , ψ in ; ψ min , ψ max ) | with normalized amplitudes in Fig. 3 and Fig. 4, respectively. It can be observed that the obtained beam patterns based on the proposed beamforming framework well fulfill the desired design. It is noteworthy that the reflection coefficients in [START_REF] Zhou | Robust beamforming design for intelligent reflecting surface aided MISO communication systems[END_REF] and [START_REF] Abeywickrama | Intelligent reflecting surface: Practical phase shift model and beamforming optimization[END_REF] are given in closed-form and are physically realizable. By direct inspection of Fig. 3 and Fig. 4, we observe that (i) the NBS design criterion allows one to obtain a small focused region, which can be useful for beamforming applications; and (ii) the SBF design criterion allows one to obtain a wide focused region, which can be useful for broadcasting applications.

C. Extension to Holographic RISs

As far as the design of DPA-based RISs is concerned, we have considered an element spacing equal to d = λ/2, which is a critical spacing based on the Nyquist sampling theorem. This assumption has been, implicitly or explicitly, adopted in previous research works [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF]- [START_REF] Abeywickrama | Intelligent reflecting surface: Practical phase shift model and beamforming optimization[END_REF], [START_REF] Mishra | Channel estimation and low-complexity beamforming design for passive intelligent surface assisted MISO wireless energy transfer[END_REF]- [START_REF] Ma | Joint channel estimation and data rate maximization for intelligent reflecting surface assisted terahertz MIMO communication systems[END_REF]. However, the design and optimization of RISs based on elements spaced at the critical distance has some inherent drawbacks. In particular, (i) due to the periodicity and the non-negligible sidelobes of the Dirichlet kernel functions, a power leakage phenomenon [START_REF] Wan | Compressive sensing based channel estimation for millimeter-wave full-dimensional MIMO with lens-array[END_REF] is usually observed, which may result in inter-beam interference; (ii) the energy received and reflected by the RIS highly depends on its effective reflection area [START_REF] Tang | Wireless communications with reconfigurable intelligent surface: Path loss modeling and experimental measurement[END_REF]. More specifically, the use of critically-spaced RISs usually degrades the effective reflection area, which may result in a reduced energy efficiency.

Since the physical size of an RIS is limited by several practical factors, an effective solution to overcome the just mentioned drawbacks is to increase the number of reflecting elements and to reduce their spacing (d < λ/2) while keeping A x and A y unchanged, as illustrated in Fig. 2(b). A DPAbased RIS whose elements have an element spacing d < λ/2 is referred to as ultra-dense RIS. It is worth mentioning that the analysis reported in the previous sub-section can be applied to critically-spaced and ultra-dense RISs by appropriately choosing the value of d. In Figs. 5 and6, to elucidate these aspects, we report the beam patterns of DPA-based RISs for different values of d (only the beam patterns along the azimuth direction, i.e., with respect to ψ azi out , are reported for ease of illustration). It can be observed that the beam patterns of an ultra-dense RIS with d = {λ/4, λ/8} are quite similar to those of a critically-spaced RIS (d = λ/2), but the power leakage is suppressed. By increasing the number of reflecting elements, more importantly, the effective reflection area of an ultra-dense RIS is expected to increase as compared with critically-spaced RIS, and thus a larger fraction of the energy of the incident EM signal can be steered towards the desired direction. This point is elaborated in detail in Section V, where the numerical results are presented.

The performance improvement from critically-spaced RISs to ultra-dense RISs naturally motivates us to ask: what is the performance of DPA-based RISs in the asymptotic regime d → 0 and N x , N y → ∞ (while keeping N x d = A x and N y d = A y fixed)? An RIS that is obtained by letting d → 0 is referred to as CMS. In the following text, we show that the beam patterns and the design of the reflection coefficients of a CMS can be obtained by extending the analysis of DPA-based RISs. Consider the CMS illustrated in Fig. 2(c), where each point (x, y), x ∈ [0, A x ], y ∈ [0, A y ] is capable of manipulating the phase and amplitude of the incident signal. The reflection coefficients of a CMS, which are denoted by Φ (x, y), are continuously distributed within [0, A x ] × [0, A y ]. Based on the properties of the Fourier transform, the DTFT in [START_REF] Wu | Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network[END_REF] and ( 12) tends to the continuous-time Fourier transformation (CTFT).

In particular, we have

Φ (x, y) ∝ lim d→0 Φ (m, n) | dm=x,dn=y , x ∈ [0, A x ] , y ∈ [0, A y ] . ( 25
)
The corresponding beam pattern of a CMS, g (ψ out , ψ in ) can be calculated as

g (ψ out , ψ in ) ∝ lim d→0 g (ψ out , ψ in ) . ( 26
)
In this sub-section, the parameters related to CMSs are topmarked with ∼ (•) in order to distinguish them from DPA-based RISs. We note, in particular, that the constraints N x d = A x and N y d = A y need to be inherently enforced in [START_REF] Gradoni | End-to-end mutual coupling aware communication model for reconfigurable intelligent surfaces: An electromagnetic-compliant approach based on mutual impedances[END_REF] and [START_REF] Mishra | Channel estimation and low-complexity beamforming design for passive intelligent surface assisted MISO wireless energy transfer[END_REF]. Based on [START_REF] Gradoni | End-to-end mutual coupling aware communication model for reconfigurable intelligent surfaces: An electromagnetic-compliant approach based on mutual impedances[END_REF] and [START_REF] Mishra | Channel estimation and low-complexity beamforming design for passive intelligent surface assisted MISO wireless energy transfer[END_REF], the next two corollaries report the beamforming design of NBS and SBF for application to CMSs.

Corollary 1. Consider a CMS based on the NBS-based beamforming design. The reflection coefficients and the beam pattern, which are denoted by ΦNBS (x, y; ψ opt ) and gNBS (ψ out , ψ in ; ψ opt ), respectively, can be formulated as

ΦNBS (x, y; ψ opt ) ∝ e j(xkopt+ylopt) , (27)
gNBS (ψ out , ψ in ; ψ opt ) ∝ e -j Ax 2 ∆x e -j Ay 2 ∆y

× sinc A x 2 (k opt -∆ x ) sinc A y 2 (l opt -∆ y ) . ( 28
)
Proof. See Appendix A.

where x = x -A x /2 and ȳ = y -A y /2. The corresponding beam pattern, gSBF (ψ out , ψ in ; ψ min , ψ max ), can be formulated as follows

gSBF (ψ out , ψ in ; ψ min , ψ max ) ∝ e -j Ax 2 ∆x e -j Ay 2 ∆y × lmax lmin kmax kmin sinc A x 2 (k -∆ x ) sinc A y 2 (l -∆ y ) dkdl. (30)
Proof. The proof of Corollary 2 is similar to that of Corollary 1 in Appendix A.

In Figs. 5 and6, we report the beam patterns of a CMS and compare them with those of a critically-spaced RIS and an ultra-dense RIS. As far as CMSs are concerned, we observe that the Dirichlet kernel functions that characterize the beam patterns of critically-spaced and ultra-dense RISs are replaced by the "sinc" functions, which have relatively small sidelobes and no periodicity compared with Dirichlet kernel functions. This is in agreement with the far-field results in [START_REF] Danufane | On the path-loss of reconfigurable intelligent surfaces: An approach based on Green's Theorem applied to vector fields[END_REF].

Although CMSs are not realizable in practice, since it is not possible to build surfaces with an infinite number of reflecting elements, Figs. 5 and6 show that the beam pattern of a CMS can be well approximated by an ultra-dense RIS, which is a practical extension of conventional critically-spaced RISs, especially in the THz band. Based on the examples illustrated in Figs. 5 and6, we observe that a spacing d ≤ λ/4 makes the beam patterns of ultra-dense RISs and CMSs almost indistinguishable from each other (see also Fig. 6 in [START_REF] Danufane | On the path-loss of reconfigurable intelligent surfaces: An approach based on Green's Theorem applied to vector fields[END_REF]). In this context, a holographic RIS can be defined as a CMS in theory and as an ultra-dense RIS in practice.

As illustrated in Fig. 5, the minimum width of the mainlobe of the NBS beam pattern is mainly determined by the physical size, i.e., A x and A y , of the surface. We refer to this minimum width as spatial resolution. This parameter characterizes the ability of RISs of distinguishing different UEs that are closely located. Based on Fig. 5, we observe that the spatial resolution is the same for the three implementations, i.e., critically-spaced RISs, ultra-dense RISs, and CMSs, as long as their physical size is kept unchanged [START_REF] Rajatheva | White paper on broadband connectivity in 6G [White paper[END_REF].

IV. CLOSED-LOOP CHANNEL ESTIMATION SCHEME

In this section, we investigate the CE problem for the holographic RIS described in the previous section and for application to THz massive MIMO systems, as illustrated in Fig. 1. The transmission frame structure is illustrated in Fig. 7. During the whole CE stage, short-length OFDM symbols, which are also referred to as unique words (UWs) [START_REF] Huemer | Classical and Bayesian linear data estimators for unique word OFDM[END_REF], are transmitted as pilot signals. Each UW consists of N CP subcarriers in the frequency domain and has length (duration) equal to 2N CP T s (half of which constitutes the mainbody of N CP -point OFDM symbol while the rest constitutes the N CPpoint CP) in the time domain, as elaborated in Fig. 7. For each RFC at the BS, the beamforming is assumed to be the NBS in order to obtain the high beamforming gain, and it is obtained by assuming that the LoS direction of the BS-RIS channel is known. In particular, the beamforming design of each RFC at the BS is2

Φ B r (m, n) = Φ B NBS m, n; ψ B , ∀r = 1, ..., N RF . (31)
The specific expression of Φ B NBS is obtained from [START_REF] Zhou | Robust beamforming design for intelligent reflecting surface aided MISO communication systems[END_REF]. It is worth noting that the constraint of constant modulus in the hybrid analog-digital architecture is implicitly satisfied when the NBS beamforming is considered. The signals of all the RFCs are added together and a single data stream is obtained. Thus, by dropping the index r in (5), we re-write it as

h DL (τ ) = α DL G B h DL,LoS (τ ) + h DL,NLoS (τ ) , (32)
where G B = g B NBS ψ B ; ψ B is the beam pattern corresponding to Φ B NBS m, n; ψ B . The proposed CS scheme consists of two phases, which are applied to the downlink and uplink transmissions. The two phases are described in the following two sub-sections, respectively.

A. Downlink CE Stage

During the downlink CE stage, the UE coarsely estimates µ LoS and ν LoS by constraining their values within smaller ranges. First, the whole range of AoD at the RIS is divided into G x groups along the azimuth direction and into G y groups along the elevation direction. The range of azimuth-AoD in the g x -th azimuth group, g x = 1, ..., G x , is [ψ azi min,gx , ψ azi max,gx ] with

ψ azi min,gx = 2π λ -1 + 2 (g x -1) G x , ψ azi max,gx = 2π λ -1 + 2g x G x - λ A x , (33)
where λ Ax is the azimuth resolution of the RIS. Without loss of generality, a gap equal to λ Ax between two adjacent groups is assumed. Similarly, the range of elevation-AoD in the g y -th elevation group, g y = 1, ..., G y , is [ψ ele min,gy , ψ ele max,gy ], where

ψ ele min,gy = 2π λ -1 + 2 (g y -1) G y , ψ ele max,gy = 2π λ -1 + 2g y G y - λ A y . (34)
Then, we aim to find the group that contains the LoS AoD of the RIS-UE channel, µ LoS . In particular, the downlink CE stage can be formulated as follows Find (ĝ x , ĝy ) = (g x , g y ) ψ azi min,gx ≤ µ azi LoS ≤ ψ azi max,gx and ψ ele min,gy ≤ µ ele LoS ≤ ψ ele max,gy . [START_REF] Rajatheva | White paper on broadband connectivity in 6G [White paper[END_REF] To this end, for the (g x , g y )-th group, the RIS uses the SBF beamforming with the cut-off angles given in [START_REF] Qiao | Secure transmission for intelligent reflecting surface-assisted mmWave and terahertz systems[END_REF] and [START_REF] Pizzo | Spatially-stationary model for holographic MIMO small-scale fading[END_REF], that is Φ(x, y) = ΦSBF (x, y; ψ min,gx,gy , ψ max,gx,gy ),

where ψ min,gx,gy = ψ azi min,gx , ψ ele min,gy T , ψ min,gx,gy = ψ azi max,gx , ψ ele max,gy T , and the specific expression of ( 36) is obtained from (24) (for the DPA-based RIS) or from (29) (for the CMS).

At the UE, given that the dimension of the antenna array at the UE is relatively small, the NBS beamforming towards M U desired directions is employed to coarsely estimate ν LoS .

Specifically, for the (n x , n y )-th desired direction, 1 ≤ n x ≤ M U

x , 1 ≤ n y ≤ M U y , the UE uses the beamforming design

Φ U (m, n) = Φ U NBS m, n; ψ opt,nx,ny , (37)
where ψ opt,nx,ny = 2π λ -1 + 2(nx-1)

M U x , -1 + 2(ny-1) M U y T
. In an OFDM-based system, the effective baseband channel in the delay domain [START_REF] Ma | Joint channel estimation and data rate maximization for intelligent reflecting surface assisted terahertz MIMO communication systems[END_REF] can be transformed to the frequency domain as follows

h DL,Fd k = 1 √ N CP NCP-1 d=0 h DL (dT s ) e -j 2πd N CP k = γ DL,LoS k g µ LoS , ψ R g U ν LoS + 1 √ LK f L l=1 γ DL l,k g µ l , ψ R g U (ν l ) , ( 38
)
where

1 ≤ k ≤ N CP , h DL,Fd k is the frequency- domain (Fd) channel in the k-th subcarrier, γ DL,LoS k = α DL β DL,LoS G B √ NCP NCP-1 d=0 p dT s -τ DL,LoS e -j 2πd N CP k , and γ DL l,k = α DL β DL l G B √ NCP NCP-1 d=0 p dT s -τ DL l e -j 2πd N CP
k . Substituting the beamforming designs ( 36) and ( 37) into ( 38), the received pilot signal at the UE in the k-th subcarrier can be expressed as

y k,
≤ g x ≤ G x , 1 ≤ g y ≤ G y , 1 ≤ n x ≤ M U
x and 1 ≤ n y ≤ M U y in successive UWs, the UE conducts a search to decide which groups µ LoS and ν LoS belong to, i.e., (ĝ x , ĝy , nx , ny ) = arg max (gx,gy,nx,ny) NCP k=1 y k,gx,gy,nx,ny . [START_REF] Huemer | Classical and Bayesian linear data estimators for unique word OFDM[END_REF] Considering the energy focusing property of NBS (as shown in Fig. 3) and the bandpass property of SBF (as shown in Fig. 4), we can expect that ν LoS ≈ ψ opt,nx,ny , ψ azi min,ĝx,ĝy ≤ µ azi LoS ≤ ψ azi max,ĝx,ĝy , and ψ ele min,ĝx,ĝy ≤ µ ele LoS ≤ ψ ele max,ĝx,ĝy . In other words, after the downlink CE stage, the range of possible values for µ LoS is narrowed down to ψ azi max,ĝxψ azi min,ĝx = 2π 2: Set the beamforming design for each RFC at the BS as in (31); 3: for the (g x , g y )-th group at the RIS,

1 ≤ g x ≤ G x , 1 ≤ g y ≤ G y , do 4:
Set the beamforming design at the RIS as in (36); 5:

for the (n x , n y )-th group at the UE,

1 ≤ n x ≤ M U x , 1 ≤ n y ≤ M U y , do 6:
Set the beamforming design at each UE as in (37); 7:

The BS broadcasts one UW as pilot signals;

8:

Each UE receives the pilot signals y k,gx,gy,nx,ny as in (39); 9: end for 10: end for 11: Each UE obtains its index of optimal groups (ĝ x , ĝy , nx , ny ) via [START_REF] Huemer | Classical and Bayesian linear data estimators for unique word OFDM[END_REF]; direction. This significantly reduces the search space of the uplink finer-grained CE discussed in the next sub-section. The proposed downlink CE scheme is summarized in Algorithm 1.

B. Uplink CE stage

After the downlink CE stage, each UE obtains the indices of the optimal groups (ĝ x , ĝy , nx , ny ), and this information is fed back to the BS via the control links for UE scheduling. Specifically, the BS schedules the UEs having the same (ĝ x , ĝy ) into the same uplink group. Each scheduled group performs the finer-grained uplink CE and the subsequent payload data transmission. In this sub-section, therefore, we consider a generic group as follows G ĝx,ĝy = All the UEs that have the same (ĝ x , ĝy ) obtained from the donwlink CE , [START_REF] Tang | Wireless communications with reconfigurable intelligent surface: Path loss modeling and experimental measurement[END_REF] and we assume card(G ĝx,ĝy ) = N UE . During the uplink CE and payload data transmission stages, the beamforming design at the u-th UE (1 ≤ u ≤ N UE ) in G ĝx,ĝy is chosen as follows

Φ U u (m, n) = Φ U NBS m, n; ψ opt,n u x ,n u y , ( 42
)
where {n u x , nu y } is the index of the optimal group obtained by the u-th UE during the downlink CE stage.

As mentioned in the previous sub-section, after the downlink CE stage, the LoS angle of the u-th UE in G ĝx,ĝy , i.e., µ LoS u , is coarsely estimated and its value is confined to the range between ψ min,ĝx,ĝy = ψ azi min,ĝx , ψ ele min,ĝy T and ψ max,ĝx,ĝy = ψ azi max,ĝx , ψ ele max,ĝy

T

. Hence, during the uplink CE stage, we only need to search within this range in order to determine the LoS angle with a finer resolution. Given that the spatial resolution of the considered RIS is λ/A x along the azimuth direction and λ/A y along the elevation direction, we define the search space for the uplink CE stage as

ζ azi bx = ψ azi min,ĝx + 2π λ (b x -1)λ A x , 1 ≤ b x ≤ B x , ( 43
)
ζ ele by = ψ ele min,ĝy + 2π λ (b y -1)λ A y , 1 ≤ b y ≤ B y , ( 44
)
where . Instead of executing an exhaustive beam scanning over the B x B y directions in ( 43) and [START_REF] González-Coma | Channel estimation and hybrid precoding for frequency selective multiuser mmWave MIMO systems[END_REF], during the i-th time slot of the uplink CE stage, the RIS employs the overlapped NBS beamforming towards all B x B y directions with different random phases. This can be formulated as follows

B x ∆ = λ 2π ψ azi max,
Φ(x, y) = 1 B x B y Bx bx=1
By by=1 e jθ i,bx ,by ΦNBS x, y; ζ bx,by , ( 46) where θ i,bx,by ∼ U[0, 2π).

By using the beamforming designs in ( 42) and ( 46), and by capitalizing on the channel reciprocity between the downlink and uplink transmissions, the uplink effective baseband channel related to the u-th UE in the i-th time slot can be formulated as Further, h UL,NLoS i,u (τ ) in ( 47) can be written as follows

h UL i,u (τ ) = h UL,LoS i,u (τ ) + h UL,NLoS i,u (τ ), (47)
h UL,NLoS i,u (τ ) = w T i L l=1 h UL,NLoS u,l p τ -τ UL u,l , (50)
where h UL,NLoS u,l is the effective angular-domain channel that corresponds to the l-th NLoS path. The formulation of h UL,NLoS u,l is similar to (49) and thus it is omitted for brevity. In order to simultaneously perform the uplink CE for all the UEs in G ĝx,ĝy and to avoid the interference of different UEs' pilot signals at the BS, we consider the set of dedicated subcarriers (DSCs) for the u-th UE as follows

UW i

K u = {k u,n |n = 1, ...., N used } , (51)
where Taking into account the DSC allocation above, the uplink channel in the delay domain [START_REF]Attenuation by atmospheric gases[END_REF] can be transformed to the frequency domain, which is similar to [START_REF] Lee | Channel estimation via orthogonal matching pursuit for hybrid MIMO systems in millimeter wave communications[END_REF]. In particular, the frequency-domain representation of the channel related to the u-th UE corresponding to its n-th DSC in the i-th time slot can be written as

N used ∆ = N CP /N UE (N used ∈ Z without loss of generality), K u ⊆ {1, ..., N CP }, card(K u ) = N used
h UL,Fd i,u,n = 1 √ N CP NCP-1 d=0 h UL i,u (dT s )e -j 2π N CP ku,n = w T i H AdDd u f ku,n , (52)
where f ku,n is the k u,n -th column vector of the N CP × N CP DFT matrix F NCP , and

H AdDd u = [h 1 , ..., h NCP ] ∈ C BxBy×NCP (53) with h t = h UL,LoS u p (t -1) T s -τ UL,LoS u + L l=1 h UL,NLoS u,l
p (t -1) T s -τ UL u,l , t = 1, .., N CP , being the effective angular-domain and delay-domain (AdDd) uplink channel to be estimated.

We assume that each UE transmits the pilot signal P UL Tx N used over its DSCs during the uplink CE stage, where P UL Tx is the total transmit power of the UE. The received pilot signal in the i-th time slot at the BS can be expressed as

y i,u,n = P UL Tx N used h UL,Fd i,u,n +n i,u,n = P UL Tx N used w T i H AdDd u f ku,n +n i,u,n , (54)
where n i,u,n ∼ CN (0, σ 2 n ) is the AWGN. By collecting the received pilot signals of all the DSCs for the u-th UE {y i,u,n } N used n=1 , we have

y T i,u = P UL Tx N used w T i H AdDd u F u + n T i,u , (55)
where

y i,u = [y i,u,1 , ..., y i,u,N used ] T ∈ C N used ×1 , n i,u = [n i,u,1 , ..., n i,u,N used ] T ∈ C N used ×1
, and

F u = [F NCP ] Ku ∈
C NCP×N used is the partial DFT matrix. For N P successive time slots, we aggregate the channel observations of the u-th UE y T i,u NP i=1 into the matrix Y u ∈ C NP×N used , which can be formulated as follows

Y u = [y 1,u , ..., y NP,u ] T = WH AdDd u F u + N u , (56)
u is referred to as the dual sparsity in both the angular domain and delay domain, which can be exploited to solve the CE problem.

Based on these considerations, Algorithm 2 reports the details of the proposed uplink CE scheme, where the CS-based orthogonal matching pursuit (OMP) algorithm is adopted to recover sparse channels based on (56). Once the estimated channel ĤAdDd u is obtained, interpolation-based methods can be applied to reconstruct the channels in all K subcarriers [START_REF] Zheng | Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization[END_REF].

Remark: The proposed closed-loop CE framework can be extended to the case where the LoS paths between the BS and the UEs exist. To this end, the RIS can be configured in order to reduce the scattering from it (i.e., absorption state [START_REF] Huang | Holographic MIMO surfaces for 6G wireless networks: Opportunities, challenges, and trends[END_REF]). Under this configuration, the UEs do not receive the signal reflected by the RIS. Therefore, the direct channel between the BS and the UEs can be estimated by using state-of-the-art CE algorithms (see, e.g., [START_REF] Lee | Channel estimation via orthogonal matching pursuit for hybrid MIMO systems in millimeter wave communications[END_REF], [START_REF] Gao | Spatially common sparsity based adaptive channel estimation and feedback for FDD massive MIMO[END_REF]- [START_REF] Liao | Closed-loop sparse channel estimation for wideband millimeter-wave full-dimensional MIMO systems[END_REF]). It is worth mentioning that the proposed closed-loop CE scheme can be applied to conventional MIMO (in the absence of RISs) communications in order to estimate the direct channels between the BS and the UEs. Once the direct link is estimated, it can be removed from the signal received in the presence of the RIS, and the proposed CE algorithm can be applied to estimate the channel between the RIS and the UEs.

C. Pilot Overhead and Computational Complexity Analysis

In this subsection, we analyze pilot overhead and the computational complexity of the proposed closed-loop CE scheme. As far as the pilot overhead is concerned, we evince from Algorithm 1 and Algorithm 2 that

G x G y M U x M U y Algorithm 2
The uplink CE scheme 1: Determine G ĝx,ĝy , and the DSC allocation scheme K u , ∀u ∈ {1, ..., N UE } for N UE UEs; 2: Set the beamforming design at each UE as in (42); 3: for the i-th time slot, 1 ≤ i ≤ N P , do 4:

Generate w i and accordingly set the beamforming design at the RIS as in (46); The BS receives the pilot signals y i,u from each UE in (55); 7: end for 8: Compute W, Y u , and F u in (56), ∀u ∈ {1, ..., N UE }; 9: % CS-based CE algorithm below. Take a certain u-th UE as an example. 10: Initialization: iter = 0, I = empty set, r = vec(Y u ), Ĥ is an all-zero matrix of size B x B y × N CP , N max is the maximum number of iterations; 11: while iter < N max , do and N P UWs are required for the downlink and uplink CE stages, respectively. Therefore, the total required pilot overhead of the proposed CE scheme is

T DL + T UL = 2N CP T s (G x G y M U x M U y + N P )
, where 2N CP T s is the length (duration) of one UW according to Fig. 7.

The computational complexity of the proposed scheme consists of the following two parts:

1) Downlink computational complexity. This is determined by the search at each UE in order to determine the angular groups that the LoS angles belong to, as detailed in [START_REF] Huemer | Classical and Bayesian linear data estimators for unique word OFDM[END_REF]. In this case, each UE needs to find the index of the maximum among G x G y M U

x M U y signals, and, therefore, the computational complexity at each UE is O(G x G y M U

x M U y ), where O(N ) stands for "of the order of N ". It is noteworthy that this complexity is affordable even for energy-constrained UEs, since M U

x , M U y , G x and G y are much smaller than the numbers of elements available at the RIS and BS.

2) Uplink computational complexity. This is mainly determined by the uplink CE algorithm. Since the CE problem has been formulated in (56) as a sparse signal recovery problem via under-determined measurements, various off-theshelf algorithms, such as greedy algorithms (e.g., [START_REF] Wan | Broadband channel estimation for intelligent reflecting surface aided mmWave massive MIMO systems[END_REF], [START_REF] Ma | Joint channel estimation and data rate maximization for intelligent reflecting surface assisted terahertz MIMO communication systems[END_REF]), Bayesian algorithms (e.g., [START_REF] He | Cascaded channel estimation for large intelligent metasurface assisted massive MIMO[END_REF]), and deep learning methods (e.g., [START_REF] Liu | Deep denoising neural network assisted compressive channel estimation for mmWave intelligent reflecting surfaces[END_REF]), can be used for CE. The corresponding computational complexity may, therefore, vary significantly. As far as the OMP algorithm adopted in Algorithm 2 is concerned, we consider the total number of complex-valued multiplications to evaluate the computational complexity, as listed in [START_REF] Wan | Compressive sensing based channel estimation for millimeter-wave full-dimensional MIMO with lens-array[END_REF]Table I]. The specific results are presented in the next section.

V. SIMULATION RESULTS

In this section, we present numerical results to evaluate the different types of RISs, and the performance of the proposed CE scheme.

A. Experimental Setting

We consider the system model as shown in Fig. 9. The BS and the RIS serve the active UEs distributed within a sector of radius R and central angle 120 • . The BS and the RIS are elevated to the height of h 1 and the UEs have the height of h 2 . We assume that the normal directions of the arrays at the BS and RIS point towards each other, which yields As for the downlink channel model in [START_REF] Sarieddeen | Terahertz-band ultra-massive spatial modulation MIMO[END_REF], we set L = 1, K f = 30 dB. The angle µ LoS is calculated according to the position of the UE, and the angles ν LoS , µ l and ν l are randomly generated. The delay offsets τ UL,LoS and τ DL l (τ UL,LoS < τ DL l ) follow a uniform distribution U(0, (N CP -1)T s ). The channel coefficients α DL and β DL,LoS (β DL l can be similarly modeled) are modeled as follows [START_REF] Tang | Wireless communications with reconfigurable intelligent surface: Path loss modeling and experimental measurement[END_REF]Eq. (20)]

ψ B = ψ R = [0, 0] T .
α DL = e jθα G Tx S eff 4πR 2 A abs (f c , R) , ( 57
) BS RIS UE R 1 h 2 h RIS UE
β DL,LoS = e jθ β G RIS G Rx A abs (f c , d RIS→UE ) λ 4πd RIS→UE , (58)
where θ α , θ β ∼ U[0, 2π) are the phase shifts introduced by the channels, S eff is the effective reflection area of the RIS, d RIS→UE is the distance between the UE and the RIS, G Tx , G RIS , and G Rx are the array gains of the BS, the RIS, and the UE, respectively 3 , and A abs (f c , d Tx→Rx ) is the attenuation caused by molecular absorption [START_REF] Elayan | Terahertz band: The last piece of RF spectrum puzzle for communication systems[END_REF], [START_REF] Rappaport | Wireless communications and applications above 100 GHz: Opportunities and challenges for 6G and beyond[END_REF]. A abs (f c , d Tx→Rx ) is related to the carrier-frequency f c and the transmission distance d Tx→Rx , and its specific values are obtained based on the recommendations of the International Telecommunications Union (ITU) [START_REF]Attenuation by atmospheric gases[END_REF]. The effective reflection area S eff of the RIS can be modeled as its whole physical area (aperture) [START_REF] Tang | Wireless communications with reconfigurable intelligent surface: Path loss modeling and experimental measurement[END_REF] S eff = A x A y , for CMS, 47)-( 50) can be similarly modeled and thus are omitted for brevity.

B. Numerical Results

Fig. 10 shows the required pilot overhead and computational complexity versus the number of groups {G x , G y } when N P = 40. It is observed that, by setting different groups {G x , G y }, the proposed closed-loop scheme provides a trade-off between the pilot overhead and the computational complexity. Two cases deserve further attention in Fig. 10: (i) the case of G x = G y = 1 refers to estimating the complete channels in the uplink (in an open-loop manner). In this case, the uplink channel matrix to be estimated in (53) has a total size of 40, 000 × 64, which causes an unaffordable computational complexity and storage burden [START_REF] Liao | Closed-loop sparse channel estimation for wideband millimeter-wave full-dimensional MIMO systems[END_REF]; and (ii) the case of G x = 2A x /λ = 200 and G y = 2A y /λ = 200 refers to acquiring the complete CSI only in the downlink by exhaustive beam scanning. This case would suffer from a long T DL (about 0.6 sec as shown in Fig. 10), which may degrade the net spectral efficiency.

Fig. 11 shows the accuracy of the downlink CE stage by investigating the probability of grouping failure versus the total downlink transmission power at the BS. The indices obtained by Algorithm 1 are compared with the oracle LoS angles µ LoS and ν LoS to decide whether the downlink CE succeeds or not. Three different types of RISs, i.e., criticallyspaced RIS (d = λ/2), ultra-dense RIS (d < λ/2) and CMS (d → 0) are considered. We observe that the CMS provides the best performance among the three types of RISs. The performance of ultra-dense RISs is significantly better than traditional critically-spaced RISs, since ultra-dense RISs have larger effective reflection area and smaller sidelobes (i.e., higher array gain in the mainlobe) which bring a better receive signal-to-noise ratio (SNR). the performance of a practical ultra-dense RIS can well approach that of an ideal (but unrealistic) CMS. A practical d = λ/8 renders only a minor performance gap compared with the CMS, so it would be sufficient to treat the ultradense RIS as the holographic RIS. Moreover, it is observed that the downlink CE performance improves as the number of groups {G x , G y } increases. This is because with a larger {G x , G y } the SBF beam pattern has a narrower passband, which enhances the amplitude of the beam pattern (i.e., array gain) in the passband and yields a better receive SNR. Next, we investigate the performance of the uplink finergrained CE stage that is obtained by using Algorithm 2, and compare it with some existing schemes. We set G x = G y = 10, i.e., B x = B y = 20. In Figs. 12-14, the normalized mean square error (NMSE) is adopted as the performance metric of interest, which is given by E

ĤAdDd u -H AdDd u 2 F H AdDd u 2 F
. Fig. 12 depicts the uplink NMSE performance versus the uplink pilot overhead N P and the different DSC allocation schemes as illustrated in Fig. 8 for N UE = 4. It can be observed that the random DSC allocation scheme achieves better CE performance, while the uniform allocation scheme fails to work properly. If the random DSC allocation scheme is considered, in addition, a sufficiently high CE accuracy can be ensured even with a low compression ratio (N P N used )/(B x B y N CP ) in the range {0.0063, 0.0187, 0.0313, 0.0437, 0.0563, 0.0688} in Fig. 12. Due to its superior performance, the random DSC allocation strategy is adopted to obtain the rest of the results.

In Fig. 13, we plot the NMSE performance of the proposed CE scheme as a function of the number of simultaneouslyserved UEs N UE for N P = 40. As the benchmark, we adopt the simultaneous weighted OMP (SW-OMP) [START_REF] González-Coma | Channel estimation and hybrid precoding for frequency selective multiuser mmWave MIMO systems[END_REF] algorithm which estimates the channels in the frequency domain (rather than the delay domain) via well-determined measurements. When SW-OMP is considered, N P UWs are equally divided into N UE parts, each of which is dedicated for one UE. Therefore, only N P /N UE UWs are available for each UE to conduct CE. It can be observed that the proposed CE scheme outperforms the considered benchmark even if a larger number of UEs are served simultaneously. This is obtained because the proposed CE scheme exploits the dual sparsity of THz MIMO channels in both the angular domain and delay domain, while the CE scheme based on SW-OMP only utilizes the sparsity in the angular domain. Fig. 14 compares the NMSE performance of different CE schemes against the uplink transmission power P UL Tx for N UE = 4. As a benchmark, we consider the LS estimator [START_REF] Mishra | Channel estimation and low-complexity beamforming design for passive intelligent surface assisted MISO wireless energy transfer[END_REF], [START_REF] Zheng | Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization[END_REF] with well-determined measurements in both the angular domain and delay domain, which requires a large number N P = B x B y N UE = 1600 of UWs as pilot signals. By leveraging the dual sparsity of THz MIMO channels in both the angular domain and delay domain, the proposed CS-based CE scheme outperforms the LS estimator even if a much smaller number of pilot signals (N P = {40, 80}) is used. As a second benchmark scheme, we analyze the openloop CE scheme that implements the uplink CE stage without using the downlink grouping. In such a case, the UEs do not have any prior information of the coarsely-estimated LoS angles and thus the NBS beamforming in (42) is unavailable. Instead, random pilot signals [START_REF] Gao | Spatially common sparsity based adaptive channel estimation and feedback for FDD massive MIMO[END_REF], [START_REF] Gao | Channel estimation for millimeter-wave massive MIMO with hybrid precoding over frequency-selective fading channels[END_REF] are employed at the UEs to realize an omni-directional beam pattern. Compared with the NBS beamforming towards the coarsely-estimated LoS direction, the omni-directional beam pattern disperses the transmit energy towards several directions. This results in the poor CE performance, as demonstrated in Fig. 14.

In Fig. 15, we compare the average spectral efficiency (ASE) performance that is achieved by different schemes in order to further evaluate the accuracy of the CE schemes. The ASE of the u-th UE is defined as

ASE = 1 - T DL + T UL T coh × E 1 K K k=1 log 2 1 + P UL Tx h Fd u,k /(Kσ 2 n ) (60)
where h Fd u,k is the frequency-domain baseband chanof the k-th subcarrier by the NBS beamforming towards the estimated LoS direction of the RIS. The NBS beamforming towards the LoS direction is illustrated as an upper-bound (i.e., perfect CSI without pilot overhead). We observe that the ASE obtained by the proposed CE scheme has good tightness with the upper-bound, while the ASE obtained by the well-determined LS estimator is worse due to the long time duration that is required for pilot transmission. It is also observed that the ASE improves as the element spacing d decreases, and it can approach the performance of the ideal CMS with practical d (e.g., d = λ/8), which further verify that the ultra-dense RIS is a good realization of the holographic RIS. Further, the ASE that is obtained in the absence of an RIS is illustrated. In this case, we assume that the UE communicates with the BS via the NLoS link with perfect CSI and no pilot overhead. We observe that the virtual LoS link provided by an RIS in THz massive MIMO systems significantly increases the ASE.

In addition, we investigate the impact of quantization error, i.e., only a finite number of phase shifts can be realized in practice, on the performance of the proposed schemes, similar to [START_REF] Qian | Beamforming through reconfigurable intelligent surfaces in single-user MIMO systems: SNR distribution and scaling laws in the presence of channel fading and phase noise[END_REF], [START_REF] Ma | Joint channel estimation and data rate maximization for intelligent reflecting surface assisted terahertz MIMO communication systems[END_REF]. Note that some state-of-the-art architectures of reflecting element, such as semiconductor diodes in [START_REF] Abeywickrama | Intelligent reflecting surface: Practical phase shift model and beamforming optimization[END_REF], are not suited to much higher frequencies (e.g., THz) [START_REF] Carrasco | Tunable graphene reflective cells for THz reflectarrays and generalized law of reflection[END_REF]. Therefore, the analysis of the impact of the interplay between the phase and the amplitude of the reflecting elements of the THz holographic RIS is postponed to a future search contribution. We assume that the phase of each reflecting element of the RIS is quantized by using B quantization bits. This implies that the phases of the reflection coefficients are drawn from the finite set B = 2π 1-2 B-1

2 B , 2π 2-2 B-1 2 B , ..., 2π 2 B-1 2 B
. As far as the proposed beamforming designs (NBS and SBF) are concerned, this correspond to quantizing the reflection coefficients Φ(m, n) in [START_REF] Zhou | Robust beamforming design for intelligent reflecting surface aided MISO communication systems[END_REF] or [START_REF] Abeywickrama | Intelligent reflecting surface: Practical phase shift model and beamforming optimization[END_REF] as Φ Q (m, n) = |Φ(m, n)| e jφm,n , where φ m,n ∈ B is the phase element in B that is closest to the phase of Φ(m, n). Then, the corresponding beam patterns are computed via [START_REF] Basar | Wireless communications through reconfigurable intelligent surfaces[END_REF] by replacing Φ(m, n) with Φ Q (m, n). In Fig. 16, we illustrate some examples of the SBF beam patterns as a function of the number of quantization bits. If B = 1, the quantization error is too large and the accuracy of the beam patterns degrades significantly. If B = {2, 3}, on the other hand, the obtained beam patterns after quantization exhibit the desired band-pass properties, with an acceptable degradation compared with the ideal case without quantization error (B = ∞). In Fig. 17, we evaluate the impact of quantization error on the ASE. The numerical results confirm that the ASE performance degradation is not significant if B = {2, 3}. In addition, this holds true for any values d ≤ λ/2.

VI. CONCLUSIONS

Motivated by the concept of holographic communications, we studied the physical layer transmission of holographic RISs in the THz band, where a large number of sub-wavelength reconfigurable elements are densely integrated into a compact space to approach a spatially continuous aperture. We derived the beam pattern of a holographic RIS and proposed two beamforming designs that are formulated in closed-form expressions. Based on the proposed beamforming designs, we proposed a closed-loop broadband CE scheme for the RISaided THz massive MIMO systems. The proposed CE scheme encompasses a grouping stage and an uplink finergrained CE stage. In order to reduce the pilot overhead, a CS-based CE algorithm was proposed that exploits the dual sparsity of THz MIMO channels in both the angular domain and delay domain. Simulation results showed that holographic RISs are able to outperform traditional designs based on nonholographic RISs as well as communication schemes that do not use RISs. Possible future research directions based on the results obtained in this paper include the practical amplitude and phase shift model of reflecting elements, the analysis of hardware impairments, the development of robust signal processing methods, analysis of near-field and proof-of-concept validations and field experiments.

1 √

 1 ) LoS angle of the RIS-UE channel seen by the RIS (UE) µ l (ν l ) The l-th NLoS angle of the RIS-UE channel seen by the RIS (UE) τ DL,LoS (τ DL l ) Delay offsets of the LoS (l-th NLoS) path in the downlink effective baseband channel α DL , β DL,LoS , β DL l Channel coefficients of the BS-RIS-UE channel accounting for large-scale fading g B r (ψout) / g U (ψout) General beam pattern of the MIMO systems at the r-th RFC of the BS / at the UE g (ψout, ψ in ) General beam pattern of the RIS corresponding to an incident angle ψ in *: The uplink version of these symbols is obtained by replacing the superscript "DL" with "UL". The user index u may be added as a subscript if necessary. unit. C and Z are the sets of complex-valued numbers and integers, respectively. a ∝ b denotes a = Cb with C being a non-zero constant. [•] m and [•] m,n represent the m-th element of a vector and the m-th row, n-th column element of a matrix, respectively. [A] I denotes the submatrix consisting of the columns of A indexed by the ordered set I. diag(•), • F , and ⊗ represent the diagonalization, Frobenius norm, and Kronecker product, respectively. vec (•) is the vectorization operation according to the columns of the matrix, and vec -1 (•) is the corresponding inverse operation. card(I) is the cardinality of the set I. Ξ N (x) is the N -order Dirichlet kernel function given by Ξ N (x) = sin(N x/2) N sin(x/2) for x = 2kπ, and Ξ N (2kπ) = (-1) k(N -1) , k ∈ Z. The "sinc" function is defined as sinc (x) = sin x x for x = 0, and sinc (0) = 1. δ(x) is the Dirac function. F N is an N × N normalized discrete Fourier transformation (DFT) matrix whose elements are [F N ] m,n = N e -j 2π(m-1)(n-1) N . U(a, b) denotes the uniform distribution within (a, b). Important system parameters and channel variables used in this paper are collected in Table

Fig. 1 .

 1 Fig. 1. The model of an RIS-aided THz massive MIMO system. (a) A multiantenna UE is served by the BS with the help of an RIS when the LoS path is blocked by possible obstacles. (b) The hardware architectures at the BS and the UE.

Fig. 2 .

 2 Fig. 2. Illustration of different types of RIS. (a) Critically-spaced RIS with d = 0.5λ; (b) ultra-dense RIS with d < 0.5λ, which has a spatial quasi-continuous aperture; (c) CMS with d → 0 (ignoring the physical size of a single reflecting element), which has an ideal reconfigurable continuous aperture.

  accounting for the AoA of the incident signal1 . Similarly, for the AoD denoted by θ azi out and θ ele out , we define ψ out ∆

g 2 e

 2 (ψ out , ψ in ) = d 2π dnl-yn∆y) dkdl = A x A y 4π 2 e j d-Ax 2 ∆x e j d-Ay 2 ∆y 2π d

Fig. 3 .

 3 Fig. 3. An example of beam pattern of NBS based on (16). Nx = Ny = 64, d = λ/2. The desired angle is λ 2π ψopt = [0.6, -0.2] T .

d

  for k and equal to b, b + 2π d for l, we obtain the following beam pattern g SBF (ψ out , ψ in ; ψ min , ψ max ) ∝ e j d-Ax 2 ∆x e j d-Ay 2

Φ

  SBF (m, n; ψ min , ψ max ) Nx+1 2 dk e -j Ny +1 2 dl e j(dmk+dnl) dkdl = e -jd( mψ azi in +nψ ele in ) e jd mψ azi max -e jd mψ azi min d m e jdnψ ele max -e jdnψ ele min dn ,

Fig. 4 .

 4 Fig. 4. An example of beam pattern of SBF based on (23). Nx = Ny = 64, d = λ/2. The cut-off angles are λ 2π ψ min = [-0.2, 0.2] T and λ 2π ψmax = [0.2, 0.6] T .

Fig. 5 .

 5 Fig. 5. The beam patterns of NBS with normalized amplitudes when Ax = 16λ. (a) λ 2π ψ azi opt = 0.2; (b) λ 2π ψ azi opt = 0.9. Power leakage is observed in (b) for d = λ/2.

Fig. 6 .

 6 Fig. 6. The beam patterns of SBF with normalized amplitudes when Ax = 16λ. (a) λ 2π ψ azi min = -0.2, λ 2π ψ azi max = 0.2; (b) λ 2π ψ azi min = -1, λ 2π ψ azi max = -0.6. Power leakage is observed in (b) for d = λ/2.

Corollary 2 .

 2 Consider a CMS based on the SBF-based beamforming design. The reflection coefficients, ΦSBF (x, y; ψ opt ), can be formulated as ΦSBF (x, y; ψ min , ψ max ) ∝ e -j(xψ azi in +ȳψ ele in ) e j xψ azi max -e j xψ azi min x e j ȳψ azi max -e j ȳψ azi min ȳ ,

Fig. 7 .

 7 Fig. 7. The transmission frame structure of the system. T DL (T UL ): Pilot overhead of downlink (uplink) CE stage. T coh : Channel coherence time, during which the CSI is assumed to be invariant.

λ 2 Gx 2 Gy -λ Ay along the elevation Algorithm 1

 221 λ Ax along the azimuth direction, and to ψ ele max,ĝy -ψ ele min,ĝy = 2π λ The downlink CE scheme 1: Determine ψ B , ψ R , G x and G y .

= 1 √

 1 which consists of the LoS part and the NLoS part. h UL,LoS u,i (τ ) can be expressed as (48), as shown at the top of the next page, where w i BxBy e jθi,1,1 , ..., e jθ i,1,By , e jθi,2,1 , ..., e jθ i,Bx ,By T ∈ C BxBy×1 , and h UL,LoS u , as shown in (49) at the top of the next page, is the effective angular-domain LoS channel with finer angular resolution. The elements in h UL,LoS u show the relation between the LoS angle µ LoS u and all the pre-defined codewords ζ bx,by .

FrequencyFig. 8 .

 8 Fig. 8. Examples of the received UW at the BS with DSCs, where N CP = 16 and N UE = 4. Three different DSC allocation schemes are illustrated. (a) Block allocation; (b) Uniform allocation; (c) Random allocation.

  , and K u ∩ K u = ∅ for ∀u = u . During the uplink CE stage, the u-th UE transmits its own pilot signals by using only the N used subcarriers indexed by K u out of the N CP available subcarriers. On the other hand, no signals are transmitted in other subcarriers. This makes easier to separate the pilot signals of different UEs at the BSs, since different DSCs are used. Three examples of the possible structure of an UW with DSCs are illustrated in Fig. 8.

5 :

 5 All N UE UEs transmit the UWs with pilot signals by using the assigned DSCs; 6:

r

  = vec(Y u ) -F T u ⊗ W I ĥtemp ; 16: iter = iter + 1; 17: end while 18: ĥ = vec( Ĥ); 19: [ ĥ] I = ĥtemp ; 20: ĤAdDd u = vec -1 ( ĥ); Output: The estimated channels ĤAdDd u , ∀u ∈ {1, ..., N UE };

  The simulation setup is detailed as follows unless stated otherwise: R = 20 m,h 1 = 10 m, h 2 = 1.5 m, M B x = M B y = 64, N RF = 4, M U x = M U y = 8, A x = A y = 0.2 m, f c = 0.15 THz (λ ≈ 2 mm), T s = 2 × 10 -9 sec (bandwidth BW = 1/T s = 500 MHz), N CP = 64, K = 256.A raised cosine filter p(τ ) with roll-off factor 0.8 is employed. The noise power spectrum density at the receiver is σ 2 NSD = -174 dBm/Hz. Thus, the power of the AWGN σ 2 n is σ 2 n = σ 2 NSD × BW ≈ -87 dBm. The number of iterations N max in Algorithm 2 is set to 20.

Fig. 9 .

 9 Fig. 9. RIS-aided THz massive MIMO systems. (a) side view; (b) top view.

AxAy d 2 S

 2 ele , for DPA-based RIS with spacing d, (59) where S ele ≤ d 2 is the physical size of a single reflection element in a DPA-based RIS. We assume S ele ≈ 200 µm × 200 µm as demonstrated in [32, Fig. 1]. The parameters of the uplink channel in (

Fig. 11

 11 Fig.10shows the required pilot overhead and computational complexity versus the number of groups {G x , G y } when N P = 40. It is observed that, by setting different groups {G x , G y }, the proposed closed-loop scheme provides a trade-off between the pilot overhead and the computational complexity. Two cases deserve further attention in Fig.10:(i) the case of G x = G y = 1 refers to estimating the complete channels in the uplink (in an open-loop manner). In this case, the uplink channel matrix to be estimated in (53) has a total size of 40, 000 × 64, which causes an unaffordable computational complexity and storage burden[START_REF] Liao | Closed-loop sparse channel estimation for wideband millimeter-wave full-dimensional MIMO systems[END_REF]; and (ii) the case of G x = 2A x /λ = 200 and G y = 2A y /λ = 200 refers to acquiring the complete CSI only in the downlink by exhaustive beam scanning. This case would suffer from a long T DL (about 0.6 sec as shown in Fig.10), which may degrade the net spectral efficiency.Fig.11shows the accuracy of the downlink CE stage by investigating the probability of grouping failure versus the total downlink transmission power at the BS. The indices obtained by Algorithm 1 are compared with the oracle LoS angles µ LoS and ν LoS to decide whether the downlink CE succeeds or not. Three different types of RISs, i.e., criticallyspaced RIS (d = λ/2), ultra-dense RIS (d < λ/2) and CMS (d → 0) are considered. We observe that the CMS provides the best performance among the three types of RISs. The performance of ultra-dense RISs is significantly better than traditional critically-spaced RISs, since ultra-dense RISs have larger effective reflection area and smaller sidelobes (i.e., higher array gain in the mainlobe) which bring a better receive signal-to-noise ratio (SNR). Fig.11also reveals that3 The array gain can be calculated as G = 4π

Fig. 10 .Fig. 11 .

 1011 Fig.10. The trade-off between the total pilot overhead and computational complexity of the proposed closed-loop CE scheme, where N P = 40 is considered.

Fig. 12 .

 12 Fig. 12. NMSE performance versus both the uplink pilot overhead N P and DSC allocation schemes when N UE = 4 and P UL Tx = 23 dBm.

Fig. 13 .

 13 Fig.[START_REF] Li | Reconfigurable intelligent surface assisted UAV communication: Joint trajectory design and passive beamforming[END_REF]. NMSE performace veruse the number of UEs N UE when N P = 40 and P UL Tx = 23 dBm. The SW-OMP algorithm[START_REF] González-Coma | Channel estimation and hybrid precoding for frequency selective multiuser mmWave MIMO systems[END_REF] is adopted as the benchmark.

Fig. 14 .

 14 Fig. 14. NMSE performance of different CE schemes against the uplink transmission power P UL Tx when N UE = 4. The CMS is considered.

Fig. 15 .

 15 Fig. 15. Comparison of ASE performance when N UE = 4 and T coh = 5 ms. N P = 40 is considered for the proposed CS-based CE scheme.

Fig. 16 .

 16 Fig. 16. SBF beam patterns used in the proposed CE scheme with different quantization bits. The ultra-dense RIS with elements spacing d = λ/4 is considered.

Fig. 17 .

 17 Fig. 17. ASE performance with different quantization bits under the same scenario as Fig. 15. Perfect CSI is considered.

1 ΞΞ 2 (Ξ

 12 The result in[START_REF] Zheng | Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization[END_REF] can be obtained by letting dm = x and dn = y in[START_REF] Zhou | Robust beamforming design for intelligent reflecting surface aided MISO communication systems[END_REF]. As far as[START_REF] He | Cascaded channel estimation for large intelligent metasurface assisted massive MIMO[END_REF] is concerned, in particular, we have gNBS (ψ out , ψ in ; ψ opt )∝ lim d→0 g NBS (ψ out , ψ in ; ψ opt ) = lim d→0 e j d-Ax2 ∆x e j d-Ay2 ∆y Ξ Nx [d (k opt -∆ x )] Ξ Ny [d (l opt -∆ y )] = e -j Ax 2 ∆x e -j Ay 2 ∆y lim d→0 Ξ Nx [d (k opt -∆ x )] Ξ Ny [d (l opt -∆ y )] . (61)Given the symmetry, we only need to prove limd→0 Nx [d (k opt -∆ x )] = sinc Ax 2 (k opt -∆ x )in order to obtain[START_REF] He | Cascaded channel estimation for large intelligent metasurface assisted massive MIMO[END_REF]. When k opt = ∆ x , based on the definition of the function Ξ Nx , we havelim d→0 Nx [d (k opt -∆ x )] = lim d→0 sin Nxd 2 (k opt -∆ x ) N x sin d 2 (k opt -∆ x ) Ax 2 (k opt -∆ x ) A x sin d 2 (k opt -∆ x ) k opt -∆ x ) Ax 2 (k opt -∆ x ) cos d 2 (k opt -∆ x ) = sin Ax 2 (k opt -∆ x ) Ax 2 (k opt -∆ x ) = sinc A x 2 (k opt -∆ x ) ,(62)where the equality (a) follows from N x d = A x , and the equality (b) is obtained by applying the De l'Hôpital rule with respect to d. In addition, when k opt = ∆ x , lim d→0 Nx (d • 0) = 1 = sinc Ax 2 • 0 . Thus, we obtain (62). This completes the proof of Corollary 1.

  Number of NLoS paths of the RIS-UE channel, and the corresponding Rican factor ψ B (ψ R ) LoS angle of the BS-RIS channel seen by the BS (RIS) µ LoS (ν LoS

			TABLE I
		SOME IMPORTANT SYSTEM PARAMETERS AND CHANNEL VARIABLES
	Notation	Definition	Notation	Definition
	M B	Number of antennas at the BS	M U	Number of antennas at the UE
	N RF	Number of RFCs at the BS	Ax × Ay	Physical size of the RIS
	d	Element spacing of the DPA-based RIS Φ(m, n)	Reflection coefficient of the (m,n)-th RIS element
		Notation *		Definition
	L, K f

  gx,gy,nx,ny LoS ; ψ opt,nx,ny , G U l,nx,ny = g U NBS ν l ; ψ opt,nx,ny , GLoS LoS , ψ R ; ψ min,gx,gy , ψ max,gx,gy , and Gl,gx,gy = gSBF µ l , ψ R ; ψ min,gx,gy , ψ max,gx,gy . In particular, gSBF and g U SBF are the beam patterns corresponding to ΦSBF in[START_REF] Hu | Beyond massive MIMO: The potential of data transmission with large intelligent surfaces[END_REF] and Φ U NBS in[START_REF] Huang | Holographic MIMO surfaces for 6G wireless networks: Opportunities, challenges, and trends[END_REF], respectively. After collecting y k,gx,gy,nx,ny for 1

	=		P DL Tx N CP	γ DL,LoS k	GLoS gx,gy G U,LoS nx,ny
	+	√	1 LK f	L l=1	γ DL l,k Gl,gx,gy G U l,nx,ny	+ n k,gx,gy,nx,ny , (39)
	where P DL Tx	is the total downlink transmit power,
	n k,gx,gy,nx,ny is the additive white Gaussian noise
	(AWGN) with distribution CN (0, σ 2 n ), G U,LoS nx,ny	=
	g U NBS ν gx,gy	= gSBF µ

  ,B y ∈ Z without loss of generality). In particular, the target of the uplink finer-grained CE is to solve the following optimization problem ( bx , by ) = arg max

			gx -ψ azi min,gx λ/Ax	+ 1 = 2Ax λGx , B y	∆ =
	λ 2π	ψ ele max,gy -ψ ele min,gy λ/Ay	+ 1 = λGy (B (bx,by) 2Ay gNBS µ LoS u , ψ R ; ζ bx,by ,	(45)
	where ζ bx,by = ζ azi bx , ζ ele by	T

x

The proposed beamforming design can be applied to the active beamforming for MIMO systems as well. In this case, we can set ψ in = [0, 0] T , and Φ(m, n) can be considered as the phase shift induced by the phase shifter.

where

The objective of the uplink CE stage is therefore, to estimate H AdDd u by exploiting the knowledge of W, F u , and the noisy matrix Y u . Usually, we have N P < B x B y due to the limited channel coherence time, and N used < N CP due to the allocation of different DSCs among multiple UEs. These constraints make (56) an underdetermined system, which is not possible to solve by using traditional estimation techniques such as the LS estimator [START_REF] Mishra | Channel estimation and low-complexity beamforming design for passive intelligent surface assisted MISO wireless energy transfer[END_REF], [START_REF] Zheng | Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization[END_REF]. Fortunately, due to the strong LoS link between the UE and the RIS, the channel matrix H AdDd u is expected to be sparse, and, in particular, only the elements whose indices fulfill the conditions µ LoS