# Time-Domain Channel Estimation for
Extremely Large MIMO THz Communications
with Beam Squint

Evangelos Vlachos , Aryan Kaushik ,
Yonina C. Eldar , and George C. Alexandropoulos E. Vlachos is with Industrial
Systems Institute, ATHENA Research and Innovation Centre, 26504 Rio-Patras,
Greece (e-mail: evlachos@athenarc.gr).A. Kaushik is with School of Engineering
and Informatics, University of Sussex, UK (e-mail:
aryan.kaushik@sussex.ac.uk).Y. C. Eldar is in the Faculty of Mathematics and
Computer Science, Weizmann institute of Science, Israel
(e-mail:yonina.eldar@weizmann.ac.il).G. C. Alexandropoulos is with the
Department of Informatics and Telecommunications, National and Kapodistrian
University of Athens, 15784 Athens, Greece (e-mail: alexandg@di.uoa.gr).

###### Abstract

In this paper, we study the problem of extremely large (XL) multiple-input
multiple-output (MIMO) channel estimation in the Terahertz (THz) frequency
band, considering the presence of propagation delays across the entire array
apertures, which leads to frequency selectivity, a problem known as beam
squint. Multi-carrier transmission schemes which are usually deployed to
address this problem, suffer from high peak-to-average power ratio, which is
specifically dominant in THz communications where low transmit power is
realized. Diverging from the usual approach, we devise a novel channel
estimation problem formulation in the time domain for single-carrier (SC)
modulation, which favors transmissions in THz, and incorporate the beam-squint
effect in a sparse vector recovery problem that is solved via sparse
optimization tools. In particular, the beam squint and the sparse MIMO channel
are jointly tracked by using an alternating minimization approach that
decomposes the two estimation problems. The presented performance evaluation
results validate that the proposed SC technique exhibits superior performance
than the conventional one as well as than state-of-the-art multi-carrier
approaches.

###### Index Terms:

Channel estimation, beam squint, extremely large MIMO, THz, alternating
minimization, single-carrier transmission, sparse estimation.

##  I Introduction

Terahertz (THz) communications (in the range of $0.1-10$ THz) have recently
received noticeable attention within the global wireless community due to
their increased potential for: seemless data transfer, wide bandwidth that can
theoretically reach up to hundreds of Gigahertz (GHz), data rates of the order
of terabits per second enabling ultra–fast downloading for immersive
applications, and latency in the order of microseconds ($\mu$sec) [1].
Therefore, it has been recently recognised as one of the promising candidate
technology for future sixth generation (6G) wireless networks [2].

To confront with the high penetration loss at the THz frequency band,
extremely large (XL) multiple-input multiple-output (MIMO) are being
considered [3], capable of realizing highly directive beamforming. However,
due to the ultra-high bandwidth employed in THz communications, the
propagation delay across the large antenna arrays at the communication
terminals can exceed the sampling period. In addition, as the antenna size
increases, its near-field range gets expanded, and thus, the phase of the
array steering vector becomes nonlinear to the antenna index. This spatial-
wideband effect causes the so-called beam squint in the frequency domain [4],
according to which the direction-of-arrival (DoA) varies with frequency, and
the array gain becomes frequency selective.

The predominant literature in channel estimation, with a specific focus on the
beam-squint effect, revolves around Orthogonal Frequency Division Multiplexing
(OFDM) schemes. OFDM, however, often grapples with the challenge of a high
peak-to-average power ratio (PAPR), a predicament exacerbated in the context
of ultra-high-frequency transmissions in the Terahertz (THz) range due to low
transmit power levels. Furthermore, THz-specific channel-induced impairments
and the presence of phase noise have been documented, posing additional
hurdles for multi-carrier (MC) transmission strategies. On the other side, the
single-carrier (SC) waveform is known for having a lower PAPR compared to the
OFDM waveform and being robust to system impairments and phase noise.
Especially for THz communications, due to the low output power and the non-
linearity effect induced by the available THz power amplifiers (PAs), it is
preferable to use SC transmissions rather than OFDM [5].

###  I-A Literature Review

In massive MIMO systems, accurate channel state information plays a vital role
to achieve capacity-approaching performance. To leverage the sparsity of the
channel in both the angle and delay domains, several channel estimation
techniques have been devised specifically for millimeter-wave (mmWave) [6, 7]
and THz communications [8]. In [8], the angle-delay sparsity of mmWave
transmissions have been exploited to design a compressive sensing based
channel estimation algorithm. In [7], a framework for exploiting both the
sparsity and low-rank properties was proposed and the beam-squint effect was
discussed. A range-azimuth Doppler recovery algorithm exploiting channel
sparsity, as well as compressive-sensing-based sparse recovery for MIMO radar
applications were presented in [9] and [10], respectively. The authors in [11]
presented a mmWave MIMO channel estimation approach relying on a non-linear
model and sparse recovery, which was based on expectation-maximization
generalized approximate message passing.

The utilization of compressive sensing methods to recover sparse channels in
THz channel estimation formulations [12, 12, 13, 14] is motivated by their
effectiveness in mmWave communications, where the channels exhibit similarly
lower sparsity levels. Despite the sparsity of the channel, the real-time
estimation of THz channels can be computationally intensive in dense multi-
user wideband scenarios with numerous channel paths. Compressive-sensing-based
estimation may require substantial numbers of measurements leading to
extensive latency and computational complexity. In such cases, traditional
methods like minimum mean square error (MMSE) and least square (LS) channel
estimation can be employed to estimate the second-order statistics of THz
channels [15]. Additionally, joint activity detection and channel estimation
are effective techniques to reduce the reliance on pilots and the
computational complexity in wideband random massive access THz systems [16].
Fast channel tracking presents an alternative approach to minimize the
overhead of channel estimation in high mobility scenarios, as demonstrated in
[17] for THz beamspace massive MIMO. Ultimately, efficient hardware
implementations are necessary to validate the practicality of THz channel
estimation algorithms, as exemplified in [18] for mmWave beamspace channel
estimation algorithms that exploit sparsity [19].

Channel estimation along with the beam-squint effect have been recently
investigated in [4] for frequency-division duplex (FDD) mmWave massive MIMO
OFDM systems with hybrid analog/digital precoding. A technique based on super-
resolution compressive sensing was designed to extract the frequency-dependent
parameters (angle-of-arrivals and path delays) as well as the frequency-
independent complex gains. The case of reconfigurable intelligent surfaces
(RIS) aided communication systems [20] has been investigated for cascaded
channel estimation in wideband scenarios considering the effect of beam squint
[21]. Specifically, it has been observed that the mutual correlation function
between the spatial steering vectors and the two-hop cascaded channel
reflected by the RIS has two peaks, which leads to a pair of estimated angles
for a single propagation path, due to the effect of beam squint.

To extend the communication distance in THz wireless systems, several advanced
solutions have been proposed, e.g., adaptive designs at the physical layer
[22, 23, 24] and time reversal [25], ultra-massive multiple-input multiple-
output (MIMO) [26, 27, 28], and multi-functional RISs [29, 30, 31, 32, 33, 34,
35]. At this frequency band, due to the small form factor and inter-element
spacing, it is feasible to design XL antenna arrays with very large numbers of
antenna elements, which enables highly directive beamforming that can combat
the high propagation loss. Typically, in Uniform Planar Arrays (UPAs),
hundreds or even thousands of densely packed antennas are being considered [1,
28]. These systems are of special interest as they can effectively increase
communication range and further enhance capacity in THz wireless
communications. On another direction, RISs have recently emerged as a
promising new paradigm to achieve smart and reconfigurable wireless
propagation environments [36, 37, 38], and are lately being studied for THz
communications [39].

The beamforming gain enhancement offered by the aforementioned technologies
can effectively combat stronger atmospheric attenuation and higher free space
pathloss appearing at THz frequencies. In conventional narrowband systems,
hybrid analog and digital beamforming schemes can achieve the full array gain
by aligning the beams generated by the analog beamformer with the physical
directions of the channel path components [3]. Theoretically, the far-field
beam squint indicates that the beams received at different frequencies squint
toward different directions, while the near-field beam squint indicates that
the beams at different frequencies aim at different locations and that most of
beams fail to focus on a targeted user [40].

###  I-B Motivation and Contributions

To the best of our knowledge, the effect of beam squint has not yet been
studied for SC schemes. Although the problems of beam squint and channel
estimation have been investigated for MC systems in THz communications, this
is the first work that investigates the joint problem of beam-squint and
channel estimation for SC transmissions. As previously mentioned, the low
output power and the non-linearity induced by available THz PAs motivates the
adoption of SC transmissions rather than OFDM [5]. Moreover, the first sub-THz
standard (IEEE 802.15.3d [41]) describes a SC modulation mode to support long
range and high rate wireless applications (such as $100$ Gbps). On this
premise, this work studies the problem of SC XL MIMO channel estimation in the
time domain, and presents a novel joint beam-squint and channel estimation
algorithm.

