Ambit-Process Based Channel Model for Urban Microcellular Communication at 140 GHz

The design and development of Terahertz (THz) and sub-Terahertz (sub-THz) communication systems entail the need for new channel models that can precisely predict channel attributes at such frequencies ( $\geq 100$ GHz) in outdoor and dynamic environments. This work proposes a novel hybrid-stochastic ultra-wideband channel model for sub-THz bands, developed within the framework of a class of spatio-temporal stochastic processes called the ambit-process. The proposed model is capable of supporting bandwidths of upto 1GHz. The spatio-temporal evolution of the ambit framework allows for a spatially consistent, reasonably accurate and tractable characterization of the fading statistics and multipath propagation of the cellular channels. We leverage a recently proposed convolution-based low-complexity algorithm with necessary modifications to study key features of the microcellular sub-THz channel like associated diffused reflection and scattering, molecular absorption, spatio-temporal correlations, and consistency between the time-evolving delay and Doppler of the multipaths. Simulation results on path loss, shadowing, delay spread, and channel correlations indicate that the ambit model accurately captures the typicalities of an urban microcellular sub-THz channel and agrees well with the measurement results reported in the literature.

Gigahertz with each frequency window having its own unique propagation characteristics.The World Radiocommunication Conference 2019 (WRC-19) envisioned different parts of the THz band to be used and studied for different use cases since no one band can cater to all possible scenarios.The present work deals with proposing an acceptable channel model for urban microcellular communication at sub-THz frequencies beyond 100 GHz.The D band (110 GHz-170 GHz) has been considered by WRC-19 as a promising candidate in this regard [1].Recent works, specifically, [2] have shown that sub-THz frequencies, mainly in the D-band centered around 140 GHz may be used for cellular communications for distances up to 200 m in typical urban microcells.
Existing models provide a comprehensive approach towards simulating reliable channels up to a frequency of 100 GHz [3].However, they lack certain essential attributes typical of radio propagation at sub-THz frequencies.A key feature of sub-THz propagation is the dominant role of scattering and atmospheric absorption in addition to rough-surface reflection [4].While cellular systems operating in the microwave and mm-Wave bands can neglect these effects, they cannot be ignored for the THz and sub-THz bands.
While it is essentially difficult to incorporate all the salient features mentioned above into a single channel model having a moderate complexity, we try to design a low-complexity stochastic ray-traced model that partially addresses the requirements for beyond 5G sub-THz communication.The proposed model aims to be one of the first spatially consistent stochastic-hybrid models for microcellular communication at sub-THz frequencies.
We extend a newly proposed framework based on ambitprocesses introduced in [5] to model the sub-THz wireless channel.The ambit framework is quite general in nature and capable of modeling a variety of geometry-based stochastic channel models (GBSM)s including the popular one-ring and COST models [6], [7].The advantage of this modeling approach is that its formulation allows many of the existing GBSMs to be brought under a single mathematical framework having a rich collection of analytical tools useful for studying their underlying characteristics.The proposed approach can, therefore, provide new insights into existing GBSMs.Towards this end, an equivalence proof between the ambit and the COST259 channel model [7] is given in Appendix A.
Essentially, an ambit-process models the temporal effect of stochastic fields exerted by certain random elements in the 1536-1276 © 2023 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
environment on a point located on the Euclidean plane.From a communication perspective, the stochastic field is analogous to the signal statistics at the receiver under the influence of fading, mobility and random environmental scattering.The stochastic field experienced at the receiver can be generally modeled through a combination of deterministic and stochastic functions.More on this can be found in Section III.It may be observed that many of the existing GBSMs generate large-scale parameters (LSPs) based on their probability distribution for each channel realization.Subsequently, the synthesized LSPs are used to simulate the channel impulse response (CIR) for each time instance.This approach, however, introduces modeling errors while simulating a larger movement of wireless nodes beyond the autocorrelation distance [7], [8].For cellular sub-THz systems, the correlation distance has been measured to be sufficiently small [9].Thus, significant modeling discrepancies may arise while simulating a mobile sub-THz channel using such approaches.A special class of advanced GBSMs get around this limitation by adopting specific spatial consistency models to generate power delays and angle-of-arrivals such that the generated parameters are temporally and spatially consistent with the receiver motion at each realization.These are termed non-stationary GBSMs and are designed to accommodate wireless nodes with significant movement beyond the limits of the autocorrelation distance.The QuaDRiGa [10] and IMT-2020 [11] are two notable models belonging to this class.These models are, however, limited to frequencies up to 100 GHz.Consequently, the CIRs of a sub-THz channel generated in the ambit framework are derived from the continuous description of the Levy basis, which is also independent of the movement of the wireless nodes.The ambit model, therefore, has a spatial consistency similar to non-stationary GBSMs and can easily extend channel simulations beyond the autocorrelation distance.
Lastly, most of the existing channel models are either too simple to recreate realistic propagation scenarios or too complicated to merit parametric investigations over multiple long runs.The 3GPP channel model has been shown to compromise on channel properties when improving on simulation time and vice-versa [12].The ambit framework on the other hand does not suffer from this shortcoming and guarantees a polynomial complexity for simulating time-varying channels [5].Moreover, [5] has shown that the ambit framework is capable of handling modelling scenarios with high-mobility-induced diffused reflections even in the microwave and mm-wave bands.Likewise, the THz and sub-THz bands experience severe scattering in addition to diffused reflection, thereby making the ambit framework a natural tool to model these scenarios quite effectively with reasonable complexity.
To summarize, our contributions to this manuscript are as follows: • We propose a novel generalized outdoor sub-THz channel model based on the theory of ambit-process encompassing diverse modeling scenarios that are able to simultaneously track the received signal amplitude, propagation delay, and phase of multipaths with polynomial complexity over different, possibly arbitrary, scatterer distributions.
• The modeling approach allows for formulating well-defined and tractable spatio-temporal correlation and Doppler functions.Furthermore, we show that the ambit channel model is equivalent to Clarke's fading model [13] under isotropic scatterer distributions.• We then present a channel simulator based on a Poisson point ambit-process by leveraging a simple two-dimensional convolution operation over a spatial grid.Additionally, we characterize the kernel functions of the ambit framework to reflect the different physical phenomenons that are associated with sub-THz propagation, namely diffused reflection, scattering, molecular absorption and vegetative losses.
• Finally, we corroborate the validity of the proposed modeling approach by comparing different metrics of the simulated channel like path-loss and shadowing, delay spread, and channel correlation with recently reported studies and measurement campaigns.Subsequent sections of the manuscript are organized as follows.Section II looks at a brief overview of related works and prior attempts to model THz and sub-THz channels for specific scenarios.Section III presents a primer on the ambit-process framework and its utility in modeling a mobile wireless channel.Section IV proposes a low-complexity channel simulator for sub-THz communication based on the framework developed in Section III.The simulation results and related discussions are presented in Section V. Parameter validation for the proposed model is discussed in Section VI followed by the conclusion in Section VII.

