Ambit-Process-Based Spatial-Wideband MIMO Channel Model for Sub-THz Urban Microcellular Communication

The design and development of sub-Terahertz (sub-THz) cellular systems entail the need for new channel models that can precisely predict channel characteristics beyond 100GHz in outdoor and dynamic environments. This work proposes a novel multiple-input and multiple-output (MIMO) channel model for cellular communication, developed within the framework of a class of spatio-temporal stochastic processes called ambit-process. The modeling methodology effectively captures the typicalities of sub-THz propagation like molecular absorption and scattering of the evolving multipaths while accounting for the propagation delay of electromagnetic waves across large array apertures deployed at the transmitter and the receiver. This allows for an accurate characterization of the spatial-wideband effect along with other relevant spatio-temporal attributes of the channel. Numerical simulations indicate a good level of agreement between the spectral efficiency and spatio-temporal correlation of the proposed model against a state-of-the-art stochastic Terahertz (THz) channel model and measurements reported in the literature.

and upcoming applications push the current state-of-the-art 5G system to its limits, increasing interests have emerged in the unused part of the radio-frequency (RF) spectrum beyond 100 GHz.More recently, the frequency band ranging from 0.1-0.5 THz has emerged as a promising candidate to enable seamless ultra-high-speed connectivity in beyond 5G systems.Additionally, the world radio-communication conference of 2019 has identified the following bands for mobile and fixed service applications -275-296 GHz, 306-313 GHz, 318-333 GHz, and 356-450 GHz [1].Design and deployment of sub-THz (0.1 THz-0.3THz) communication systems require appropriate channel models which can accurately predict signal statistics at the receiver.The sub-THz and THz propagation are characterized by huge chunks of available bandwidth, high spreading and molecular absorption losses.Presently, the absorption-defined windows around 140 GHz and 300 GHz are being investigated for cellular and point-to-point communications respectively [2].Specifically, the present work aims to provide a comprehensive channel model to address cellular communication in the 140 GHz band.

A. Spatial-Wideband Effect
The sub-THz band (100 GHz -300 GHz) has approximately 21.2 GHz of bandwidth available for wireless broadband [3].Ultra-massive antenna arrays and ultrawide bandwidths are envisioned as key enablers for sub-THz communication.The small propagation wavelength allows ultra-massive arrays to pack hundreds or even thousands of antennas at the basestation [4].Nevertheless, modeling the wireless channel with such a huge number of transmit and receive antennas has not been thoroughly investigated in literature, with most of the massive MIMO channels being direct extensions of the existing MIMO channel models, albeit with an increased channel dimensionality.Physical propagation over a massive array, however, introduces certain effects that can't typically be ignored.If the propagation delay across the array is comparable to or even exceeds the symbol period, different antenna elements would receive different symbol amplitudes and phases or even completely different symbols at the same sampling instance, thereby rendering the conventional MIMO model inadmissible.This inherent property of an ultra-massive array is called the spatial-wideband effect [5].Unlike conventional narrowband beamforming, where the beamforming weights are frequency-independent, the spatial-wideband effect results in frequency-dependent array steering vectors.A key feature of the spatial-wideband effect is its independence of array geometry and orientation.Consequences of the spatial-wideband effect have not been typically addressed in communications since most of the massive MIMO test-beds are not massive enough in any of the dimensions to warrant a consideration of the effect, particularly on beamforming and channel estimation.However, with sub-THz systems employing truly massive arrays, the spatial-wideband effect can no longer be ignored.

B. Related Work
As far as cellular communication is concerned, it is interesting to note that received signal statistics in an outdoor scenario at 140 GHz is highly dependent on the array configurations, propagation bandwidth and geometry of the environment.There are several approaches in the literature to modeling sub-THz channels [3], [6], [7].Based on the modeling methodology, sub-THz channel models may broadly be classified into deterministic, stochastic, and quasi-deterministic or hybrid models [3], [6].Deterministic modeling [8] provides accurate site-specific simulation of a channel based on ray-tracing (RT) but involves complex computations, whereas stochastic modeling offers relatively lower computational complexity; albeit with a compromise on modeling accuracy [9], [10].In deterministic ray-traced models, reflection, scattering [11], [12], and diffraction mechanisms are approximated with geometric optics and propagation losses are computed using simple Friss free space equations.Static indoor multipath channel models based on the ray-tracing techniques for the THz band can be found in [13] and [14].Reference [14], in particular addresses wireless communication over the entire THz range using graphene-based massive MIMO arrays but does not address the spatial-wideband effect dominant in such arrays.Reference [15] reports an extensive measurement campaign to characterize a 2×2 MIMO wideband channel at 275−325 GHz over a bandwidth of 50 GHz.The model, however, caters to indoor scenarios involving propagation distances of just a few meters with no elements of mobility or spatio-temporal broadening.Additionally, RT models are scenario-specific and typically require information on the exact geometry of the propagation environment.
Stochastic models, on the other hand, can be effectively used to statistically model the sub-THz propagation characteristics, based on empirical channel measurements.A stochastic indoor channel model for characterizing propagation attributes at 0.3 THz is developed in [16].This model uses a Saleh-Valenzuela (SV) framework to account for the complete characterization of multipath parameters of a static indoor channel but does not support mobility.Reference [17] extends this to present a hybrid beamformed channel model over short distances using frequency-dependent phase shifters but does not explicitly address the dual-wideband characterization.These models, though fairly accurate in indoor or static settings, fall short in the presence of massive transmitter-receiver arrays or when either the transmitter or the receiver or both are in motion, thus making them unsuitable for urban device-todevice (D2D) or cellular communication.
To trade off computational complexity and modeling accuracy, a hybrid modeling approach may be adopted for simulating sub-THz channels.The stochastic scatterer placement and ray-tracing hybrid approach, simply called SSRTH, is one such method [10], [18], [19], [20].In this setting, scatterers, in the form of discrete points are stochastically placed over the Euclidean plane and the resulting channel impulse response (CIR) and multipath propagation is tracked via simple raytracing techniques.The locations of scatterers are established on the basis of predetermined statistical distributions derived from empirical measurements.References [21] and [22] are two of the recently reported hybrid stochastic models addressing microcellular communication at sub-THz frequencies.While [21] and [23] use highly directive horn antennas to characterize the power-delay-profiles of urban microcellular and D2D channels at 140 GHz, [22] looks at the spectral efficiency of a hybrid beamformed sub-THz channel derived form measurements.The present work adopts a similar hybrid methodology to designing sub-THz cellular channels with the added assumptions of user mobility and the dominant spatialwideband effect.Owing to a more realistic distribution of scatterers [19], SSRTH models usually yield channel statistics that seem to be in better agreement with measurement data [24].