The contributions of this paper are summarized as follows:

  * •

We introduce an innovative formulation for the beam-squint effect in the
context of SC transmissions for THz wireless communications. In particular, we
represent this effect as a unknown binary vector which is included in the
considered channel estimation problem.

  * •

We present a novel time-domain joint beam-squint and channel estimation
problem, which serves as a more generalized formulation of a mixed-integer
non-linear problem, beyond those incorporating typical quadratic programming
expressions. To tackle this problem, we develop a specialized alternating
minimization solver.

  * •

Given the inherent complexity of the considered non-linear estimation problem,
our proposed technique leverages the potential availability of position
information of a user equipment (UE) at the base station (BS) to enhance the
accuracy of the estimation process.

  * •

We explore the robustness of the proposed technique in the presence of
deviations from the true UE position, which result in erroneous Angle of
Arrival (AoA) derivations, thus, in erroneous partial composition of the line-
of-sight (LOS) component of the unknown XL MIMO matrix.

Extensive simulation results are conducted to evaluate the proposed estimation
approach and compare it with other SC and MC benchmark techniques that are
available in the open technical literature.

###  I-C Notation and Organization

A summary of the notation used throughout this manuscript is included in Table
I.

TABLE I: The Notations of this Paper. $a,\mathbf{a}$, and $\mathbf{A}$ | Scalar, vector, and matrix
---|---
$j\triangleq\sqrt{-1}$ | The imaginary unit
$\mathbf{A}^{\text{T}}$ and $\mathbf{A}^{\text{H}}$ | Matrix transpose and Hermitian transpose
$\mathbf{A}^{-1}$ | Matrix inverse and pseudo-inverse
$[\mathbf{A}]_{i,j}$ | Matrix $\mathbf{A}$ element at the $i$-th row
| and $j$-th column
$[\mathbf{a}]_{i}$ | The $i$-th element of vector $\mathbf{a}$
$\mathbf{I}_{N}$ |  $N\times N$ identity matrix
$\mathbf{1}_{N}$ |  $N\times 1$ vector containing only $1$’s
$\times$ | Scalar multiplication
$\otimes$ | Kronecker product
$\textrm{blkdiag}(\mathbf{A}_{1},\mathbf{A}_{2},\ldots)$ | Places the input matrices
| on the main diagonal
$\mathbb{R},\mathbb{C}$ | The sets of real and complex numbers
$x\sim\mathcal{CN}(0,\sigma^{2})$ |  $x$ is a zero-mean complex Gaussian random
| variable with variance $\sigma^{2}$
$\mathcal{E}\\{x\\}$ | Expectation of random variable $x$

The remainder of the paper is organized as follows: Section II presents the
considered channel and system models, while Section III includes the proposed
problem formulation and the position-aided channel estimation approach.
Section IV verifies the proposed estimation framework through simulation
results. Section V sums up the outcomes of the proposed framework and sketches
directions for future research.

##  II Considered Models

We consider the downlink of a MIMO THz communication system comprising a
multi-antenna BS and a multi-antenna UE. Both nodes are assumed equipped with
uniform linear arrays (ULAs) consisting of $N$ and $M$ antenna elements,
respectively.

###  II-A Channel Model

To model the THz massive MIMO wireless channel, we adopt the Saleh-Valenzuela
(S-V) channel model that is based on the time-cluster spatial-lobe approach
[42]. According to this approach, time clusters are composed of multipath
components traveling close in time that arrive from potentially different
directions in a short propagation time window. The spatial lobes represent the
main directions of arrival (or departure) where energy arrives (or departs)
over several hundreds of nanoseconds. The propagation channel between the BS
and the UE is represented by the channel response matrix
$\mathbf{H}\in\mathbb{C}^{M\times N}$ comprising $L$ channel paths, which is
mathematically modeled as follows:

$\mathbf{H}\triangleq\sum_{\ell=1}^{L}\mathbf{H}(\ell)=\sum_{\ell=1}^{L}h_{\ell}\mathbf{a}_{\text{R}}(\theta_{\ell})\mathbf{a}_{\text{T}}^{\text{H}}(\phi_{\ell}),$
(1)

where $h_{\ell}$ is the random complex channel gain, and
$\mathbf{a}_{\text{R}}(\theta_{\ell})\in\mathbb{C}^{M\times 1}$ is the $M$-th
dimensional array response vector at the receiving UE, expressed as:

$\mathbf{a}_{\text{R}}(\theta_{\ell})=\frac{1}{\sqrt{M}}\left[1,e^{j2\pi\theta_{\ell}},\ldots,e^{j2\pi(N-1)\theta_{\ell}}\right]^{\text{T}},$
(2)

where the normalized AoA is defined as
$\theta_{\ell}\triangleq\frac{d\sin\vartheta_{\ell}}{\lambda}$ with
$\lambda\triangleq c/f_{c}$ being the wavelength of the carrier, having the
frequency $f_{c}$, and $c$ denotes the speed of light. The array response
vector $\mathbf{a}_{\text{T}}(\phi_{\ell})\in\mathbb{C}^{N\times 1}$ at the
transmitting BS for the planar wave departing from the normalized angle of
departure (AoD) $\phi_{\ell}\triangleq\frac{d\sin\varphi_{\ell}}{\lambda}$ in
(1) is defined in a similar manner.

The variance of each $\ell$-th channel gain coefficient $h_{\ell}$ depends on
the respective propagation distance $d_{\text{BS-UE}}$ between the BS and the
UE as well as the carrier frequency $f_{c}$ via the following expression:

$\mathcal{E}\\{|h_{\ell}|^{2}\\}=P_{t}\sqrt{\frac{NM}{L}}\frac{1}{d_{\text{BS-
UE}}^{\xi_{\ell}}}e^{-\frac{1}{2}\mathcal{K}(f_{c})},$ (3)

where $P_{t}$ is the BS transmit power, $\xi_{\ell}$ represents the pathloss
exponent which is equal to $\xi_{1}=2$ for the LoS path, i.e., $\ell=1$, and
$\xi_{\ell}=3$ for $\ell=2,\ldots,L$. $\mathcal{K}(f)$ represents the function
of the molecular absorption losses that depends on the carrier frequency [43].

###  II-B The Beam-Squint Effect

We assume that the considered point-to-point MIMO communication takes place on
a frame-by-frame basis, where the wireless channel remains constant during
each frame, but it may change independently from one frame to another. Every
frame consists of $T$ time slots dedicated for channel estimation, whereas the
rest of the frame is used for data communication. To estimate the intended THz
MIMO channel, the $M$-antenna BS utilizes training symbols for each of the $T$
slots used for channel estimation.

The sampled signal at the $m$-th ($m=1,2,\ldots,M$) receive antenna of the UE
at each $t$-th discrete time slot can be mathematically expressed as follows
[4]:

$y_{m}(t)\triangleq\sum_{n=1}^{N}\sum_{\ell=1}^{L}[\mathbf{H}]_{m,n}(\ell)q_{n}(t-\tau_{\ell,m})+n_{m}(t),$
(4)

where $\tau_{\ell,m}$ is the time delay of the $\ell$-th channel path for the
$m$-th receiving antenna, $q_{n}$ denotes the training symbol transmitted from
each $n$-th ($n=1,2,\ldots,N$) transmit (TX) antenna element, and
$n_{m}(t)\sim\mathcal{N}(0,\sigma_{n}^{2})$ is the complex additive white
Gaussian noise (AWGN) having zero mean and variance $\sigma_{n}^{2}$.

Following the far-field assumption, according to which the antenna array sizes
are much smaller than the distance between the TX and the receiver (RX), it is
known that [4]:

$\tau_{\ell,m}=\tau_{\ell,1}+(m-1)\frac{d\sin(\vartheta_{\ell})}{c},$ (5)

where $\tau_{\ell,1}$ is the time delay at the first receive antenna element
due to the THz channel characteristics. Considering the case where the inter-
element antenna spacing is $d=\lambda/2$, yields $d/c=1/(2f_{c})$, resulting
in the following expression:

$\tau_{\ell,m}=\tau_{\ell,1}+\frac{(m-1)}{2f_{c}}\sin{\vartheta_{\ell}},$ (6)

while the difference of the time delay between the first and each $m$-th
antenna becomes:

$\tau_{\ell,m}-\tau_{\ell,1}=\frac{(m-1)}{2f_{c}}\sin(\vartheta_{\ell}).$ (7)