II. RELATED WORK
There are several approaches in the literature to modeling THz and sub-THz channels [4], [14], [15].Based on the modeling scenario, a channel model may typically be classified as an indoor-channel [16], an outdoor-channel, a short-range channel [17] or a nanoscale channel [18].A large number of channel models are available for indoor, static or short-range communication scenarios [19].Outdoor models, on the other hand, are mostly limited to higher frequencies (≥300 GHz) and point-to-point backhaul links, with no significant advances having been made to date on modeling channels for cellular sub-THz communication [20].Measurement campaigns at 300 GHz [21] showed that outdoor communication in this band is strictly limited to about 30 m and therefore, may not be suitable for cellular communications.Furthermore, reliable communication at 300 GHz depends on the availability of lineof-sight (LoS) links [21].This condition, though is violated in urban microcells where propagation may be dominated by non-line-of-sight (NLoS) components [9], [22].The 140 GHz band, however, doesn't suffer from these limitations [9], [22].Notably, the empirical 140 GHz channel model developed by the European HEXA-X project [23] is currently restricted to short-range communications at over 10-15 m with long-range (100-150 m) measurements underway.The lack of spatially consistent stochastic channel models for cellular communication at 140 GHz motivates our present effort.
Based on the modeling approaches proposed in the literature, channel models for sub-THz communications may broadly be classified into deterministic, stochastic, and quasi-deterministic or hybrid models [4], [14].Deterministic modeling [24] 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 [25], [26].In deterministic ray-traced models, reflection, diffraction and scattering mechanisms are approximated with geometric optics and propagation losses are computed using simple Friss free space equations [27].A static indoor multipath channel model based on ray-tracing techniques for the entire THz band can be found in [19].An extensive measurement-based RT channel model for high-speed scenarios has been proposed in [28] as part of a study on the viability of THz communication-enabled smart-rail mobility.The model however is designed for the 300 GHz band.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 different propagation attributes at sub-THz frequencies, based on empirical channel measurements.A stochastic indoor channel model for characterizing radio propagation at 300 GHz is developed in [29].This model accounts for the complete characterization of multipath parameters of a static indoor channel but does not support mobility.Additionally, a static cluster-based stochastic indoor channel model parameterized through extensive measurements at 28 GHz and 140 GHz has been developed in [30].These models, though significantly accurate in indoor or static settings, fall short when either the transmitter or the receiver or both are in motion, thus making them unsuitable for urban 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 [26], [31], [32], [33], [34].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 ray-tracing techniques [14].The locations of scatterers are determined on the basis of predetermined statistical distributions derived from empirical measurements.Owing to a more realistic distribution of scatterers [32], SSRTH models usually yield channel statistics that seem to be in better agreement with measurement data [35].This is apparent from the channel results presented in [36] where the RT-stochastic model is parameterized using measurements conducted in an indoor office building at 130-140 GHz.The SSRTH approach has been extensively used in the COST-IRACON channel models to effectively capture mm-Wave V2X propagation in an urban environment [37].Coincidentally, the COST-IRACON model for THz communication is presently limited to high data-rate long-range point-to-point backhaul links at 275-450 GHz or short-range indoor communications [38].A variation of the SSRTH approach is the statistical and ray-tracing hybrid (SRH) approach [14], in which dominant multipath components such as LoS and specular reflections can be simulated using ray-tracing techniques and the remaining multipaths via a stochastic approach [39].The ambit model falls within this category.For a more holistic overview of the modeling aspects and measurement methodologies of THz wireless channels in general, the readers are referred to [27] and [34].
Consequently, existing models may not be trivially extendable to higher frequencies [4].Most stochastic hybrid raytracing models are designed for frequencies below 100 GHz.Therefore, they do not incorporate the significant effects that scattering or molecular absorption has on signal propagation at higher bands.Alternatively, hybrid ray-tracing models in the THz and sub-THz range have mostly been designed either for different frequency bands, point-to-point links or short-range and indoor scenarios.Therefore, although they incorporate the effects of scattering and molecular absorption in calculating the channel impulse response, the absence of mobility makes them unsuitable for modeling cellular channels.
Extensive indoor and outdoor measurement campaigns are currently ongoing in the THz and sub-THz bands to determine the necessary modifications needed for the existing models to capture the unique propagation attributes at these frequencies.A comprehensive overview and comparison of the large-scale indoor path-loss models at mm-Wave and sub-THz bands can be found in [40].It highlights the key similarities and differences between mm-Wave and sub-THz propagation in terms of the path-loss exponent and multipath time dispersion.References [9], [22], and [41] have recently reported their initial results on channel measurements which are specifically aimed towards conducting a feasibility study for urban microcellular and D2D communication at 140 GHz.This merits the timely development of new spatially consistent hybrid stochastic models to complement existing state-of-theart empirical models in this frequency band.