C. Contributions
State-of-the-art channel models provide a comprehensive and reliable approach to simulating cellular channels up to a frequency of 100 GHz.However, existing models lack certain essential features typical of radio propagation at sub-THz frequencies.For example, designing hybrid multipath channels with a proper characterization of the wideband spatial and temporal broadening effects across massive arrays to be deployed at the transmitter and the receiver, has not been fully addressed in the literature.A second essential feature of wideband propagation, mostly absent in existing models is the dominance of the spatial-wideband effect, which, if not addressed properly would lead to severe inter-user interference.While it is essentially difficult to incorporate every salient feature mentioned above into a single channel model, we try to present a novel spatial-wideband massive MIMO channel model which partially addresses the requirements for beyond 5G sub-THz communication in a tractable manner.The proposed model aims to be one of the first stochastic-hybrid MIMO models for microcellular communication at sub-THz frequencies with good spatio-temporal consistency.To summarize, our contributions to this manuscript are as follows: • We propose a novel and generalized massive MIMO channel model for sub-THz cellular communication based on the theory of ambit-process encompassing diverse modeling scenarios over different and possibly arbitrary, scatterer distributions.• The ambit-framework allows for formulating a well-defined and tractable spatio-temporal channel correlation and Doppler that converges to classical scenarios under appropriate assumptions.
• We show that the ambit kernel function may be leveraged to accurately model the spatial-wideband effect experienced by an ultra-wideband sub-THz signal when propagating over a large array aperture.
• Subsequently, we propose a discrete channel simulator that utilizes 2D convolution over a rectangular spatio-temporal grid to simultaneously track the received power, propagation delay, and phase of multipaths of a spatial-wideband MIMO channel in a dynamic environment.
• Finally, we corroborate the validity of the proposed modeling approach by comparing the statistical parameters of the simulated channel with recently reported studies and measurement campaigns.Subsequent sections of the manuscript are organized as follows.Section II presents a primer on the ambit-process framework and its utility in modeling a mobile wireless channel.Section III looks at the kernel characterization of an ambitprocess-based MIMO channel.Section IV characterizes the spatial-wideband effect for sub-THz propagation over large antenna arrays.Section V looks into the spatio-temporal attributes of the proposed model.Section VI proposes a discrete-time channel simulator for sub-THz MIMO channels based on the framework developed in Sections III and IV.The simulation results and related discussions are presented in Section VII followed by the conclusion in Section VIII.

II. A PRIMER ON AMBIT-PROCESS CHANNEL MODELING
We present in this section a brief overview of the ambit-process framework and its utility in modeling mobile wireless channels.A detailed treatment on the general characterization of wireless channels using ambit-processes may be found in [25].
The classical time-varying SISO channel impulse response h(t, τ ) of a wireless channel is given by - The equivalent ambit representation of the classical channel is given by - In ambit-process literature, h(t, y i ) denotes the stochastic field experienced by a receiver at a point (v 0 t, y i ) on the Euclidean plane at time t when the receiver is moving with velocity v 0 .A(t) denotes the subset of dominant scatterers present around the receiver at time t.The stochastic field is nothing but the time-varying impulse response of the mobile sub-THz channel experienced at the receiver.α L and α r (r ∈ {1, 2, .., N R (t)}) denote the time-variant complex amplitudes of the line-ofsight (LoS) and non-line-of-sight (NLoS) components respectively, while τ L (t) and τ r (t) denote the corresponding LoS and single-bounce delays.N R (t) represents the number of scatterers inside the scattering region at time t.k represents the wave number defined by k = 2π/λ c , where λ c is the signal wavelength.The complex amplitudes of the different multipath components in (1a) may be expressed as - where Γ L and Γ r are distance-dependent reference pathgains accounting for various signal power attenuation factors and antenna gains for the LoS and scattered components respectively and are typically obtained from channel measurements.µ(f ) is the coefficient of absorption of the propagation medium at frequency f which has been subsequently defined in (7).κ(t, y i ; u, y j ) represents a deterministic kernel function to capture the received signal power of individual multipath components at a particular time instant t.L(du, dy j ) is called the Levy base of the ambit channel model h(t, y i ; u, y j : (u, y j ) ∈ A(t)).It captures the distribution of scatterers around the receiver on the two-dimensional Euclidean plane.In this work, L(du, dy j ) is assumed to be a homogeneous Poisson Point Process with density λ r = Rr 2Rv0 where R r is the arrival-rate of new multipaths, R is the radius of the scattering zone centralized on the instantaneous location of the moving receiver, and v 0 is the velocity of the MU [25].R can equivalently be called the radius of the ambit set A(t). σ(u, y j ) denotes the stochastic volatility of a scatterer associated with the ambit channel model at time u and position y j .Stochastic volatility captures the small-scale fading along with the localized and random motion of the scatterers about their mean position.In this work, we assume the scatterers to be static, and thus σ(u, y j ) = χ r (u = t) at time t is simply taken to be a complex-valued zero-mean unit variance Gaussian random variable ∀ (u, y j ) ∈ R + × R [25].The stochastic volatility is essentially the ambit equivalent of the small-scale fading statistics χ r (t) defined in the classical model.