The right hand term in (7), which represents the time delay due to the beam-
squint effect, becomes non-negligible for higher sampling frequencies,
$f_{s}\triangleq 1/T_{s}$, as it can also be seen in Fig. 1. The following
proposistion presents the bounded worst-case time delay due to the beam-squint
effect.

Figure 1: Propagation time delay as a function of the carrier frequency
$f_{c}$ in GHz due to the beam-squint effect, considering different sizes $M$
of the receiving UE antennas and different values for the sampling frequency
$f_{s}$.

###### Proposition 1

The $T_{s}$-sampled training signal from the $n$-th TX antenna element, which
is represented by $q_{n}(t-\tau_{\ell,1}-\tau_{\ell,m})$, is shifted in time
due to the beam squint effect at most by the following quantity for all AoAs
$\vartheta_{\ell}$’s and all RX antenna indices $m$:

$|\tau_{\ell,m}-\tau_{\ell,1}|\leq\frac{M-1}{2f_{c}T_{s}}.$ (8)

The received training signal in (4) can be expressed in the discrete time
domain as $q_{n}(k-\kappa_{\ell,m})$, where $k$ represents the time slots and
$\kappa_{\ell,m}$ denotes the time shift that occurs at each $m$-th RX antenna
terminal due to each $\ell$-th signal propagation path. Then, the discrete
training signal can be written as an inner product, as follows:

$q_{n}(k-\kappa_{\ell,m})\equiv\mathbf{e}_{\ell,m}^{\text{T}}\mathbf{q}_{n}(k),$
(9)

where
$\mathbf{q}_{n}(k)\triangleq[q_{n}(k-\tau_{\text{max}}),\ldots,q_{n}(k+\tau_{\text{max}})]^{\text{T}}$
and $\mathbf{e}_{\ell,m}\in\\{0,1\\}^{M\times 1}$ is a binary vector with
zeros everywhere except the $\kappa_{\ell,m}$-th position.

By leveraging Proposition 1, we can establish that in the most challenging
scenario, when the product of $f_{c}$ and $T_{s}$ is approximately equal to 1,
i.e., $f_{c}T_{s}\approx 1$, the maximum allowable shift is given by
$\tau_{\text{max}}\triangleq\lceil\frac{M-1}{2}\rceil$, defined as the
smallest integer greater than or equal to $\frac{M-1}{2}$. To simplify our
analysis, we make the assumption that $\tau_{\text{max}}$ must be an integer,
leading to the relationship:

$2\tau_{\text{max}}+1=M.$

Depending on the severity of the beam-squint effect, one or more UE antennas
may be affected from it. We make use of $\mathbf{e}_{\ell,m}$ to model this
effect at each $m$-th receiving UE antenna element via its non-zero entries.
To this end, for each $k$-th time slot, $y_{m}(k)$ in (4) can be expressed by
the more general, beam-squint-including, expression:

$y_{m}(k)=\sum_{n=1}^{N}\sum_{\ell=1}^{L}\mathbf{H}_{m,n}(\ell)\mathbf{e}^{\text{T}}_{\ell,m}\mathbf{q}_{n}(k)+n_{m}(k).$
(10)

By introducing the $N\times M$ matrix $\mathbf{Q}(k)$, with its $n$-th row
including the vector $\mathbf{q}_{n}^{\text{T}}(k)$, and the notation
$\mathbf{\bar{h}}_{m}(\ell)\in\mathbb{C}^{1\times N}$ to represent each $m$-th
row of $\mathbf{H}(\ell)$, expression (10) can be written more compactly in
matrix form, as follows:

$y_{m}(k)=\sum_{\ell=1}^{L}\mathbf{\bar{h}}_{m}(\ell)\mathbf{Q}(k)\mathbf{e}_{\ell,m}+n_{m}(k).$
(11)

##  III Joint Beam-Squint and Channel Estimation

In this section, we first present the proposed optimization formulation for
the joint beam-squint and channel estimation problem. Then, we detail an
efficient solution for this problem that is based on an alternating
minimization approach. We also discuss a position-aided initialization for the
proposed iterative estimation approach as well as present bounds on its
estimation performance.

###  III-A Problem Formulation

In this work, we consider the case of the same beam-squint effect in all
channel propagation paths, i.e.:

$\kappa_{1,m}=\kappa_{2,m}=\ldots=\kappa_{L,m}\triangleq\kappa_{m}.$ (12)

This case represents the scenario where all paths originate from the same far-
field source, thus, AoA $\vartheta_{\ell}$ are equal for $\ell=1,\ldots,L$,
and the time delay for the different antenna terminals in (7) does not depend
on $\ell$. Therefore, (11) can be re-expressed, using the notation
$\mathbf{h}_{m}\triangleq\sum_{\ell=1}^{L}\mathbf{h}_{m}(\ell)$, as follows:

$\displaystyle y_{m}(k)$
$\displaystyle=\sum_{\ell=1}^{L}\mathbf{\bar{h}}_{m}(\ell)\mathbf{Q}(k)\mathbf{e}_{m}+n_{m}(k)$
(13) $\displaystyle=\mathbf{h}_{m}\mathbf{Q}(k)\mathbf{e}_{m}+n_{m}(k),$ (14)

which coincides with (11) when $\mathbf{e}_{m}$ has its only non-zero element
in the position $\kappa_{m}$.

Gathering the received training signals of $T$ time instances at each $m$-th
UE antenna element, and using (11) and the notation
$\mathbf{n}_{m}\triangleq\left[n_{m}(1)\,\cdots\,n_{m}(T)\right]^{\text{T}}\in\mathbb{C}^{1\times
T}$, the following expression is deduced :

$\displaystyle\left[\begin{array}[]{c}y_{m}(1)\\\ \vdots\\\
y_{m}(T)\end{array}\right]\\!=\\!\left[\begin{array}[]{ccc}\mathbf{h}_{m}&&\\\
&\ddots&\\\
&&\mathbf{h}_{m}\end{array}\right]\left[\begin{array}[]{ccc}\mathbf{Q}(1)&&\\\
&\ddots&\\\
&&\mathbf{Q}(T)\end{array}\right]\left[\begin{array}[]{c}\mathbf{e}_{m}\\\
\vdots\\\ \mathbf{e}_{m}\end{array}\right]\\!+\\!\mathbf{n}_{m}.$ (27)

Making use of the notation
$\mathbf{y}_{m}\triangleq\left[y_{m}(1)\,\cdots\,y_{m}(T)\right]^{\text{T}}\in\mathbb{C}^{1\times
T}$ and the all-ones vector $\mathbf{1}_{T}\in\\{1\\}^{T\times 1}$, the latter
expression can be compactly re-written as follows:

$\mathbf{y}_{m}=(\mathbf{I}_{T}\otimes\mathbf{h}_{m})\mathbf{\bar{Q}}(\mathbf{1}_{T}\otimes\mathbf{e}_{m})+\mathbf{n}_{m},$
(28)

where the $T\times MT$ block diagonal matrix $\mathbf{\bar{Q}}$ is obtained
as:

$\mathbf{\bar{Q}}\triangleq\text{blkdiag}(\mathbf{Q}(1),\ldots,\mathbf{Q}(T)).$
(29)

Capitalizing on (28) for each $m$-th UE antenna element, we formulate the
following problem for the joint beam-squint and channel estimation objective:

$\displaystyle\min_{\mathbf{h}_{m},\mathbf{e}_{m}}$
$\displaystyle\|\mathbf{y}_{m}-(\mathbf{I}_{T}\otimes\mathbf{h}_{m})\mathbf{\bar{Q}}(\mathbf{1}\otimes\mathbf{e}_{m})\|_{2}^{2}$
s.t. $\displaystyle\mathbf{e}_{m}\in\\{0,1\\}^{M\times 1},$ (30)

where the unknowns $\mathbf{h}_{m}$ and $\mathbf{e}_{m}$ need to be jointly
recovered. It is worth emphasizing that (III-A) represents a broader class of
problems compared to the traditional mixed-integer quadratic problems (MIQP)
involving vector-valued variables [44]. Consequently, solving it, is not as
straightforward as applying conventional solvers designed for MIQP [45].

###  III-B Channel Sparsity Exploitation

The expression of the channel to the beamspace domain essentially provides its
decomposition into the discrete Fourier transform (DFT) domain. Specifically,
if $\mathbf{F}_{\rm RX}\in\mathbb{C}^{M\times M}$ and $\mathbf{F}_{\rm
TX}\in\mathbb{C}^{N\times N}$ are the DFT matrices, then it holds:

$\mathbf{Z}\triangleq\mathbf{F}_{\rm RX}\mathbf{H}\mathbf{F}_{\rm
TX}^{\text{H}},$ (31)