III. A PRIMER ON AMBIT-PROCESS BASED 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 characterization of wireless channels using ambit-processes has been thoroughly discussed in [5].
Considering both specular reflection and diffuse scattering to be dominant phenomenons, the classical time-varying channel impulse response h(t, τ) of a sub-THz wireless channel is given by - Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The equivalent ambit representation of the classical channel is given by - κ r (t, y i ; u, y j )σ (u, y j )L(du, dy j ) κ s (t, y i ; u, y j )σ (u, y j )L(du, dy j ) 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 along the x-coordinate.A s (t)(A r (t)) denotes the subset of dominant scatterers (reflectors) present around the receiver at time t.Variables u and y j in the expression h(t, y i ; u, y j : (u, y j ) ∈ A r (t) ∪ A s (t)) are dummy variables that capture the time and spatial variation of the ambit channel imposed by the subset of dominant scatterers A s (t) and reflectors A r (t) around the receiver at any time t.The stochastic field describes the time-varying channel impulse response of the mobile sub-THz channel experienced at the receiver.The first term in (1a) corresponds to the LoS component of the time-varying channel while the subsequent terms correspond to the single bounce multipaths from the surrounding scatterers and reflectors.α L , α r (r ∈ {1, 2, .., N r (t)}), and α s (s ∈ {1, 2, .., N s (t)}) denote the time-variant complex amplitudes of the LoS, reflected and scattered components respectively, while τ L (t), τ r (t), and τ s (t) denote the corresponding delays.N s (t)(N r (t)) represents the number of scatterers (reflectors) inside the ambit-set at time t.k represents the wave number defined by k = 2π/λ , where λ is the signal wavelength.The LoS, reflected and scattered delays are computed as τ L (t) = d L (t)/c, τ r (t) = d r (t)/c, and τ s (t) = d s (t)/c where c is the velocity of light in free space, d L (t) is the LoS distance between the transmitter and receiver and d s (t)(d r (t)) is the single bounce scattered (reflected) distance between the transmitter and the receiver via the s-th scatterer (r-th reflector) at time t.The complex amplitudes of the different multipath components in (1a) may be expressed as - where Γ L (t, f ), Γ r (t, f ), and Γ s (t, f ) are distance, time, and frequency-dependent path gains accounting for various signal power attenuation factors and antenna gains for the LoS, reflected, and scattered components respectively and are typically obtained from channel measurements.We translate the classical channel model defined in (1a) to its equivalent ambit representation in (1b) using the following characterizations.κ r (t, y i ; u, y j )(κ s (t, y i ; u, y j )) represents a deterministic spatio-temporal kernel function to capture the received signal amplitude of individual reflected (scattered) multipaths at a particular time instant t.Choosing an appropriate kernel function is essentially dictated by the modeling scenario.The indices ((t, y); (u, y j )) of κ (•) (•) denote the kernel value at the x and y coordinates and their respective temporal variations along the Euclidean plane.As an example, consider a simple wireless propagation scenario in which a MU travels along a straight line.This scenario is modelled using an ambit-process defined on the time-space domain R + × R. Scatterers in the time-space domain located at (t; y) ∈ R + × R can be represented completely in the space domain with equivalent points (x; y) ∈ R 2 .The mapping, t → x is achieved via the transformation x = x 0 + v 0 t, where v 0 denotes the uniform velocity of MU and x 0 is a scalar depending on the initial distance between the BS and the MU.This mapping is adopted in order to obtain further modeling insights using the ambit framework.
L(du, dy j ) is called the Levy base of the ambit channel h(t, y i ; u, y j : (u, y j ) ∈ A (i) (t), i ∈ {r, s}).It captures the large-scale fading statistics associated with the distribution of scatterers and reflectors on the two-dimensional Euclidean plane.L(du, dy j ) is assumed to be a homogeneous Poisson point-process with density λ p .More general distributions of L(du, dy j ) may be used to capture more realistic scatterer dispositions.λ p may be calculated using the expression is the arrival-rate of new reflected (scattered) multipaths, R p,r (R p,s ) is the radius of the reflecting (scattering) zone centered on the instantaneous location of the moving receiver, and v 0 is the receiver velocity.R p,i (i ∈ (r, s)) can equivalently be called the radius of the ambit-sets A r (t) and A s (t) respectively.A physical depiction of the ambit radius is shown in Fig. 1.For simplicity and clarity, only reflected rays are depicted in the figure with scattered rays having similar geometric interpretations.
σ (u, y j ) denotes the stochastic volatility of a reflector (scatterer) located at (v 0 u, y j ).Stochastic volatility captures the small-scale fading along with the localized and random motion of the reflectors (scatterers) about their mean position.In the present work, the reflectors (scatterers) are assumed to be static and thus σ (u, y j ) = χ(u = t) is taken to be a zero-mean unit variance complex Gaussian random variable ∀ (u, y j ) ∈ R + × R. The stochastic volatility is essentially the ambit equivalent of the small-scale fading statistics χ(t) defined in the classical model.
A (i) (t), i ∈ {r, s} denote the time-varying ambit-sets associated with the MU.It represents the collection of reflectors (scatterers) that are in close proximity to the MU and contributes to the major proportion of the multipath power captured at the receiver.Reflectors (scatterers) lying outside the ambit-set of the MU are assumed to have a negligible effect on the received power.The reflective ambit-set is taken to be of the form, A r (t) = (v 0 t, 0) + A, where A is a 2D circular disc with center (0,0) and radius R p,r where R p,r = N r λ p π , N r being the average number of reflectors present within the ambit-set.Thus, A r (t) is simply modeled as a translation of the circular disc A with respect to time t.A similar modeling assumption can be made for the scattering ambit-set A s (t) as well.