III. MIMO AMBIT CHANNEL MODEL
We extend the SISO ambit channel in (1) to its equivalent MIMO representation as follows.Let the base-station is equipped with a uniform linear array (ULA) of M equispaced dipole antennas and the receiver is equipped with a ULA of size N .Then, following (1), the time-varying classical channel impulse response between the m-th receive and n-th transmit antenna is given by - The equivalent ambit representation of the classical MIMO channel is given by - Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
= α L,m,n (t)e jkd L,m,n (t) δ(t − τ L,m,n (t)) + A(t) κ r,m,n (t, y i ; u, y j )σ(u, y j )L(du, dy j ) (3b) where κ r,m,n (t, y i ; u, y j ) represents the deterministic kernel function between the m-th receive and n-th transmit antennas at time t.α L,m,n denotes the time-variant complex amplitude of the LoS component between the m-th receive and n-th transmit antenna and τ L,m,n (t) is the corresponding delay given by - where d L (t) is the LoS distance between the transmit and receive first reference antennas, d B (d u ) and θ r,B (t)(θ r,u (t)) are the inter-antenna spacing and angle-of-arrival (AoA) (angle-of-departure (AoD)) via the r-th path respectively at the receiver (transmitter).In subsequent sections, we characterize this kernel in terms of the angle of departure, angle of arrival, frequency of operation, array size, and multipath delay.