where the sparse matrix $\mathbf{Z}$ is the expression of the channel matrix
in the beamspace. If the number of the non-LOS (NLOS) components is small,
$\mathbf{Z}$ has a small number of high amplitude values, implying a small
$\ell_{1}$ norm. Moreover, the beamspace expression of the $m$-th row of the
channel matrix $\mathbf{H}$ is given by the product $[\mathbf{F}_{\rm
RX}]_{m,m}\mathbf{h}_{m}\mathbf{F}_{\rm TX}^{\text{H}}$. Note that, since the
scalar term $[\mathbf{F}_{\rm RX}]_{m,m}$ does not affect the sparsity, it is
ommited in the following analysis.

To exploit the inevitable sparsity property of the channel in the beamspace
domain, we enrich the formulation in (III-A) as follows:

$\displaystyle\min_{\mathbf{h}_{m},\mathbf{e}_{m}}$
$\displaystyle\|\mathbf{h}_{m}\mathbf{F}\|_{1}+\|\mathbf{y}_{m}-(\mathbf{I}_{T}\otimes\mathbf{h}_{m})\mathbf{\bar{Q}}(\mathbf{1}\otimes\mathbf{e}_{m})\|_{2}^{2}$
s.t. $\displaystyle\mathbf{e}_{m}\in\\{0,1\\}^{M\times 1}.$ (32)

###  III-C Problem Decomposition

To establish performance bounds for the estimation performance of the proposed
technique, we decompose problem (III-A) into two ideally separate problems
that are solved independently [46]:

  * •

Minimization over $\mathbf{e}_{m}$ assuming that the channel for the $m$-th
receiving UE antenna is perfectly known ($\mathbf{h}_{m}^{*}$), i.e.:

$\displaystyle\mathbf{e}_{m}^{*}=\arg\min_{\mathbf{e}_{m}}$
$\displaystyle\|\mathbf{y}_{m}-(\mathbf{I}_{T}\otimes\mathbf{h}_{m})\mathbf{\bar{Q}}(\mathbf{1}\otimes\mathbf{e}_{m})\|_{2}^{2}$
s.t. $\displaystyle\mathbf{e}_{m}\in\\{0,1\\}^{M\times 1},$ (33)

where $\mathbf{h}_{m}$ is assumed to be perfectly known.

  * •

Minimization over $\mathbf{h}_{m}$ assuming that the beam-squint vector with
respect to the $m$-th UE antenna is perfectly known ($\mathbf{e}_{m}^{*}$),
i.e.:

$\displaystyle\mathbf{h}_{m}^{*}\\!\\!=\\!\arg\min_{\mathbf{h}_{m}}$
$\displaystyle\|\mathbf{h}_{m}\mathbf{F}\|_{1}\\!+\\!\|\mathbf{y}_{m}\\!\\!-\\!(\mathbf{I}_{T}\\!\otimes\\!\mathbf{h}_{m})\mathbf{\bar{Q}}(\mathbf{1}\\!\otimes\\!\mathbf{e}_{m}^{*})\|_{2}^{2},$
(34)

where the beam-squint vector $\mathbf{e}_{m}$ is perfectly known.

The optimization problems at hand can be tackled using a variety of available
techniques, each offering distinct trade-offs between computational complexity
and solution accuracy. Some indicative approaches include Soft-Thresholding,
the Least Angle Regression (LARS) algorithm [47], and Iterative Thresholding,
among others. In our research, we opted to utilize the CVX package [48] due to
its user-friendly interface and widespread adoption within the scientific
community.

###  III-D Proposed Solution

To solve the joint optimization problem in (III-B), we utilze the alternating
method of multipliers (ADMM) according to which the estimation is decomposed
into subproblems that are solved iteratively [49]. Following the ADMM
approach, an auxialiry variable $\mathbf{c}_{m}\in\mathbb{C}^{1\times N}$ is
introduced to decompose the cost function of (III-B), as follows:

$\displaystyle\min_{\mathbf{c}_{m},\mathbf{h}_{m},\mathbf{e}_{m}}$
$\displaystyle\|\mathbf{c}_{m}\mathbf{F}\|_{1}+\|\mathbf{y}_{m}-(\mathbf{I}_{T}\otimes\mathbf{h}_{m})\mathbf{\bar{Q}}(\mathbf{1}\otimes\mathbf{e}_{m})\|_{2}^{2}$
s.t. $\displaystyle\mathbf{e}_{m}\in\\{0,1\\}^{M\times 1}\text{ and
}\mathbf{c}_{m}=\mathbf{h}_{m}.$ (35)

Then, the problem (III-D) can be decomposed into the following steps that
performed at each ADMM algorithm iteration $i$:

  * •

Derive $\mathbf{e}^{(i+1)}$ by minimizing:

$\min_{\mathbf{e}_{m}}\|\mathbf{y}_{m}-(\mathbf{I}_{T}\otimes\mathbf{h}_{m}^{(i)})\mathbf{\bar{Q}}(\mathbf{1}\otimes\mathbf{e}_{m})\|_{2}^{2}\text{s.t.
}\mathbf{e}_{m}\in\\{0,1\\}^{M\times 1}$ (36)

  * •

Derive $\mathbf{c}_{m}^{(i+1)}$ by minimizing:

$\min_{\mathbf{c}_{m}}\|\mathbf{c}_{m}\mathbf{F}\|_{1}+(\mathbf{c}_{m}-\mathbf{h}_{m}^{(i)})\mathbf{z}^{(i)}+\frac{\rho}{2}\|\mathbf{c}_{m}-\mathbf{h}_{m}^{(i)}\|_{2}^{2}$
(37)

  * •

Derive $\mathbf{h}_{m}^{(i+1)}$ by minimizing:

$\displaystyle\min_{\mathbf{h}_{m}}\|\mathbf{y}_{m}-$
$\displaystyle(\mathbf{I}_{T}\otimes\mathbf{h}_{m})\mathbf{w}\|_{2}^{2}$
$\displaystyle+(\mathbf{c}_{m}^{(i)}-\mathbf{h}_{m})\mathbf{z}^{(i)}+\frac{\rho}{2}\|\mathbf{c}_{m}^{(i)}-\mathbf{h}_{m}\|_{2}^{2},$
(38)

where $\mathbf{w}\triangleq\mathbf{\bar{Q}}(\mathbf{1}\otimes\mathbf{e}_{m})$.

  * •

Update the dual variable:

$\mathbf{z}^{(i+1)}=\mathbf{z}^{(i)}+\rho(\mathbf{c}_{m}^{(i+1)}-\mathbf{h}_{m}^{(i+1)})$

In this iterative procedure, $\mathbf{h}_{m}^{(0)}$ represents the initial
value of the estimated channel for each $m$-th receiving UE antenna. The
auxiliary variable $\mathbf{c}_{m}^{(0)}$ and the dual variable
$\mathbf{z}^{(0)}$ are initialized to zero.

It can be easily concluded that the previous Step 1 includes an integer
programming problem. Such problems are well known to be NP-hard, thus, being
very difficult to solve in general. To obtain a solution using standard tools
(i.e., CVX [48]), this problem can be relaxed to a linear one using box
constraints. Afterwards, the solution can be rounded to the closest integer
via a simple thresholding mechanism. To this end, we focus on the following
relaxed form of the optimization problem in (36):

$\displaystyle\mathbf{\tilde{e}}_{m}^{(i)}=\arg\min_{\mathbf{e}_{m}}$
$\displaystyle\|\mathbf{y}_{m}-(\mathbf{I}_{T}\otimes\mathbf{h}_{m}^{(i-1)})\mathbf{\bar{Q}}(\mathbf{1}\otimes\mathbf{e}_{m})\|_{2}^{2}$
s.t. $\displaystyle[\mathbf{e}_{m}]_{i}\in(0,1)\,\forall i=1,\ldots,M.$ (39)

This problem can be approximated using thresholding on a LS-based estimation
of $\mathbf{\tilde{e}}_{m}^{(i)}$, i.e.,

$\mathbf{\tilde{e}}_{m}^{(i)}=\left(\mathbf{\Phi}_{m}^{(i-1)}\right)^{-1}\mathbf{y}_{m},$

where
$\mathbf{\Phi}_{m}^{(i-1)}=\left[(\mathbf{h}_{m}^{(i-1)}\mathbf{Q}(1))^{\text{T}}\,\cdots\,(\mathbf{h}_{m}^{(i-1)}\mathbf{Q}(T))^{\text{T}}\right]^{\text{T}}$.
Afterwards, the values are being thresholded using a predefined value
$\kappa$, i.e.,
$\mathbf{e}_{m}^{(i)}=\text{thres}_{\kappa}\left(\mathbf{\tilde{e}}_{m}^{(i)}\right)$,
where the thresholding function $\text{thres}_{\kappa}(\cdot)$ is applied to
each element of the vector as follows:

$[\mathbf{e}_{m}]_{i}\triangleq\left\\{\begin{array}[]{cc}1,&\frac{|[\mathbf{e}_{m}]_{i}|}{\|\mathbf{e}_{m}\|_{2}}>\kappa\\\
0,&\text{otherwise}\end{array}\right..$ (40)

Concerning Step 2, the optimization problem stated in (37), this can be
expressed as a standard LASSO problem, e.g.,

$\min_{\mathbf{c}_{m}}\|\mathbf{c}_{m}\mathbf{F}\|_{1}+\frac{1}{2}\|\mathbf{c}_{m}-\mathbf{v}\|_{2}^{2},$
(41)

where $\mathbf{v}\triangleq\mathbf{z}^{(i)}+\rho\mathbf{h}_{m}^{(i)}$. This
problem can be solved efficiently by using the soft-thresholding operator,
$\mathcal{S}_{\tau}$, where the solution is given by:

$\mathbf{\bar{c}}^{(i+1)}=\mathcal{S}_{\tau}(\mathbf{v}).$ (42)

The 3rd Step of the ADMM algorithm involves the solution of (• ‣ III-D), which
has a closed form solution, given by:

$\mathbf{h}_{m}^{(i+1)}=(\mathbf{y}_{m}\mathbf{w}^{\text{H}}+\rho\mathbf{c}^{(i+1)}+\mathbf{z}^{(i)})(\mathbf{w}\mathbf{w}^{\text{H}})^{\dagger}.$
(43)

The proposed algorithm is summarized in 1.

#### Complexity Analysis

Regarding the computational complexity of Algorithm 1, the primary workload is
concentrated within lines 4 and 8. Specifically, line 4 requires the inversion
of the $T\times M$ matrix $\mathbf{\Phi}_{m}^{i-1}$ for $I_{\text{max}}$
iterations, thus the complexity order is bounded by $\mathcal{O}(M^{3}T)$ for
$M>T$. Line 8 requires the invertion of a $T\times T$ matrix, with general
complexity order $\mathcal{O}(T^{3})$. However, the algorithm complexity can
be significantly reduced by using more efficient variants, such as gradient
descent, or Woodbury identity lemma for the inversion of the rank-1 matrix
$\mathbf{w}\mathbf{w}^{\text{H}}$.

Algorithm 1 Estimation of each $m$-th vectors $\mathbf{h}_{m}$ and
$\mathbf{e}_{m}$

0: $\mathbf{y}_{m}$, $\phi_{1}$, $\theta_{1}$, $\mathbf{Q}$, $\mathbf{F}$,
$\kappa$, and $I_{\text{max}}$

1: Perform the initialization
$\mathbf{h}_{m}^{(0)}\leftarrow\left[\mathbf{a}_{\text{R}}(\theta_{1})\right]_{m}\mathbf{a}_{\text{T}}^{\text{H}}(\phi_{1})$

2: for $i=1,2,\ldots,I_{\text{max}}$ do

3: Set
$\mathbf{\Phi}_{m}^{(i-1)}=\left[(\mathbf{h}_{m}^{(i-1)}\mathbf{Q}(1))^{\text{T}}\,\cdots\,(\mathbf{h}_{m}^{(i-1)}\mathbf{Q}(T))^{\text{T}}\right]^{\text{T}}$

4: Compute
$\mathbf{\tilde{e}}_{m}^{(i)}=\left(\mathbf{\Phi}_{m}^{(i-1)}\right)^{-1}\mathbf{y}_{m}$.

5: Perform
$\mathbf{e}_{m}^{(i)}=\text{thres}_{\kappa}(\mathbf{\tilde{e}}_{m}^{(i)})$
using (40)

6: Compute $\mathbf{v}=\mathbf{z}^{(i)}+\rho\mathbf{h}_{m}^{(i)}$

7: Soft-thresholding using (42)

8: Update $\mathbf{h}_{m}^{(i+1)}$ using (43)

9: Dual variable udpate
$\mathbf{z}^{(i+1)}=\mathbf{z}^{(i)}+\rho(\mathbf{c}_{m}^{(i+1)}-\mathbf{h}_{m}^{(i+1)})$

10: end for

###  III-E Position-Aided Initialization

It has been stated that location information can contribute in addressing
several of the key challenges in non-terrestrial and vehicular networks, such
as the increase in traffic and number of devices, robustness for mission
critical services, and the reduction in total energy consumption and latency
[50]. For instance, a point-to-point mmWave backhaul communication system with
large antenna arrays at both communication ends was considered in [51], where
position information available from displacement sensors was shared between
the nodes to devise a beam alignment approach with reduced latency. In [52],
the authors considered probabilistic LOS channel models and jointly optimized
the positions of two Unmanned Aerial Vehicles (UAVs) together with the power
control problem for maximizing their sum-rate throughput. The estimation of
the position of a UAV was also leveraged in [53] to devise a matrix completion
algorithm for the recovery of the full mmWave MIMO channel matrix.

It is thus apparent that location/position knowledge has been widely exploited
to enhance XL MIMO channel estimation and reduce the number of training
symbols. In our case, to recover $\mathbf{e}_{m}$ in Step 1 of the proposed
iterative algorithm, knowledge of the initial instantaneous channel vector
$\mathbf{h}_{m}^{(0)}$ is required. From the channel model in (1), the
instantaneous channel matrix $\mathbf{H}$ comprises both the LOS and NLOS
components, which are summed together. Therefore, its $m$-th row of
$\mathbf{H}$ can be expressed as follows:

$\mathbf{h}_{m}\\!=\\!\underbrace{h_{1}\left[\mathbf{a}_{\text{R}}(\theta_{1})\right]_{m}\mathbf{a}_{\text{T}}^{\text{H}}(\phi_{1})}_{\triangleq\mathbf{h}_{m,\rm
LOS}}+\underbrace{\sum_{\ell=2}^{L}h_{\ell}\left[\mathbf{a}_{\text{R}}(\theta_{\ell})\right]_{m}\mathbf{a}_{\text{T}}^{\text{H}}(\phi_{\ell})}_{\triangleq\mathbf{h}_{m,\text{NLOS}}}.$
(44)

In THz MIMO communications, the steering matrices can be retrieved by
leveraging the known positions of the TX and RX [54]. By representing the UE
and BS positions on the 2D-plane as $(x_{\rm UE},y_{\rm UE})$ and $(x_{\rm
BS},y_{\rm BS})$, respectively, the physical AoA and AoD are expressed as
follows:

$\displaystyle\hat{\vartheta}_{1}=\arcsin\frac{x_{\rm UE}}{d_{\text{BS-UE}}},$
(45) $\displaystyle\hat{\varphi}_{1}=\arcsin\frac{x_{\rm BS}}{d_{\text{BS-
UE}}},$ (46)

thus, the normalized AoA and AoD are
$\hat{\theta}_{1}=\frac{d\sin\hat{\vartheta}_{1}}{\lambda}$ and
$\hat{\phi}_{1}=\frac{d\sin\hat{\varphi}_{1}}{\lambda}$. Therefore, the
steering matrices for the LOS component are given by:

$\mathbf{a}_{\text{R}}(\hat{\theta}_{1})=\frac{1}{\sqrt{M}}\left[1,e^{j2\pi\theta_{1}},\ldots,e^{j2\pi(N-1)\theta_{1}}\right]^{\text{T}},$
(47)

and

$\mathbf{a}_{\text{T}}(\hat{\phi}_{1})=\frac{1}{\sqrt{N}}\left[1,e^{j2\pi\phi_{1}},\ldots,e^{j2\pi(N-1)\phi_{1}}\right]^{\text{T}}.$
(48)

Recall from (9) that the binary vector $\mathbf{e}_{m}$ represents the angular
shift due to the beam-squint shift. Hence, it is not affected by variations of
the signal amplitude. Using this observation, the following proposition holds
true.

###### Proposition 2

The estimation of the binary vector $\mathbf{e}_{m}$ is not affected by the
channel gain uncertainty.

Based on Proposition 2, the instantaneous channel gain can be disregarded in
(44), thus, for the initialization of Algorithm 1 the following approximation
can be used:

$\mathbf{h}_{m}^{(0)}=[\mathbf{a}_{\text{R}}(\hat{\theta}_{1})]_{m}\mathbf{a}^{\rm
H}_{\text{T}}(\hat{\phi}_{1}).$ (49)

##  IV Numerical Results and Discussion

In this section, we evaluate the performance of the proposed technique via
computer simulation results using MATLAB${}^{\textrm{TM}}$. The considered
downlink THz communication system consists of a large MIMO UE and an XL MIMO
BS. The distances between the BS and the UE are expected to be larger than the
Fraunhofer distance [55], thus, the THz far-field model is adopted. This
implies that plane waves reach the antenna arrays at both communication ends.

###  IV-A System Setup

TABLE II: Default setting of the simulation parameters. Setting | Value
---|---
Carrier Frequency |  $300$ GHz
System Bandwidth |  $10$ GHz
Number of BS Antennas $N$ | $64$
Number of UE Antennas $M$ | $64$
BS transmit power ($P_{t}$) |  $15$ dBm
Number of multipath components $L$ | $3$
LOS pathloss exponent $\xi_{1}$ | $2$
BS and UE distance $d_{\text{BS-UE}}$ |  $10$ m
NLOS pathloss exponent $\xi_{\ell}$ | $3$
Molecular absorption loss |  $3$dB at $300$GHz[56]

In SC modulation, the BS communicates with the UE using data-carrying frames,
and each frame is composed by $T$ time instances allocated for the training
symbols. Thus, for $t=1,2,\ldots,T$, the BS transmits the vector
$\mathbf{q}(t)=[q_{1}(t)q_{2}(t)\cdots q_{N}(t)]^{\text{T}}$, including the
training symbols, for the targeted joint beam-squint and channel estimation.
Moreover, each time frame has been considered as a new Monte Carlo realization
for all involved random variables (e.g., thermal noise and complex channel
gains). The default channel parameter settings are included in Table II.

The following cases are compared in the evaluation subsection that follows:

  * •

SC idealized case [46]). The two problems of beam-squint (• ‣ III-C) and
channel estimation (34) are being solved independently using the CVX package
[48]. When solving the channel estimation problem, it was assumed perfect
knowledge for the beam-squint vector $\mathbf{e}_{m}$. Respectively, when
solving for the beam-squint estimation, perfect knowledge for the channel
$\mathbf{h}_{m}$ was assumed.

  * •