IV. AMBIT-PROCESS BASED SISO SUB-THz CHANNEL MODEL
As previously stated, the kernel function essentially captures the received signal amplitude of the reflected multipaths.Accurate characterization of the kernel function is, thus, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. of utmost importance for precise modeling of the sub-THz wireless channel.

A. Characterization of the Kernel Function for Reflected Rays
The instantaneous received amplitude at any time t for a reflected ray in the sub-THz band can be represented by a time-frequency dependent function of the form - where 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 [18] - 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.The absorption coefficient of an outdoor propagation medium is typically dominated by the amount of water vapour molecules present within the 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 equation as follows [18] - 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.A thorough analysis of the absorption coefficient can be found in [18].The vegetative loss in urban areas for frequencies ranging from 110 GHz to 170 GHz can be described using the following empirical model.The average loss of received power per meter depth of vegetation is given by [42] - where p represents the plant area index (PAI), d ρ(s) is the depth of the vegetation in meters through which the reflected (scattered) ray travels and f is the frequency of operation.A,B,C,D are constant parameters obtained through a minimum-mean-square-error (MMSE) fit of the measured data and are given by A=20.4,B=-0.4,C=0.3, and D=0.9.
The plant area index represents the density of the vegetation and is directly proportional to it's height.PAI for different vegetation types can be obtained from [42].Now, consider a characterization of the time-variant sub-THz channel, where the base-station is located at (BS x , BS y ).Then, at any time t, the location of the receiver moving horizontally along the x-axis with velocity v is given by ) where y i is some constant.Similarly, without loss of generality, the coordinates of the r-th reflector located within the ambit-set of the receiver can be taken to be Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
) 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, d(x 1 , y 1 ; is the Euclidean distance between points (x 1 , y 1 ) and (x 2 , y 2 ).Assuming that φ r is the angle of departure and θ r is the corresponding angle of arrival via the r-th path, we have, from elementary geometry - Substituting 7(a) in 7(b), we have - which captures the double directionality of the ambit channel.
The frequency-dependent reflection coefficient defined in (3) accounts for the rough-surface reflection losses at sub-THz frequencies which is predominantly dependent on the material properties, shape and roughness of the reflecting surfaces.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 [43].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 [44] - where θ ri is the angle of incidence on the reflector at any time instance t given by - Under the proposed framework, d r 1 (t) ≃ d L (t) ≫ d r 2 (t), 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 (9) remains reasonably accurate at sub-THz frequencies.η t refers to 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 (9) is the negative sign associated with the Fresnel reflection coefficient which simply indicates a phase shift of π caused by the reflection phenomenon.
The Rayleigh factor of a rough-surface defined in (8) depends on a Gaussian distributed statistical parameter called the rough-surface height standard deviation (σ rs ) described in [43].The Rayleigh factor is characterized by this roughness effect and expressed as [45] - where J 0 (•) is the zero-order Bessel function of the first kind.
To make the model consistent with the framework proposed in (1), we define the kernel function as follows, κ r (t, y i ; u, y j ) = Γ r (t, f )e jkd r (t) δ (t − τ r (t)), where Γ r (t, f ) is the time-frequency dependent function for rough-surface reflection derived in (3) and d r (t) = d r 1 (t) + d r 2 (t) is the total distance travelled by the r-th reflected ray.Plugging in the respective values from equations ( 3)-( 11) in (12), we get the complete spatio-temporal kernel function of a reflected ray as - or equivalently, where ϒ( f ) is as defined in ( 8)- (10) and n is the path-loss exponent.

B. Characterization of the Kernel Function for the Scattered Rays
The rough-surface reflection coefficient defined in (8) only captures the amount of power lost due to scattering in the reflected direction but gives no insight into the strength of the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
scattered power in a particular direction.The diffuse scattered power in a specific direction can be captured through the Beckmann-Kirchhoff scattering coefficient [46] defined below.
Most of the assumptions taken for the reflected rays hold true for the scattered rays as well.The spatio-temporal kernel function for the scattered rays can be similarly written in the following form - or equivalently, where d s (t) = d s 1 (t)+d s 2 (t) following the definition in ( 6)- (7), is the delay incurred by the s-th scattered ray at time t and Φ( f ) is the Beckmann-Kirchhoff rough-surface scattering coefficient.The scattering coefficient has been defined in [46] as - Plugging γ( f ) from ( 9) in (17a) we have - with the following parameters - where l c is the rough-surface correlation length, σ rs is the rough-surface height standard deviation, k is the wavenumber, θ sc 1 ∈ [0, π/2] and θ sc 2 ∈ [−π/2, π/2] are the angles of incidence and scattering respectively corresponding to the s-th scatterer.Additionally, in (17b), we use Taylor's approximation for the exponential term and round-off the factorial sum to its second power.These approximations have been shown to have good agreement with experimental results on diffuse scattering for rough-surfaces at large incidence and scattered angles [46].

C. Characterization of the Ambit-Process Based Spatio-Temporal Autocorrelation Function
The ambit-process based generalized spatio-temporal autocorrelation function (ACF) of the channel coefficients at a given propagation delay τ is shown to be [5] - where w = (w t , w s ) ∈ R + × R are small increments in the x and y directions at time t respectively and A(τ) is the circular ambit-set of all reflectors with delay τ at time t centered at the receiver.Consider a scenario where the reflector r is located at (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 is taken to be deterministic, one can assume that the temporal ACF would only depend on the time-difference w t .Then, the difference in path-lengths of the signals received at positions X and Y are - where θ is the angle made by the r-th reflector with the x-axis.Now, since reflector r is located within a circle of radius R = R p,r around the MU, we can equivalently represent it in terms of its polar coordinates as (v 0 u = r cos θ , y = r sin θ ) where r ∈ [0, R] denotes the radial distance of the reflector from the center of the ambit-set and θ is the corresponding angle of arrival.The general expression for the normalized temporal ACF can then be obtained as - e ( jk[d(v 0 w t ,0;v 0 u,y)−d(0,0;v 0 u,y)]) dtdy = 1 where J 0 (•) is the zero-order Bessel function of the first kind.
A closer look at (19a) reveals that ρ h (w t ) eventually reduces to the temporal autocorrelation in Clarke's 2D isotropic scattering model [13].Furthermore, the ACF of the stationary one-ring model [6] given by r µ µ (τ, 0) = 2σ 2 J 0 (2π f max τ) is equivalent to the ambit ACF ρ h (w t ) = J 0 (kv 0 w t ) derived in (19a) scaled by a power factor of 2σ 2 due to the complex channel gain.The derivation of (19a) is provided in Appendix B. Remark: The above discussion pertains to the derivation of the theoretical temporal ACF of a SISO ambit channel and its relation to Clark'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(t, τ) is the ambit channel impulse response.
The effect of diffused reflection and scattering in the propagation environment is captured in ρ(w t ) through the kernel functions κ (•) (t, y) proposed in (13).
V. SIMULATING THE SUB-THz SISO CHANNEL Next, we develop a convolution-based channel simulation algorithm for an outdoor BS-MU sub-THz SISO 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 phasesnamely the initialization phase and the channel simulation phase.Compute 2D convolution of Y mat and row wise folded X mat .Equate the convoluted sequence to Z i;mat .

A. Initialization
The ambit channel model defined in (1b) 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 continuoustime 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.A schematic of the spatial Euclidean grid has been shown in Fig. 1.

B. Channel Simulation
The algorithm leverages two key properties of an ambit-process to efficiently simulate the channel.Firstly, it decomposes both the reflective and scattered kernel function κ(•) into two components denoted by κ(•) and κ(•) such that - where Γ (•) (t, f ) is as defined in ( 14) and ( 16).(v 0 t, y i ) is the coordinate of the MU and (v 0 u, y (•) ) is the coordinate of the corresponding reflector or scatterer within the ambit-set at time t.We note that owing to the horizontal motion of the MU, y i and y (•) are assumed to be constants and d(v 0 t, y i ; v 0 u, y (•) ) can be approximated by the average distance between the reflector or scatterer and the receiver within the ambit-set.We assume, without loss of generality, the initial location of the MU to be (0,0) (thus y i = 0).This results in model simplicity and tractability without compromising on the statistical characterization of the channel, since, irrespective of the direction of motion of the receiver, we can always define a set of appropriate coordinate systems under which the receiver is only moving along one of the axes while remaining fixed with respect to the other.Such a generalization is possible because we assume a uniform velocity of motion.The rationale behind this assumption is the overall simulation duration T max .When T max is sufficiently small, it is reasonable to assume that the velocity remains constant.The composite kernel function κ(•) can be expressed in terms of κ(•) and κ(•) as [47] - 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 (22) over each sample point of the spatial grid on the 2D Euclidean plane over which the receiver moves [47].The last part of the algorithm computes the received signal amplitude of the LoS components for each time instance using the kernel function κ L (.) which is defined as - where being the LoS distance between the transmitter and the receiver and n being the LoS path loss exponent.In summary, the proposed algorithm separately tracks and calculates the LoS and NLoS received amplitudes corresponding to different time and delay points and does an element-wise addition of the two components to obtain the final 2D channel matrix.

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.Therefore, an approximate complexity of Algorithm 1 may be given by O(MNPD 2 ).

VI. SISO SUB-THz CHANNEL SIMULATION RESULTS
This section discusses 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 MU is capable of moving at uniform velocities of 5 to 500 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 [22].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 δ τ = 0.0092 ns and the maximum delay τ max = 2.36 ns.The rough-surface standard deviation (σ rs ) 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 η t 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 [45].The atmospheric absorption coefficient µ( f ) at f = 140 GHz is 0.4 × 10 −5 / m [18].
The simulated channel with the above parameters gives a mean coherence bandwidth of roughly 13 MHz over a transmission bandwidth of 1.6 GHz and a receiver velocity of 40 kmph.This translates to having roughly 128 subcarriers if OFDM signaling is employed.A typical worst-case scenario sees the same receiver experience frequency flat fading at a reduced coherence bandwidth of 4 MHz thereby increasing the number of required subcarriers to 512.This is evident from Fig. 2(a) where we see the minimum coherence bandwidth may fall to values as low as 4 MHz for certain channel realizations.Such scenarios demand the use of highly directional ultra-massive MIMO arrays employing 512 or more antenna elements and subcarriers.Subsequent discussion regarding this follows below.
Fig. 2(a) shows the steady-state convergence of the coherence bandwidth over different iterations of averaged-out realizations of the ambit-set. 1We see that the mean coherence bandwidth becomes almost flat after approximately 32 iterations, giving an indication of the convergence rate of the proposed algorithm.The above observation holds true for the median coherence bandwidth as well.Contrary to this, the modal, minimal, and maximal coherence bandwidths vary drastically over different iterations, simply re-establishing the fact that they are not ideal measures for convergence analysis or system design.A similar analysis involving the deviation of the coherence bandwidth about its central tendencies is shown in Fig. 2(b).It may be concluded from the plots that  both the mean (median) coherence bandwidth, as well as the mean deviation of the coherence bandwidth about its mean (median), remain relatively constant over increasing iterations of the channel simulation, thereby indicating fast convergence of the proposed simulator.
Figs. 3 and 4 compare the simulated and theoretical autocorrelation functions (ACF) [5] and Dopplers of the ambit channel model with the recently proposed TeraMIMO model [48].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 [49], the simulated and theoretical ambit ACFs slightly deviate from 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 scattering, 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.The diffuse scattering is further exacerbated by the motion of the receiver and the surface-roughness of the reflectors within the ambit-set.It is to be noted 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 channel correlation time.The correlation time is the time instance at which the temporal autocorrelation falls to 1/e or 0.37.Evidently, 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, thus reinforcing the validity of the proposed modeling approach.

VII. SISO SUB-THz CHANNEL MODEL VALIDATION
This section presents a preliminary comparison of the proposed channel model with recent measurements reported in literature [9], [22].The validation parameters used in this section have been drawn from the measurement settings and values obtained in the channel sounding campaign carried out in [9] and [22].

A. Overview of the Validation Scenario
The proposed model is validated against two distinct sets of measurement data obtained from [9] and [22].The static measurement environments considered in [9] and [22] and their translation into our proposed framework are briefly discussed below.Reference [22] considers an urban microcellular environment(UMi) where the transmitter is kept fixed and the receiver moves along a rectangular trajectory of 39 m×12 m (102 m long), where each consecutive and adjacent receiver Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.location is 3 m apart.To model the statistical distribution of scatterers in the NYU campus, where the experiments were conducted [22], we utilize a Google map of the campus area and calculate the transmitter-receiver separation distances along with their respective locations.Brief descriptions of the scattering environment in [22] prompt us to keep the average number of scatterers at a given time t fixed at 10, i.e N s (t) = N r (t) = N s = N r = 10.Since the measurements were explicitly done in an almost static setting, we restrict the arrival rate of new multipaths to a small value of 0.1 (i.e.R p,r = R p,s = 0.1) signifying that a new scatterer would arrive at the measurement scenario only every 10 s.The receiver velocity is restricted to a very small value of 0.0028 m/s to mimic an almost static scenario, and the simulation time T max is taken to be 850 s.Such a setting results in a receiver displacement of 2.3 m over a transmitter-receiver separation distance of 100 m.This compares well against the actual measurement scenario where consecutive receiver locations are 3 m apart.The delay resolution is fixed at δ τ = 1.1 ns assuming a bandwidth of 1 GHz.The maximum delay τ max is kept at τ max = 18 ns.The channel is assumed to be sparse with higher path-losses associated with higher-order reflections.The sampling time interval ∆ is set to 0.83 s while the spatial sampling interval is taken to be Λ = ∆v 0 5/18.
While carrying out the validation, we manually switch between h L (t, τ) (LoS channel) and h N (t, τ) (NLoS channel) depending upon the type of propagation scenario that is being validated for each transmitter-receiver location pair.All remaining parameters in the simulation environment remain unchanged.Since each transmitter-receiver location pair explicitly corresponds to either an LoS or an NLoS scenario, a manual switching between h L (t, τ) and h N (t, τ) for validating the CIR at the corresponding transmitter-receiver location pair yields reasonably acceptable implications as far as validations are concerned.
The relevant assumptions comply fairly accurately with real case scenarios as can be seen from the subsequent validation plots.A good level of agreement is noticed between the measured and simulated channel parameters in terms of the correlation distance of the RMS delay spread, shadowing, the path-loss, and the angular distribution of the received scattered power.A similar methodology has been followed for characterizing the scattering environment and measurement setup discussed in [9].

B. Validation of Path-Loss, Scattering & Shadowing
Fig. 5 compares the measured and simulated path-losses over different receiver location indices as shown in [22].The receiver follows a trajectory identical to the one shown in [22] while the surrounding scatterers are distributed following a Poisson point-process.A similar simulation environment is used in the subsequent validation plots as well.A correlation coefficient of 0.8067 is observed between the simulated and measured path-losses for a purely reflective channel model.The LoS region displays a path-loss bias of 9-10 dB in the absence of scattering.For a purely reflective model, the kernel function is defined as - × e jkd r (t) δ (t − τ r (t)) We modify the reflective channel model to a scatteringreflective one to effectively capture the variation in path-loss and remove some of the discrepancies arising due to the modeling bias.We redefine the composite kernel function as followsκ rs (t, y i ; u, y j ) = κ r (t, y i ; u, y j ) + κ s (t, y i ; u, y j ) × e jkd r (t) δ (t − τ r (t)) where ϒ( f ) and Φ( f ) are correction factors for rough-surface reflection and scattering that have been elaborately characterized and discussed in Section IV.Note that, a purely reflective Fig. 6.Comparison between simulated and measured total received power at different incident angles at 140 GHz [45].
model simply corresponds to the specific case when ϒ( f )=1 and Φ( f ) ≈ 0. The solid curve in Fig. 5(a) shows the simulated path-loss for the modified model.Coincidentally, the path-loss bias has been effectively removed with an improved correlation of 0.8947 between the measured and simulated values.Fig. 5(b) shows the scatter plot of the omnidirectional path loss against a transmitter-receiver separation distance of up to 100 m at 140 GHz over 32 receiver locations in a typical urban environment.The path-loss exponent (PLE) for the LoS component is reported to be almost 2.01 [22] which is apparently very close to the typical free-space path loss value of n=2.Reference [22] reports the PLE for the NLoS scenario in a typical urban environment to be 3.20.It was further observed that the shadow fading in LoS conditions experienced a slight standard deviation of only 2.9 dB while for the NLoS scenario, it increases to 7.1 dB.The scatter plot of Fig. 5(b) is obtained by plugging the above values into the proposed model and the results generated thereof are compared to the measurements reported in [22].The measured values are distributed somewhat randomly from n=2 to n=3.20 while the simulated values are much more concentrated.This deviation may again be attributed to differences in the measurement environment being considered by the proposed model and the actual environment in which the measurements were carried out.
Fig. 6 compares the angular distribution of the total simulated received power to the total measured received power for different incident angles as obtained from [45].We use a ray-based approach and the Beckmann-Kirchoff scattering model to capture the average diffuse scattered power.However, it should be emphasized, that, being a ray-based approach, our proposed model does not take backscattering into consideration and only the scattered power in the reflected direction is accounted for.It has been shown in [45] that backscattered power typically remains 20 dB below the peak measured power and has therefore been excluded from the current modeling approach.The difference between the simulated and measured total power in Fig. 6 stems from the slightly different modeling approaches presented in the two works.The lower sidelobes obtained from the Beckmann-Kirchoff model indicate the absence of backscattering.Additionally, the Beckmann-Kirchoff approximation yields good results if the surface is devoid of sharp edges, spikes, and irregularities.Moreover, multiple scattering and shadowing effects are excluded from the model so that neither depolarization nor scattering at low grazing angles is accurately covered unless the model is suitably modified.The above factors contribute to mismatches between the simulated average received power and the measured total power.
However, it should be emphasized that there is an exact match between the total received power in the specular direction as obtained from both the simulation and measurements.Additionally, it can be seen that both simulated and measured power in the reflected direction is greater at larger incident angles compared to smaller incident angles.An approximate decrease of 5 dB is observed in the received power as the incident angle decreases from 60 • to 10 • .This observation is in line with the results presented in [45].Fig. 7 compares the simulated path-loss and shadowing to the values reported in [9].Reports the PLE for LoS conditions in the measurement environment to be 1.88 while for the NLoS scenario, n increases to 2.57.A close inspection of the simulated shadowing also shows a good deal of agreement with the measured values outlined in [9].For the LoS setting, the modeled shadowing has a mean of 0 dB and a variance of 2 dB corresponding to the measured shadowing which shows a mean of 0.09 dB and a variance of 1.19 dB.For the NLoS case, simulations exhibit a mean of 0 dB and a variance of 15 dB compared to measurements which indicate a mean of 0.04 dB and a variance of 10.3 dB.The measured and simulated shadowing are within 1 dB, and 5 dB of each other for LoS and NLoS scenarios respectively.This can be attributed to the fact that a dominant LoS can be more accurately modeled compared to NLoSs due to the randomness of the scattering environment.

C. Validation of Spatial Correlation of the Delay Spread
One of the critical attributes of any stochastic model should be its spatial consistency.Spatial consistency corresponds to the phenomenon whereby a moving receiver or multiple closely spaced receivers experience similar channel conditions Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.within a local area.Any acceptable channel model should be able to emulate this correlation by effectively generating spatially consistent large and small-scale channel parameters.Fig. 8 plots the spatial correlation coefficients of the simulated and measured RMS delay spread.The correlation coefficients are empirically calculated from the measurement dataset and subsequently fitted using analytical functions yielding minimum mean-square error (MMSE).Two widely-used functions that have been extensively applied in various channel models such as 3GPP and NYUSIM to characterize the spatial correlation are the exponentially decaying function and exponentially decaying sinusoidal function respectively [22].The correlation distance has been defined in literature to be the distance where the correlation coefficient first drops below 1/e(≈ 0.37).
The spatial correlation function of the RMS delay spread of a sub-THz channel can be typically modeled using a double sinusoidal function having an exponential decay of the The solid curve in Fig. 8 corresponds to the simulated MMSE correlation curve while the dashed one depicts the measured correlation.We see that the correlation distance of the simulated RMS delay spread is 12 m while it is 11.8 m for the measured channel response.A correlation of 0.9838 is observed between the fitted curve reported in [22] and the MMSE curve obtained from the proposed model.Additionally, for the measured delay spread, D 1 = 25.5 m and D 2 = 8.9 m [22] where D 1 , D 2 are as defined in (20).Comparing this to D 1 = 73.5 m and D 2 = 9.9 m for the modified model, we see an acceptable agreement between the modeled and measured values.
The exponentially decaying sinusoidal correlation function provides a good fit for both the simulated values and reported measurements, albeit with slightly different parameter values.Additionally, the oscillating nature of the spatial correlation may be due to the symmetric structure of the environment and the rectangular route over which the measurement campaign and simulations take place.A similar comparison of the spatial correlation of the delay spread as reported from the measurements carried out in [9] for different transmitter-receiver locations is shown in Fig. 9.In each of these cases, the measured and simulated correlation distances tabulated in Table I are within 0.2 m of each other.Due to a high spatial granularity (≈ 7 m), small propagation wavelength, high diffused scattering, and the sparse nature of the measurement campaign, a very low correlation distance is observed at all these locations.[22], on the other hand, conducts measurements with a spatial granularity of 3 m resulting in an improved correlation distance of almost 12 m.This is evident y-axis, we can we set w s = 0.This reduces the spatio-temporal ACF to a purely temporal ACF given by ρ h(τ) (w) = A(τ) κ(t + w t , y)κ * (t, y)dtdy A(τ) |κ(t, y)| 2 dtdy (29) Let, initially, the receiver be located at any point (v 0 t, y i ), where y i is some constant.Then, from ( 13) and ( 14), the kernel function of the receiver at any time t is given byκ(t, y) = Γ r (t, f )e jkd r (t) For a small temporal increment of w t , the corresponding kernel function is given byκ(t + w t , y) = Γ r (t + w t , f )e jkd r (t+w t ) (31) Note that Γ r (•) is a distance-dependent channel gain given by × (V r (p, f )) − 1 2 ϒ( f ) = c 4π f (d(BS x , BS y ; v 0 u, y r ) + d(v 0 (t + w t ), y i ; v 0 u, y r )) n/2 × e − 1 2 µ( f )(d(BS x ,BS y ;v 0 u,y r )+d(v 0 (t+w t ),y i ;v 0 u,y r )) The reflector locations in (32) and ( 33) are denoted by (v 0 u, y = y r ) for some constants u and y r respectively.When w t is sufficiently small, the signal amplitudes can be taken to be approximately equal i.e.Γ r (t, f ) ≈ Γ r (t + w t , f ) Plugging ( 30)-( 33) in ( 29) and using the above approximation, we haveρ h(τ) (w) = A(τ) κ(t + w t , y)κ * (t, y)dtdy A(τ) |κ(t, y)| 2 dtdy = Γ r (t + w t , f )Γ r (t, f ) A(τ) e jkd r (t+w t ) e − jkd r (t) dtdy (Γ r (t, f )) 2 A(τ) dtdy ≈ A(τ) e jkd r (t+w t ) e − jkd r (t) dtdy A(τ) dtdy ≈ 1 Leb(A(τ)) A(τ) e jkd r (t+w t ) e − jkd r (t) dtdy (34) A(τ) is the circular ambit-set of all arriving multipaths with a propagation delay of τ and Leb(A(τ))=πR 2 denotes the Lebesgue measure of the set A(τ) (definition is available in [50](Sec.2.6)).Without loss of generality, we assume that at t = 0, the MU is located at the origin (0,0) which implies y i = 0. Since A(τ) is independent of the receiver velocity and underlying propagation mechanisms, we assume it to be of the form, A(τ) = A, where A is a 2D circular disc of radius R centered at (0,0) with R given by R = R p,r = N r λ p π .The ambitset A(τ) may have any regular or irregular geometric shape depending upon the modeling scenario.A circular ambit-set has been assumed to ensure tractability and gain additional insights into the modeling framework.

Fig. 1 .
Fig. 1.A snapshot of the ambit-set based BS-MU channel model over a 2D Euclidean grid.

Fig. 2 .
Fig. 2. Asymptotic analysis and variation of the channel coherence bandwidth over different channel realizations.

Fig. 8 .
Fig. 8. Simulated vs Measured spatial correlation coefficient of the RMS Delay spread of a 142 GHz Channel [22].

TABLE I VALIDATION
OF CORRELATION DISTANCES