A. Characterization of the Kernel Function
The kernel function essentially captures the received power of the multipath components.Accurate modeling of the kernel function is, thus, of utmost importance for precise modeling of the sub-THz wireless channel.
The instantaneous signal amplitude to the m-th receive antenna from the n-th transmit antenna via the r-th path at any time t can be represented by a time-frequency-dependent function of the form - where is the delay incurred by the reflected ray at time t, Υ(f ) is the frequency-dependent reflection coefficient, d r1 (t) is the time-dependent distance between the first reference antenna of the transmitter and the reflector, and d r2 (t) is the distance between the reflector and the first reference antenna of the receiver, n is the path-loss exponent and µ(f ) is the coefficient of absorption of the propagation medium at frequency f .µ(f ) depends on the molecular composition of the propagation medium and can be essentially characterized as - where T and p denote the absolute temperature and pressure of the propagation medium in Kelvin, p 0 denotes the standard reference pressure, T ST P is the absolute temperature at standard pressure, σ i and Q i denote the frequency-dependent absorption cross-section and the number of molecules per unit volume of the i-th gaseous component present within the propagation medium.For a given frequency f , σ i (f ) can be obtained from the HITRAN database while Q i may be derived from the ideal gas law as follows - where q i is the mixing ratio of the gas i, N A stands for the Avogadro constant and R is the ideal gas constant.The mixing ratio of each gas i in the propagation medium can be taken from the HITRAN database.The absorption coefficient of an outdoor propagation medium is typically dominated by the number of water-vapour molecules present within the medium.
A thorough analysis of the absorption coefficient can be found in [26].Now, consider a characterization of the time-variant mobile sub-THz channel, where the base station is located at (BS x , BS y ) on the Euclidean Plane.Then, at any time t, the location of the receiver moving horizontally along the xaxis with velocity v is given by ( where y i is some constant.Similarly, the coordinates of a typical scatterer1 located within the ambit set of the receiver can be taken to be ( u 0 v(ζ)d(ζ), y r ) for some constants u and y r respectively.Then, The above expressions are true for any velocity v.However, when v = v 0 is uniform, they can be further reduced to - where is the Euclidean distance between points (x 1 , y 1 ) and (x 2 , y 2 ).
The frequency-dependent reflection coefficient defined in (10) accounts for the rough surface reflection loss at sub-THz frequencies which is predominantly dependent on the material properties and shape and roughness of the reflecting surface.Towards this end, Kirchhoff's electromagnetic scattering theory coupled with the Fresnel equations can be used to effectively capture rough surface losses for specular reflection [12].Kirchhoff's theory states that the reflection coefficient of a rough surface can essentially be obtained by multiplying the smooth surface reflection coefficient γ(f ) of the corresponding surface, derived from the Fresnel equations with the respective Rayleigh roughness factor, r(f ), as - The Fresnel reflection coefficient of a smooth surface can be typically approximated as - where θ r is the angle of incidence on the reflector at any time instance t given by - Under the proposed framework, thereby implying that the ambit model typically considers reflected rays with large incident angles.Under such an assumption, Taylor's approximation of the Fresnel reflection coefficient in (11) remains fairly accurate for sub-THz frequencies.η r represents the refractive index of the scattering medium, which typically varies with the frequency and material property of the reflecting surface.A key point to be noted in (11) is the negative sign associated with the Fresnel reflection coefficient which simply indicates a phase shift of π caused by the reflection phenomenon.
The Rayleigh roughness factor of a surface defined in (10) depends on a Gaussian distributed statistical parameter called the rough surface height standard deviation (σ r ) described in [12].The Rayleigh factor is characterized by this roughness effect and expressed as [27] - where J 0 (•) is the zero-order Bessel function of the first kind.
To make the model consistent with the framework proposed in (3), we reformulate the kernel function between the m-th receive and n-th transmit antenna as follows - where Γ r,m,n (t, f ) is the time-frequency-angle dependent function for rough surface reflection derived in (5) and d r,m,n (t) is given by - is the total distance travelled by the reflected ray.Therefore, plugging in the respective values from equations ( 5) -( 13) in (14), we get the complete deterministic kernel function of a reflected ray to the m-th receive antenna from the n-th transmit antenna as - Noting that κ r,m,n (•) is essentially a function of time, frequency and delay, we can, by a slight abuse of notation drop the spatial indices and rewrite (16a) as - or equivalently, where Υ(f ) is as defined in ( 10)-( 12) and n is the path-loss exponent.

IV. CHARACTERIZING THE DUAL-WIDEBAND EFFECT
Assuming that the transmission bandwidth is BW and the sampling period T is chosen such that the Nyquist sampling theorem is satisfied, i.e.T = 1 2BW , it may be observed that the maximum delay difference across the receive antenna array (M − 1)/2f c is negligible w.r.t to the sampling period T if and only if (M −1)/2f c T ≪ 1, or equivalently, iff BW/f c ≪ 1 M −1 .This condition is referred to as the narrowband criteria in literature.It has been shown in [28] that this condition is violated for mm-wave and sub-THz systems which use large transmission bandwidths and massive antenna arrays, whereby the received signal suffers delays of the order of several symbol duration between the first and last antenna elements of the transmit and receive array.The kernel function derived in ( 16) is capable of meticulously capturing this effect as will be evident from the analysis presented below.
Following the kernel function derived in (16a), consider a snapshot of the time-varying channel at a fixed time t.Then κ r,m,n (•) is no longer a function of t and can be rewritten as - where Γ r,m,n (f ) = Γ r,m,n (f )e jk(dr 1 +dr 2 ) is the equivalent complex channel gain of the r-th path.κ r,m,n (τ ) is the spatialdelay kernel function between the m-th receive and n-th transmit antenna.ψ B,r = are the normalized AoA and AoD respectively in the spatial domain.Taking the Fourier transform of (17) - where f ∈ [0, BW ] is the transmission bandwidth and f c is the carrier frequency.Analogous to (17), κ r,m,n (f ) is the spatialfrequency kernel function between the m-th receive and n-th transmit antenna.a B (ψ is the m-th element of the receive array defined by - where (•) * denotes the complex conjugate.Similarly, is the n-th element of the transmit array given by - a B (•) and a u (•) are the frequency-dependent spatial steering vectors of the base-station and the mobile transmitter respectively.
Stacking up the multipaths together, the uplink channel at a given time t between the base-station and the mobile transmitter can be arranged in the form of a matrix given by - where (•) H denotes the conjugate-transpose.The (m, n)-th element of K(f ) denotes the kernel function of the uplink spatial-wideband frequency response between the m−th antenna of the base-station and the n-th antenna of the mobile transmitter.
Consider the scenario where the uplink angular-delay spread at the base-station is observed from the first antenna of the mobile transmitter.Then, plugging n = 1 in (17), the kernel function reduces to - For notational brevity, we can drop the transmit antenna index and rewrite (22a) as - Taking the inverse Fourier transform of (22b) with respect to the spatial domain gives us a representation of the channel kernel function in angular-delay domain as follows - where ψ B,r = (m − 1)ψ B,r /f .Following the channel characterization in (18) and stacking up all the arriving multipaths from the transmitter to the base-station, the equivalent uplink spatial-frequency channel at the m-th antenna of the basestation, as seen by the first antenna of the transmitter, is given by - Assume that OFDM is used as the preferred multicarrier modulation scheme.If the sub-carrier spacing is denoted by η = BW/K, then the spatial-frequency kernel function of the m-th receive antenna for the k-th subcarrier is given by κ m−1 ((k−1)η) and the overall spatial-frequency kernel matrix across the M receive antennas may be expressed as - where Moreover, Φ(ψ B,r ) is a M × K matrix with its (m, k)-th element given by - The channel model described in ( 25) is called the spatialfrequency wideband or the dual wideband channel.It captures the frequency-dependent radiation effect of large antenna arrays operating over sufficiently wideband channels by coupling the spatial domain steering vector a B (ψ B,r ) with the frequency domain steering vector b(τ r ) and the phase-shift matrix Φ(ψ B,r ).It may be noted that the maximum phase shift induced by Φ(ψ B,r ) is approximately BW fc d B λc (M − 1).For small antenna arrays or narrow bandwidth signals, ) is an all-ones matrix, thereby reducing (25) to the classical narrowband channel.A schematic of the ambit-process-based spatialwideband channel has been shown in Fig. 1.

V. SPATIO-TEMPORAL CORRELATION OF THE SPATIAL-WIDEBAND CHANNEL
The ambit-process based spatio-temporal autocorrelation function (ACF) of the channel coefficients at a given Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.propagation delay τ is shown to be [25]- where w = (w t , w s ) ∈ R + × R are small increments in the x and y directions of the 2D Euclidean plane respectively and A(τ ) is the circular ambit set of all scatterers with delay τ at time t centered at the receiver.Equivalently, for an ambit channel characterized by an exponential kernel function over a circular ambit set of radius R, (28) can be rewritten as [25]- κ(t + w t , y + w s )κ * (t, y)dtdy (29)

A. Spatio-Temporal Correlation
Consider a scenario where the coordinates of the scatterer S is (v 0 u, y) for some constants u and y = y r , the initial location of MU at time t is X = (0, 0) and its location at time t+w t is Y = (v 0 w t , 0).Since, the kernel function, κ r,m,n is taken to be deterministic, the temporal ACF would only depend on the time-difference w t .To evaluate the spatio-temporal correlation between the transmit and receive antennas, we normalize the kernel function used in (14) as follows - where κ m,n (t, y) is the normalized kernel function between the m-th receive and n-th transmit antenna at time t.Similarly, the normalized kernel between the p-th receive and q-th transmit antenna at time t + w t is - Assuming a constant velocity, ( 30) and ( 31) can be rewritten as - Therefore, κ p,q (t + w t , y)κ * m,n (t, y) = e jk[d(vwt,0;v0u,y)−d(0,0;v0u,y)] × e jk(p−m)d B sin θ B e −jk(q−n)du sin θu ≈ e −jkv0wt cos θu e jk(p−m)d B sin θ B e −jk(q−n)du sin θu (34) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
We define, κ pq,mn (w t ) = κ p,q (t + w t , y)κ * m,n (t, y) The changes in propagation distance of an incoming ray at the receiver corresponding to an arbitrary scatterer S over a very small time instance w t can be approximated as - where θ is angle made by the S-th scatterer with the x-axis and d(x 1 , y 1 ; x 2 , y 2 ) is the Euclidean distance between points (x 1 , y 1 ) and (x 2 , y 2 ).From elementary geometry, it is clear that θ ≈ θ u , where θ u is the AoD corresponding to the S-th scatterer.Now, representing (36) in polar coordinates, we have -

B. Spatial Crosscorrelation (CCF)
Consider the scenario shown in Fig. 1 where the MU is moving with velocity v = v 0 .Without loss of generality, at any time instance t, let the location of the MU be denoted by (0, 0) and the location of a scatterer is (v 0 u, y = y r ).
To evaluate the spatial correlation between the transmit and receive antennas, we again normalize the kernel functions as follows - where κ r,m,n (t, y) is the normalized kernel for the m-th receive and n-th transmit antenna.Similarly, the normalized kernel function between the p-th receive and q-th transmit antenna is - From elementary geometry, we know that - Plugging ( 43) in (42), we have κ m,n (t, y)κ * p,q (t, y) = e jk(m−p)d B sin θ B e −jk(n−q)du sin θu = e −jk(n−q)du sin θu e jk(m−p)d B sin θu where θ u ≈ θ, v 0 u = r cos θ u , y = r sin θ u .Substituting (44) in (29), the normalized spatial crosscorrelation is - κ m,n (t, y)κ * p,q (t, y)rdrdθ (45)

C. Temporal Autocorrelation (ACF)
We define the normalized temporal kernel functions between m-th antenna of the receiver and n-th antenna of the transmitter as follows - Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
(48) which reduces the normalized temporal autocorrelation function to the following form- e jk[d(v0wt,0;v0u,y)−d(0,0;v0u,y)] Invoking the approximation in (36), (49a) reduces to - where J 0 (•) is the zero-order Bessel function of the first kind.A close inspection reveals that (49b) essentially represents the temporal autocorrelation obtained through Clarke' s 2D isotropic scattering model.Consequently, under the small-scale approximation considered in (36), the ambit 2D scattering model can be further simplified into Clarke's 2D isotropic scattering model.

Remark:
The above discussion pertains to the derivation of the normalized theoretical temporal ACF of a MIMO ambit channel and its relation to Clarke's scattering model.The simulated ACF, on the other hand, can be directly evaluated from the kernel functions using the following relation - where E[•] is the expectation operator and (•) * is the complex conjugate.h m,n (t, τ ) is the ambit channel impulse response between m-th receive and n-th transmit antenna as defined in (3b).The effect of diffused reflection and scattering in the propagation environment is captured in ρ(w t ) through the kernel function κ r,m,n (t, τ ) proposed in ( 16) using the rough surface reflection coefficient Υ(f ) defined in ( 10)- (13).

VI. SIMULATING THE 140 GHZ SPATIAL-WIDEBAND CHANNEL
In this section, we develop a simulation algorithm for an urban microcellular dual-wideband sub-THz channel based on the modeling framework and channel characterizations developed in the preceding sections.The detailed channel simulator is presented in Algorithm 1.The algorithm has been divided into two broad phases -namely the initialization phase and the channel simulation phase.

A. Initialization
The ambit channel model defined in (3b) is essentially a continuous-time representation of a non-stationary channel impulse response experienced at the receiver while it moves over the 2D Euclidean plane.Thus, to effectively simulate the proposed channel, we need to discretize the continuous-time Compute 2D convolution of Y mat and row wise folded X mat .Equate the convoluted sequence to Z i;mat .for each j ∈ {1, 2, . . ., 2D − 1} do 30: 31: end for 33: end for 34: H = [H 1 , H 2 , . . ., H 2(N+P)−1 ] 35: Output: Time-varying channel matrix H ambit-process using suitable approximations.This can be done by sampling the time, space (y-dimension) and delay dimensions into appropriate grids using the initialization parameters P, M, N, and D respectively.The above parameters depend on the maximum propagation delay τ max , delay resolution δτ , time-sampling step size ∆, spatial-sampling step size Λ and simulation time T max through the relations given in Algorithm 1. Lastly, M and N denote the ULA array sizes at the base-station and the mobile transmitter with θ B and θ u being the corresponding AoA and AoD respectively.

B. Channel Simulation
The algorithm leverages two key properties of an ambit-process to efficiently simulate the channel.Firstly, it decomposes the kernel function between the first transmit and first receive antenna κ r,1,1 (•) ≜ κ(•) into two components denoted by κ(•) and κ(•) such that - κ(u, y r ) = Γ r,1,1 (t, f )e jkd(BSx,BSy;v0u,yr) where Γ r,1,1 (t, f ) is as defined in ( 16).(v 0 t, y i ) is the coordinate of the MU and (v 0 u, y r ) is the coordinate of a scatterer within the ambit set at time t.We note that owing to the horizontal motion of the MU, y i and y r are assumed to be constants and d(v 0 t, y i ; v 0 u, y r ) can be approximated by the average distance between the scatterer and the receiver within the ambit set.σ here represents a discrete-time zero-mean unit variance complex Gaussian random variable while ▽ is a binary variable that takes the value 1 if a scatterer is present at a particular point on the spatial-grid within the ambit set and 0 otherwise.When T max is sufficiently small (of the order of a few seconds) it is reasonable to assume that the receiver velocity remains constant as well.Now, the composite kernel function κ(•) can be expressed in terms of κ(•) and κ(•) as [29] - where (*) represents the linear convolution operator.Secondly, the continuous-time ambit field representing the channel impulse response can be effectively discretized and simulated by extending the convolution operation defined in (52) to each sample point of the spatial grid on the 2D Euclidean plane over which the receiver moves [29].
The last part of the algorithm computes the received power due to the LoS component between the first transmit and first receive antenna for each time instance by using the kernel function κ L (.) which is defined as - where being the LoS distance between the first transmitter and the receiver antennas.h ≜ h(t, τ ) gives the time-varying complex baseband channel coefficients between the first transmit and receive antenna pairs.Lines 28-34 sum the arriving multipaths corresponding to different delay taps and combine them with the transmit and receive frequency-dependent spatial steering vectors to generate a dual-wideband channel matrix H i corresponding to each time index i.The final time-varying channel matrix is generated by stacking the H i s together such that H = [H 1 , H 2 , . . ., H Tmax ].

C. Complexity
The algorithm samples the Euclidean plane into an M × N spatial grid.P represents the total number of temporal points at which the channel is simulated and D is the total number of delay taps used.The CIR at each spatio-temporal point is obtained by a 2D convolution of the Levy and kernel matrices.The channel has N transmit and M receive antennas.Therefore, an approximate complexity of the algorithm is given by O((MD(N + P) + M)M N ).

VII. RESULTS & DISCUSSIONS
In this section, we discuss some of the preliminary results obtained using the above channel modeling framework and the implications thereof.A list of the simulation parameters used in the BS-MU channel simulator has been summarized below.We assume the base station to be arbitrarily located at (-100,20), while the initial position of the MU is taken to be (0,0).The mobile terminal moves at a uniform velocity of 40 kmph.The carrier frequency f c is 140 GHz, while the path loss exponent n is taken to be 2 for LoS components and 3.2 for NLoS components [21].We take the sampling time interval ∆ to be 5 µs while the sampling space interval is taken to be Λ = 5∆v 0 /18, the delay resolution δτ = 9.2 × 10 −12 s and the maximum delay τ max = 236 × 10 −11 s.The rough surface standard deviation (σ r ) is randomly taken from the set [10e −6 , 100e −6 , 300e −6 ] m depending on whether the reflecting surface is smooth, slightly rough, or extremely rough respectively.The refractive index η r is also drawn in a similar manner from the set [2.6, 2.1, 1.4] depending on the roughness characterization of the scattering surface.These are typical values that have been drawn from measurements carried out in literature for characterizing surface roughness [27].The atmospheric absorption coefficient µ(f Fig. 2 depicts the dual-wideband effect at the base-station receiver array as seen from the first transmit antenna at the mobile terminal with Fig. 2(a) being a three-dimensional representation of the channel impulse in the angular-delay domain.The effect, though, is fairly general in nature and may be observed from any antenna at the transmitter.We assume that the transmitter and receiver are equipped with ULAs of 128 dipole elements each, that are separated by halfwavelengths.We choose OFDM as the preferred multicarrier scheme and assume a transmission bandwidth of BW ≈ 0.2f c with a subcarrier spacing of ∆f = BW/K, K = 128 being the number of subcarriers.The above parameters have been used for better visual representation as well as to stress the severity of the dual-wideband effect in ultra-wideband systems over large antenna arrays.Therefore, we choose to keep the fractional bandwidth (BW/f c ) at a high value of 0.2.It can be shown that this effect is equally pronounced in systems Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.with "somewhat more realistic" configurations.We elaborate on this below.Recall that Eq. ( 25) represents the channel in the spatial-frequency domain.The channel kernel function can be equivalently represented in the angular-frequency domain through the transformed kernel F H M K as in Fig. 2(c), spatial-delay domain through the transformed kernel KF H K as in Fig. 2(d) and in the angular-delay domain through the transformed kernel F H M KF H K as in Fig. 2(b), where F M and F K are normalized DFT matrices of dimensions M and K respectively.Such a transformation maps the AoA and AoD belonging to the continuous intervals [−π/2, π/2] to the discrete interval [0, M − 1] and the path delay belonging to the continuous interval [0, KT ] to the discrete interval [0, K − 1] with M and K representing the angular and delay domain indices respectively.Fig. 2(d) illustrates the corresponding path delays as seen by the different antennas in the spatial domain for a sparse channel having just four multipaths.The antennas see a distinct set of progressive delays for a given multipath.Similarly, different frequency sub-bands see a distinct sequence of progressive angle-ofarrivals corresponding to a given multipath.These effects have been termed as delay squints and beam squints respectively in literature [31].When the beam and the delay squints are coupled together and observed in the angular-delay domain, we get the square regions shown in Fig. 2(b).Each of the four squares corresponds to the angular-delay diffusion of each of the multipaths arriving at the receiver.
Without the spatial-wideband effect, the phase-shift matrix Φ(•) defined in (27) vanishes and each path manifests as an impulse in the angular-delay domain.In presence of the spatial-wideband effect, however, the energy of each of these impulses spills over and diffuses into a square grid as shown in Fig. 2(b).Moreover, the angular and delay squint in Figs.2(c) and 2(d) can be proven to be identical, with each being the side length of the square or equivalently the propagation delay of each path in time-domain samples.The diffusion due to the r-th multipath in the angular-delay domain may be approximately expressed as When the physical parameters like signal wavelength, inter-element spacing and AoA remain unchanged, the diffusion due to the r-th multipath in the angular-delay domain at the receive array is effectively given by d r which is a function of the array size and the fractional bandwidth.Thus, it may be inferred that a channel with a high fractional bandwidth (≈0.2) and moderate array size (of ≈128 elements) will have the same diffusion as a channel with low fractional bandwidth (≈0.025) and massive array size (of ≈1024 elements).Sub-THz systems are envisioned to have bandwidths of upto 4 GHz ( BW f C ≈ 0.03) and employ massive arrays of 1024 × 1024 configuration [4] thereby resulting in  severe angular-delay spreads in both the azimuth and elevation domains at the transmitter and the receiver.To put this into perspective, consider a scenario for LTE systems with a typical transmission bandwidth of 20 MHz over a carrier frequency of 2.6 GHz using a 128-element ULA with half-wavelength spacing.For such a system, the maximum propagation delay across the array comes out to be 0.489T where T is the symbol duration.For a typical sub-THz system using 1024 element ULA over 4 GHz bandwidth at a carrier frequency of 140 GHz, this delay increases to 14.614T , thereby causing significant inter-user-interference.Further, it is interesting to note that the multipath channels are asymptotically orthogonal to each other in the angular-delay domain [32] even if their diffusion regions partially overlap with each other.Exploiting this orthogonality may allow the base-station to distinguish between users at different locations or even users at the same location with different delays by separating their signals through filtering in the angular-delay domain thereby enabling simultaneous communication with multiple users without inter-user-interference.
Figs. 3 and 4 compare the simulated and theoretical autocorrelation functions (ACF) and Doppler of the ambit channel model with the TeraMIMO channel model recently proposed in literature [30].The solid curve shows the channel autocorrelation at the first antenna of the receiver as seen from the first antenna of the transmitter while the dashed curve captures the channel autocorrelation at the last receive antenna.The slight bumps in the dashed autocorrelation curve are due to the additional distance travelled by an impinging signal from the first receive antenna all the way upto the last one.This added distance, given by (M − 1)d B sin θ B , is significant when the narrowband assumption is violated.
It has been shown in (49b) that the normalized ACF of the ambit channel can be reduced to Jake's classical model under certain small-scale approximations.The normalized temporal ambit ACF is roughly equal to the temporal autocorrelation obtained through Clarke's 2D isotropic scattering model.Therefore, by extension, the corresponding ambit Doppler should be approximately the same as Jake's spectrum obtained from Clarke's model.Figs. 3 and 4 essentially capture this similarity in terms of the ACF and Doppler plots respectively.The simulated ambit ACFs closely follow Clarke's theoretical ACF but deviate slightly from the TeraMIMO model.The autocorrelation functions exhibit an approximate cyclic decay behaviour.This decay is due to the decorrelation of the arriving multipaths.However, since at 140 GHz, the time-varying channel experiences excessive diffuse reflection and scattering [11], the simulated ambit ACF slightly deviate from the TeraMIMO model.A similar conclusion can be drawn from Ambit and TeraMIMO Doppler plots shown in Fig. 4. Although, both the plots show a maximum Doppler of f d = 5.3 KHz, the TeraMIMO Doppler may be thought of as an asymptotic envelope to the Ambit Doppler i.e. if we average out the Doppler over all possible scatterer distributions, the corresponding ambit envelope would then closely follow the TeraMIMO model.This deviation is further augmented by the approximate assumptions on the scatterer distribution, the discretization of the ambit-process, and the 2D approximation of the scattering environment.
As a result of diffuse reflection, an arbitrary and excessive number of multipaths arrive at the receiver from different directions even if the number of reflectors within the ambit set remains typically low.A significant portion of these paths does not contribute any meaningful power to the receiver by virtue of the high path-loss and molecular absorption associated with sub-THz propagation.The effect of diffused reflection is further augmented by the motion of the receiver and the surface roughness of the reflectors within the ambit set.The simulated ACFs capture the effect of diffused reflection and scattering through the rough surface reflection coefficient Υ(f ) defined within the kernel function.It is to be pointed out that although there might be slight deviations between the absolute values of the ambit and TeraMIMO ACFs, a good agreement between the locations of their successive nulls and peaks is sufficient to characterize the correlation time of the channel.The correlation time is the time at which the temporal autocorrelation falls to 1/e or 0.37.As can be seen from the plots, the first null, locations of successive nulls and peaks, and the correlation time of both the ambit and TeraMIMO ACFs are well in agreement at low and high receiver velocities thereby reinforcing the validity of the proposed modeling approach.
Fig. 5 compares the simulated and theoretical spatial CCF at the receive antenna array as seen from an antenna at the transmitter.Note that the spatial correlation decays exponentially as the inter-element spacing increases.An important characteristic of the channel here is the correlation distance which is defined as the distance across the array at which the spatial crosscorrelation falls to 1/e ≈ 0.37.The coherence distance may be taken to be slightly less than λ c /2.While the theoretical correlation decays to almost zero, the simulated CCF still shows some residual correlation.This can be attributed to excessive diffused reflection and scattering, the motion of the mobile transceiver, non-stationary properties of the simulation environment and discretization of the ambitprocess.However, it is to be noted that although there might be slight deviations between the simulated and theoretical CCFs, a good agreement between their first trough is sufficient to characterize the correlation distance of the channel.
Figs. 6(a) and 6(b) compare the spectral efficiency as a function of average SNR per receive antenna of a sub-THz microcellular channel under different beamforming and spatial multiplexing schemes.In this work, we use the terms spectral-efficiency and channel capacity interchangeably [22].We compare the information capacity of the proposed model to the ones derived from measurements and reported in [22].The mutual information in a time-varying MIMO channel with N transmit and M receive antennas at any time t can be expressed by- where P is the total transmit power, ℵ 0 is the total noise power, H t is the M × N normalized channel matrix, I is the N s × N s identity matrix with N s being the number of transmitted spatial data streams.Strictly speaking, since the input source distribution is not optimized, C does not truly represent the channel capacity.However, for closed-loop spatial multiplexing using eigen-beamforming, where the precoding matrix is chosen from the right singular vectors of H t , the source distribution is actually optimized such that C becomes the channel capacity.
Employing eigen-beamforming and assuming perfect channel knowledge to evaluate the spectral efficiency of any channel model is an accepted practice in literature [22], [33].Nevertheless, for more practical scenarios where perfect channel state information (CSI) cannot be acquired, training-based channel estimation schemes may be preferred due to their innate ability to track the temporal variation of the wireless channel in any sub-THz radio access system.
To generate the plots in Fig. 6, we assume a 16 × 8 uniform rectangular array (URA) at the base-station and a 4 × 2 URA at the mobile transmitter.The array sizes and configurations are chosen to mimic the measurement scenario depicted in [22].The total transmit power is taken to be 0 dBm which remains constant irrespective of the number of spatial data streams.We further assume perfect channel knowledge at the transmitter such that the transmit precoding and receive combining vectors are chosen from the right and left singular vectors of the channel matrix H t obtained via singular value decomposition (SVD).In presence of rich multipaths, spatial multiplexing with multiple data-streams should result in increased spectral efficiency as compared to beamforming.This, however, may not always be the case for outdoor mm-wave and sub-THz communication owing to a sparse and ill-conditioned channel matrix where transmit power can only be focused in a few specific directions.The lower degrees of freedom of the outdoor sub-THz channel may restrict the usage of higher-order spatial multiplexing schemes.Both the simulated and measured plots indicate that beamforming outperforms spatial multiplexing at low SNR regions.The SNR thresholds at which spectral efficiency for spatial multiplexing improves are 17 dB for the NLoS and 20 dB for the LoS scenarios respectively.This indicates that the LoS component is significantly stronger than the reflected multipaths and unless SNR is above a certain threshold, it's  beneficial to direct the data-stream along the strongest path via beamforming.Additionally, the high beamforming gains of the massive antenna arrays at the transmitter and receiver would sufficiently help in combating the significant path-loss and molecular absorption associated with sub-THz propagation.
Finally, it is interesting to note that the spectral efficiency of a spatial multiplexing scheme with two data streams is higher than three or four streams.This may be accorded to the fact that although a typical outdoor sub-THz propagation environment may have four or five spatial clusters, the negligible amount of power associated with some of these clusters may not always allow for reliable data transmission via higher-order spatial-multiplexing. Lastly, it is to be noted that although (54) represents the spatial-wideband capacity of the channel, it shows good agreement with the narrowband capacity reported in [22] due to the narrow transmission bandwidth (≈1 GHz) and small transmit and receive antenna arrays, thereby, reinforcing the assumptions made in Sections III and IV.
Figs. 7(a) and 7(b) compare the information capacity of the proposed channel model against spatial correlation at the transmit and receive antenna arrays.We employ the widely accepted exponential correlation model [34] to understand the capacity degradation of a typical outdoor sub-THz channel under significant spatial correlation between adjacent antenna elements.The receive correlation matrix is defined as -R rx [i, j] = ρ j−i , i ≤ j ρ i−j , i > j , i, j ∈ {1, 2, . . ., M } (55) where R rx [i, j] is the (i, j)-th element of the receive correlation matrix R rx , and ρ ∈ [0, 1] is the correlation coefficient between adjacent antenna elements.The transmit correlation matrix R tx can be characterized in a similar manner.
As the correlation coefficient ρ increases from 0 to 1, a noticeable degradation in the spatially multiplexed channel capacity is observed for both LoS and NLoS scenarios.This degradation may be accorded to the fact that a high spatial correlation between the transmit and receive antennas reduces the rank of the channel matrix thereby reducing the number of available eigen-channels which results in reduced channel capacity through spatial multiplexing.The capacity of the beamformed channel is mostly unaffected since data streams are directed along the strongest available path and the reduced channel rank has little to no effect on this path.Although  spatial multiplexing achieves higher capacity in an uncorrelated channel, beamforming may eventually outperform spatial multiplexing in terms of spectral efficiency as the spatial correlation increases.Thus, even when the SNR per receive antenna is high, if the propagation scenario suffers from high spatial correlation, beamforming may be preferred to spatial multiplexing.
Figs. 8(a) and 8(b) show the degradation in spectral efficiency of the spatial-wideband channel vis-a-vis the number of receive antennas when conventional narrowband beamforming is employed at the receiver in both high and low SNR scenarios.The capacity is evaluated using (54) with the added assumption that the BS antenna array vector is frequency-dependent while the UE steering vector is frequency-independent.The solid diamond curve represents the increase in channel capacity as a function of the number of receive antennas when using eigen-beamforming via SVD assuming perfect channel knowledge.However, for a spatial-wideband channel, it can be seen that spectral efficiency remains almost unchanged when narrowband beamforming is used even if the number of receive antennas is increased.This is due to the detrimental effect of ignoring the spatial-wideband phenomenon and combining the frequencydependent spatial steering vectors of the base-station with the frequency-independent array vector at the receiver and viceversa.In such a scenario, following (21), the narrowbandbeamformed spatial-wideband channel at time t would take the following form - Neglecting the spatial-wideband effect at the receiver array vector a u (•) results in significant capacity degradation.Moreover, under the same channel conditions, as the fractional bandwidth increases, the spectral efficiency deteriorates further.
Figs. 9(a) and 9(b) show the degradation in spectral efficiency of the spatial-wideband channel against an increasing fractional bandwidth (BW/f c ) when using conventional narrowband beamforming for both high and low SNR values.With eigen-beamforming and perfect CSI, the spectral efficiency remains unchanged even for high fractional bandwidths.This, however, is not the case for narrowband beamforming, which suffers heavily in terms of capacity degradation.For example, a receiver employing a massive ULA with 1024 elements and BW/f c = 0.5 almost degrades to a SISO channel in terms of capacity, if the spatial wideband effect is not properly addressed at the base-staion and user terminals.In conclusion, the higher the number of array antenna elements and fractional bandwidth of the propagating waveform, the higher the loss in spectral efficiency of the narrowband-beamformed spatialwideband channel.

VIII. CONCLUSION
In this work, a novel stochastic hybrid MIMO channel model for addressing sub-THz cellular communication in urban microcells with a specific focus on the 140 GHz band has been developed.The model leverages the spatio-temporal framework of a class of stochastic processes called ambit-process to track the random evolution of a cellular channel in an urban environment.We derived a modified ambit kernel function that meticulously captures the physical propagation delay of the electromagnetic wave across large array apertures to be deployed for sub-THz propagation resulting in the spatial-wideband effect.A complete statistical characterization of the proposed MIMO channel in terms of the spatio-temporal auto and cross-correlation functions has been performed.Subsequently, we employed tools from the ambit-process literature to design a discrete convolutionbased spatial-wideband channel simulation algorithm.Finally, we concluded by presenting and discussing some of the preliminary results on the spectral efficiency and spatio-temporal correlation of the proposed channel and benchmarking the results against a state-of-the-art stochastic THz channel model as well as measurements reported in the literature.
T and (•) denotes the Hadamard product.b(•) is called the frequency domain steering vector pointing towards the r-th delay τ r at the receiver as observed from the first antenna of the transmitter and is given byb(τ r ) = [1, e −j2πητr , . . ., e −j2π(K−1)ητr ] T ∈ C K×1(26)

Fig. 1 .
Fig. 1.A snapshot of the ambit-process-based spatial-wideband channel model over a 2D Euclidean grid.

Fig. 2 .
Fig. 2. Manifestations of the dual-wideband effect in a four-path sparse sub-THz MIMO channel in different domains: (a) magnitude-angular-delay domain (b) angular-delay domain (c) angular-frequency domain, and (d) spatial-delay domain.

Fig. 6 .
Fig.6.Spectral efficiency for different number of data streams in LoS and NLoS scenarios[20].

Fig. 8 .
Fig. 8. Capacity degradation of the spatial-wideband channel as a function of receive antenna index when employing conventional narrowband beamforming.

Fig. 9 .
Fig. 9. Capacity degradation vs fractional bandwidth of the spatial-wideband channel for different MIMO configurations when using Eigen (with perfect CSI) and narrowband beamforming respectively.