MC idealized case. This case is similar to the previous idealized case, but
now for MC transmissions. Specifically, due to the beam-squint effect, the
channel changes over different MC frequencies, thus, we focused on the
estimation of the channel for a single subcarrier [57]. This restricts the
number of the available training symbols to the size $M$ of the receiving UE
array.

  * •

The proposed SC approach that solves the joint beam-squint and channel
estimation problem via the iterative Algorithm 1. For the initialization of
the beam-squint problem, the AoA information given by (49) is used. the prior
information on the AoA was considered.

  * •

SC non-idealized case. In this scenario, AoA information is harnessed to
initialize the beam-squint estimation process. The binary vector derived from
this estimation is subsequently applied to solve the channel estimation
problem using CVX package. As a consequence, any inaccuracies or errors
stemming from the beam-squint estimation phase are propagated to the channel
estimation process.

Delving deeper into the initialization process, we make the assumption that
prior knowledge exists concerning the UE position. This prior knowledge
provides an estimate for the products of the steering vectors, specifically
$[\mathbf{a}_{\text{R}}(\hat{\theta}_{1})]_{m}\mathbf{a}^{\rm
H}_{\text{T}}(\hat{\phi}_{1})$ for all values of $m$ in the equation (49). In
an ideal scenario, this information is assumed to be precise and devoid of any
errors. However, in practical situations, it is anticipated that this
information may be subject to noise due to various factors, including
inaccuracies in position estimation and other disturbances, such as phase
noise. An illustration of a XL MIMO THz channel in the beamspace domain is
depicted in Fig. 2 for $L=3$ multipath components. Specifically, the channel
gain in the beamspace is depicted over the dimensions of the $M\times N$
channel matrix $\mathbf{H}$, i.e.,
$\|\mathbf{Z}\|=\|\mathbf{F}_{\text{RX}}\mathbf{H}\mathbf{F}_{\text{TX}}^{\text{H}}\|$.

Figure 2: Illustration of a THz channel with $L=3$ channel paths and $K_{f}=3$
dB.

To evaluate the proposed technique, we employ the normalized mean square error
(NMSE) performance of the estimation of the channel, respectively, i.e.:

$\text{NMSE}\triangleq\sum_{r=1}^{R}\sum_{m=1}^{M}\frac{\|\mathbf{h}_{m}-\mathbf{\hat{h}}_{m}\|}{\|\mathbf{h}_{m}\|},$
(50)

where $\mathbf{\hat{h}}$ is the estimated vector, while $R$ represents the
total number of MC realizations.

###  IV-B Convergence Analysis

The convergence behavior of the proposed technique outlined in Algorithm 1 is
next analyzed. We will evaluate its performance via the NMSE metric across
iterations of proposed Algorithm 1, considering various scenarios with
different Signal-to-Noise Ratios (SNRs). Additionally, we include a benchmark
comparison with the CVX solution, which utilizes the CVX package [48] to solve
the idealized SC optimization problems given by (• ‣ III-C) and (34). Note
that, CVX has constant performance, as it is an one-shot, non-iterative
algorithm.

In Fig. 3, we have set $T=2M$ to investigate the convergence speed of the
proposed channel estimation algorithm, while the values for the other
simulation parameters are given in Table II. The initialization of the
proposed algorithm used the AoA information from (49). In contrast, the CVX
approach has the advantage of being initialized under ideal conditions, where
perfect knowledge of the channel or the binary beam-squint vector is assumed.
This results in swift convergence across all simulated SNRs, necessitating
only a few iterations of the ADMM algorithm with a step size of $\rho=6$.

Significantly, at SNR levels of $10$ and $20$dB, the proposed estimation
method outperforms the CVX solution, firmly establishing itself as a more
resilient approach, especially in low SNR scenarios. This heightened
performance can be attributed to the inherent robustness of the ADMM
algorithm, which consistently converges effectively even when the subproblem
solutions are not perfectly accurate [49].

In contrast, the conventional CVX approach, while benefiting from ideal
initialization, tackles the joint problem in a one-shot manner, making it
unable to effectively handle the MIQP and the introduced channel impairments
stemming from THz communications, such as the molecular absorption loss, the
energy loss due to extended distance between the transmitter and receiver.

Figure 3: Convergence of the NMSE performance for the proposed channel
estimation algorithm with $T=2M$.

###  IV-C Performance Evaluation

Figure 4 illustrates the NMSE performance of channel estimation as a function
of the SNR and the training lengths $T$. In the SNR comparison, we maintained
a fixed training length at $T=3M$, while in the training length comparison a
fixed SNR at $10$dB. The MC-idealized case yields the worst performance, since
for the multi-carrier case the training length is restricted to $M$ symbols,
which is the size of the antenna array. Both the SC idealized and non-
idealized cases showcase identical performance, thereby corroborating the
validity of Proposition 2. This proposition asserts that the channel gain
ambiguity in (49) does not exert any influence on the channel estimation
performance. Remarkably, in low SNR regimes, the proposed SC method even
surpasses the SC-idealized case, aligning with the anticipated outcomes
derived from the convergence analysis.

Figure 4: Channel estimation performance versus the SNR and the training
length $T$.

Next, we consider the distance between the TX and the RX of a THz
communication system. The training length was fixed to $T=2M$ and the SNR to
$10$dB. As shown in Fig. 5, when transitioning from a distance of 10 meters to
30 meters, the NMSE may experience a significant reduction for all the SC
cases. This outcome aligns with expectations, as greater distances between the
TX and RX inherently lead to a reduction in the received signal strength. The
proposed technique is the least affected by the increase of the TX-RX
distance. The training length for the idealized MC is not adequate for
distance greater than 16 meters.

Figure 5: Channel estimation performance versus the BS-UE distance for the SNR
value $30$ dB and the training length $T=2M$. Figure 6: Channel estimation
performance versus training length $T$ with $L=3$ and SNR$=10$dB.

In Fig. 6 we illustrate the NMSE performance as it varies with the training
length $T$ for fixed SNR at $10$dB. In this plot, we include a noisy version
of the initializer in (49), so as to include $1^{\circ}$ noise in the AoA for
$M=64$. Obviously, the proposed SC technique outperforms the CVX-SC case for
all training lengths. Notably, the proposed SC with AoA noise exhibits the
ability to approach convergence levels similar to those without noise,
particularly for larger training lengths exceeding $T>180$ symbols.

It’s worth noting that the impact of AoA noise on the performance of the
proposed algorithm is contingent on several factors. These factors encompass
the dimensions of the antenna array and the number of multipath components.
Specifically, with larger antenna arrays, finer spatial resolution becomes
attainable. Consequently, greater degrees of angle noise can exert a more
substantial influence on the initialization phase of the channel estimation
algorithm. Moreover, the presence of multiple multipath components introduces
a multitude of AoAs, intensifying the susceptibility of the system to the
adverse effects of noise. As depicted in Fig. 6, it’s evident that the
proposed algorithm exhibits the capacity to mitigate the impact of AoA noise.
However, this enhanced noise resilience comes at the cost of requiring a
larger number of training symbols.

##  V Conclusions

In this paper, we addressed XL MIMO channel estimation in the THz frequency
band, considering array-wide propagation delays causing frequency-selective
beam squint. Traditional frequency modulation exhibit high peak-to-average
power ratios, exacerbated by the low THz transmit powers. To confront this, we
presented a novel time-domain SC-based estimation approach, treating beam
squint through sparse vector recovery via optimization. Our technique deploys
alternating minimization to jointly handle the beam-squint effect and the MIMO
channel sparsity. Our performance evaluations confirmed the robustness of the
proposed XL MIMO THz channel estimation approach, as compared to conventional
SC and MC estimation methods. It was also showcased that the proposed
technique achieves comparable performance to the idealized scenario, by solely
utilizing UE prior position information to partially recover the LOS component
of the unknown XL MIMO matrix.

##  VI Acknowledgments

The work of E. Vlachos and G. C. Alexandropoulos has been supported by the
Smart Networks and Services Joint Undertaking (SNS JU) TERRAMETA project under
the European Union’s Horizon Europe research and innovation programme under
Grant Agreement No 101097101, including top-up funding by UK Research and
Innovation (UKRI) under the UK government’s Horizon Europe funding guarantee.
The work of A. Kaushik has been supported by the UKRI Higher Education
Innovation Fund projects “Net Zero and Sustainable 6G: Communications, Sensing
and Computing” and “Green and Autonomous UAVs for Advanced Airborne
Ecosystem”. For the purpose of open access, the authors have applied a
Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript
version arising.

##  VII Data Statement

Data supporting this study are included within the article and/or supporting
materials. The research materials supporting this publication can be accessed
by contacting Evangelos Vlachos (evlachos@athenarc.gr).

## References

  * [1] C. Han, Y. Wang, Y. Li, Y. Chen, N. A. Abbasi, T. Kürner, and A. F. Molisch, “Terahertz Wireless Channels: A Holistic Survey on Measurement, Modeling, and Analysis,” IEEE Commun. Surveys & Tuts, vol. 24, no. 3, pp. 1670–1707, 2022.
  * [2] IMT-2030 (6G) Promotion Group, “6G Vision and Candidate Technologies,” 2021.
  * [3] A. F. Molisch, V. V. Ratnam, S. Han, Z. Li, S. L. H. Nguyen, L. Li, and K. Haneda, “Hybrid beamforming for massive MIMO: A survey,” IEEE Communications Magazine, vol. 55, no. 9, pp. 134–141, 2017.
  * [4] B. Wang, M. Jian, F. Gao, G. Y. Li, and H. Lin, “Beam Squint and Channel Estimation for Wideband mmWave Massive MIMO-OFDM Systems,” IEEE Transactions on Signal Processing, vol. 67, no. 23, pp. 5893–5908, 2019.
  * [5] S. Tarboush, H. Sarieddeen, M.-S. Alouini, and T. Y. Al-Naffouri, “Single- versus multicarrier terahertz-band communications: A comparative study,” IEEE Open Journal of the Communications Society, vol. 3, pp. 1466–1486, 2022.
  * [6] E. Vlachos, G. C. Alexandropoulos, and J. Thompson, “Massive MIMO channel estimation for millimeter wave systems via matrix completion,” IEEE Signal Process. Lett., vol. 25, pp. 1675–1679, Nov. 2018.
  * [7] E. Vlachos, G. C. Alexandropoulos, and J. Thompson, “Wideband MIMO channel estimation for hybrid beamforming millimeter wave systems via random spatial sampling,” IEEE J. Sel. Topics Signal Process., vol. 13, pp. 1136–1150, Sep. 2019.
  * [8] M. Wang, F. Gao, N. Shlezinger, M. Flanagan, and Y. C. Eldar, “A block sparsity based estimator for mmwave massive MIMO channels with beam squint,” IEEE Transactions on Signal Processing, vol. 68, pp. 49–64, 2020.
  * [9] D. Cohen, D. Cohen, Y. C. Eldar, and A. M. Haimovich, “Summer: Sub-nyquist MIMO radar,” IEEE Transactions on Signal Processing, vol. 66, pp. 4315–4330, Aug 2018.
  * [10] M. Rossi, A. M. Haimovich, and Y. C. Eldar, “Spatial compressive sensing for MIMO radar,” IEEE Transactions on Signal Processing, vol. 62, pp. 419–430, Jan 2014.
  * [11] A. Kaushik, E. Vlachos, J. Thompson, and A. Perelli, “Efficient channel estimation in millimeter wave hybrid MIMO systems with low resolution ADCs,” in IEEE European Sig. Process. Conference (EUSIPCO), pp. 1825–1829, Sept. 2018.
  * [12] V. Schram, A. Moldovan, and W. H. Gerstacker, “Compressive Sensing for Indoor THz Channel Estimation,” in 2018 52nd Asilomar Conference on Signals, Systems, and Computers, pp. 1539–1546, 2018.
  * [13] P. Singh, B.-W. Kim, and S.-Y. Jung, “Compressed Detection for Pulse-Based Communications in the Terahertz Band,” Wireless Communications and Mobile Computing, vol. 2018, pp. 1–10, 09 2018.
  * [14] R. Ghods, A. Gallyas-Sanhueza, S. H. Mirfarshbafan, and C. Studer, “BEACHES: Beamspace Channel Estimation for Multi-Antenna mmWave Systems and Beyond,” in 2019 IEEE 20th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC), pp. 1–5, 2019.
  * [15] C. Lin and G. Y. L. Li, “Terahertz Communications: An Array-of-Subarrays Solution,” IEEE Communications Magazine, vol. 54, no. 12, pp. 124–131, 2016.
  * [16] X. Shao, X. Chen, C. Zhong, and Z. Zhang, “Joint Activity Detection and Channel Estimation for mmW/THz Wideband Massive Access,” in ICC 2020 - 2020 IEEE International Conference on Communications (ICC), pp. 1–6, 2020.
  * [17] X. Gao, L. Dai, Y. Zhang, T. Xie, X. Dai, and Z. Wang, “Fast Channel Tracking for Terahertz Beamspace Massive MIMO Systems,” IEEE Transactions on Vehicular Technology, vol. 66, no. 7, pp. 5689–5696, 2017.
  * [18] S. H. Mirfarshbafan, A. Gallyas-Sanhueza, R. Ghods, and C. Studer, “Beamspace Channel Estimation for Massive MIMO mmWave Systems: Algorithm and VLSI Design,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 67, no. 12, pp. 5482–5495, 2020.
  * [19] S. H. Mirfarshbafan and C. Studer, “Sparse Beamspace Equalization for Massive MU-MIMO MMWave Systems,” in ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 1773–1777, 2020.
  * [20] C. Huang, A. Zappone, G. C. Alexandropoulos, M. Debbah, and C. Yuen, “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,” IEEE Trans. Wireless Commun., vol. 18, pp. 4157–4170, Aug. 2019.
  * [21] S. Ma, W. Shen, J. An, and L. Hanzo, “Wideband Channel Estimation for IRS-Aided Systems in the Face of Beam Squint,” IEEE Transactions on Wireless Communications, vol. 20, no. 10, pp. 6240–6253, 2021.
  * [22] A. Kaushik, J. Thompson, E. Vlachos, C. Tsinos, and S. Chatzinotas, “Dynamic RF Chain Selection for Energy Efficient and Low Complexity Hybrid Beamforming in Millimeter Wave MIMO Systems,” IEEE Transactions on Green Communications and Networking, vol. 3, no. 4, pp. 886–900, 2019.
  * [23] I. Gavras, M. A. Islam, B. Smida, and G. C. Alexandropoulos, “Full duplex holographic MIMO for near-field integrated sensing and communications,” in European Signal Proces. Conf., (Helsinki, Finland), Sep. 2023.
  * [24] A. Kaushik, E. Vlachos, C. Tsinos, J. Thompson, and S. Chatzinotas, “Joint bit allocation and hybrid beamforming optimization for energy efficient millimeter wave mimo systems,” IEEE Transactions on Green Communications and Networking, vol. 5, no. 1, pp. 119–132, 2021.
  * [25] G. C. Alexandropoulos, A. Mokh, R. Khayatzadeh, J. De Rosny, M. Kamoun, A. Ourir, A. Tourin, M. Fink, and M. Debbah, “Time reversal for 6G wireless communications: Novel experiments, opportunities, and challenges,” IEEE Veh. Technol. Mag., vol. 17, pp. 74–82, Dec. 2022.
  * [26] A. Faisal, H. Sarieddeen, H. Dahrouj, T. Y. Al-Naffouri, and M.-S. Alouini, “Ultra-massive MIMO Systems at Terahertz Bands: Prospects and Challenges,” IEEE Vehicular Technology Magazine, vol. 15, no. 4, pp. 33–42, 2020.
  * [27] B. Ning, Z. Tian, W. Mei, Z. Chen, C. Han, S. Li, J. Yuan, and R. Zhang, “Beamforming Technologies for Ultra-Massive MIMO in Terahertz Communications,” IEEE Open Journal of the Communications Society, vol. 4, pp. 614–658, 2023.
  * [28] N. Shlezinger, G. C. Alexandropoulos, M. F. Imani, Y. C. Eldar, and D. R. Smith, “Dynamic Metasurface Antennas for 6G Extreme Massive MIMO Communications,” IEEE Wireless Communications, vol. 28, no. 2, pp. 106–113, 2021.
  * [29] G. C. Alexandropoulos and E. Vlachos, “A hardware architecture for reconfigurable intelligent surfaces with minimal active elements for explicit channel estimation,” in Proc. 2020 IEEE Int. Conf. Acoust., Speech, Signal Process.(ICASSP), (Barcelona, Spain), May 2020.
  * [30] G. C. Alexandropoulos, N. Shlezinger, I. Alamzadeh, M. F. Imani, H. Zhang, and Y. C. Eldar, “Hybrid reconfigurable intelligent metasurfaces: Enabling simultaneous tunable reflections and sensing for 6G wireless communications,” arXiv preprint arXiv:2104.04690, 2021.
  * [31] M. Jian, G. C. Alexandropoulos, E. Basar, C. Huang, R. Liu, Y. Liu, and C. Yuen, “Reconfigurable intelligent surfaces for wireless communications: Overview of hardware designs, channel models, and estimation techniques,” Intel. Converged Netw., vol. 3, pp. 1–32, Mar. 2022.
  * [32] R. A. Tasci, F. Kilinc, E. Basar, and G. C. Alexandropoulos, “A new RIS architecture with a single power amplifier: Energy efficiency and error performance analysis,” IEEE Access, vol. 10, pp. 44804–44815, May 2022.
  * [33] R. Singh, A. Kaushik, W. Shin, G. C. Alexandropoulos, M. Toka, and M. Di Renzo, “Indexed multiple access with reconfigurable intelligent surfaces: The reflection tuning potential,” IEEE Communications Magazine, early access, 2023.
  * [34] B. Al-Nahhas, Q.-U.-A. Nadeem, A. Kaushik, and A. Chaaban, “Reconfigurable Intelligent Surfaces-Assisted Communication Under Different CSI Assumptions,” in IEEE International Mediterranean Conference on Communications and Networking (MeditCom) Workshop, pp. 1–6, 2023.
  * [35] E. Vlachos and A. Kaushik, “Subset Selection Based RIS-Aided Beamforming for Joint Radar-Communications,” in 2023 IEEE Wireless Communications and Networking Conference (WCNC), pp. 1–6, 2023.
  * [36] E. Calvanese Strinati, G. C. Alexandropoulos, H. Wymeersch, B. Denis, V. Sciancalepore, R. D’Errico, A. Clemente, D.-T. Phan-Huy, E. De Carvalho, and P. Popovski, “Reconfigurable, Intelligent, and Sustainable Wireless Environments for 6G Smart Connectivity,” IEEE Communications Magazine, vol. 59, no. 10, pp. 99–105, 2021.
  * [37] G. C. Alexandropoulos, K. Stylianopoulos, C. Huang, C. Yuen, M. Bennis, and M. Debbah, “Pervasive machine learning for smart radio environments enabled by reconfigurable intelligent surfaces,” Proc. IEEE, vol. 110, pp. 1494–1525, Sep. 2022.
  * [38] G. C. Alexandropoulos, M. Crozzoli, D.-T. Phan-Huy, K. D. Katsanos, H. Wymeersch, P. Popovski, P. Ratajczak, Y. Bénédic, M.-H. Hamon, S. Herraiz Gonzalez, R. D’Errico, and E. Calvanese Strinati, “RIS-enabled smart wireless environments: Deployment scenarios, network architecture, bandwidth and area of influence,” EURASIP J. Wireless Commun. Netw., to appear, 2023.
  * [39] “TERahertz ReconfigurAble METAsurfaces for ultra-high rate wireless communications (TERRAMETA) SNS JU Project.” https://terrameta-project.eu/.
  * [40] W. Hao, X. You, F. Zhou, Z. Chu, G. Sun, and P. Xiao, “The Far-/Near-Field Beam Squint and Solutions for THz Intelligent Reflecting Surface Communications,” IEEE Transactions on Vehicular Technology, pp. 1–12, 2023.
  * [41] “IEEE Standard for High Data Rate Wireless Multi-Media Networks–Amendment 2: 100 Gb/s Wireless Switched Point-to-Point Physical Layer,” IEEE Std 802.15.3d-2017 (Amendment to IEEE Std 802.15.3-2016 as amended by IEEE Std 802.15.3e-2017), pp. 1–55, 2017.
  * [42] H. Sarieddeen, M.-S. Alouini, and T. Y. Al-Naffouri, “An overview of signal processing techniques for terahertz communications,” Proceedings of the IEEE, vol. 109, no. 10, pp. 1628–1665, 2021.
  * [43] J. M. Jornet and I. F. Akyildiz, “Channel Modeling and Capacity Analysis for Electromagnetic Wireless Nanonetworks in the Terahertz Band,” IEEE Transactions on Wireless Communications, vol. 10, no. 10, pp. 3211–3221, 2011.
  * [44] A. Ben-Tal and A. S. Nemirovskiaei, Lectures on Modern Convex Optimization: Analysis, Algorithms, and Engineering Applications.  USA: Society for Industrial and Applied Mathematics, 2001.
  * [45] S. Diamond and S. Boyd, “CVXPY: A Python-embedded modeling language for convex optimization,” Journal of Machine Learning Research, vol. 17, no. 83, pp. 1–5, 2016.
  * [46] W. P. Adams and H. D. Sherali, “Mixed-integer bilinear programming problems,” Mathematical Programming, vol. 59, no. 1-3, pp. 279–305, 1993.
  * [47] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, “Least angle regression,” The Annals of Statistics, vol. 32, no. 2, pp. 407 – 499, 2004.
  * [48] M. Grant and S. Boyd, “CVX: Matlab Software for Disciplined Convex Programming,” 2014.  Version 2.1.
  * [49] S. P. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers.,” Foundations and Trends in Machine Learning, vol. 3, no. 1, pp. 1–122, 2011.
  * [50] R. D. Taranto, S. Muppirisetty, R. Raulefs, D. Slock, T. Svensson, and H. Wymeersch, “Location-aware communications for 5G networks: How location information can improve scalability, latency, and robustness of 5G,” IEEE Signal Processing Magazine, vol. 31, pp. 102–112, November 2014.
  * [51] G. C. Alexandropoulos, “Position aided beam alignment for millimeter wave backhaul systems with large phased arrays,” in 2017 IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), pp. 1–5, IEEE, 2017.
  * [52] B. Wang, F. Gao, S. Jin, H. Lin, and G. Li, “Spatial- and frequency-wideband effects in millimeter-wave massive MIMO systems,” IEEE Trans. Signal Process., vol. 66, pp. 3393–3406, Jul 2018.
  * [53] G. C. Alexandropoulos, E. Vlachos, and B. Smida, “Position aided beam alignment for millimeter wave backhaul systems with large phased arrays,” in IEEE Asilomar Signals, Systems and Computer Conference, (Pacific Grove, USA), pp. 1–5, IEEE, 2029.
  * [54] A. M. Elbir, W. Shi, A. K. Papazafeiropoulos, P. Kourtessis, and S. Chatzinotas, “Terahertz-band channel and beam split estimation via array perturbation model,” IEEE Open Journal of the Communications Society, vol. 4, pp. 892–907, 2023.
  * [55] C. A. Balanis, Antenna theory: analysis and design.  Wiley-Interscience, 2005.
  * [56] J. Kokkoniemi, J. Lehtomäki, and M. Juntti, “Simple molecular absorption loss model for 200–450 gigahertz frequency band,” in 2019 European Conference on Networks and Communications (EuCNC), pp. 219–223, 2019.
  * [57] E. Vlachos, C. Mavrokefalidis, and K. Berberidis, “Channel estimation for UAV-based mmwave massive MIMO communications with beam squint,” in 2022 30th European Signal Processing Conference (EUSIPCO), pp. 1696–1700, 2022.