The Scattering Channel Model for Terahertz Wireless Communications

Terahertz (THz) communications systems have been considered as one of the enabling technologies for beyond 5G mobile networks, with the ability to offer higher capacity and lower latency compared to the current systems. In this paper, we develop a new THz channel model that captures rough surface scattering while considering attenuation due to molecular absorption of the atmosphere. The main advantage of the proposed model is that it conserves the energy of the channel impulse response, regardless of the number of launched rays in the ray-tracing simulator, or the size of the tiles for the scattering surface. Using the proposed model, a simple propagation scenario for the indoor environment has been analyzed, and the improvement over other models discussed in detail.


I. INTRODUCTION
T HE EXPONENTIAL growth in data demand by endusers and the emergence of new applications in areas like the Internet of Everything (IoE) require a substantial increase in wireless network capacity [1]. However, the limited available bandwidth in the mm-Wave and lower bands makes it difficult to support near-future tera-bit per second demands in beyond 5G communication systems [2]. Within this framework, the terahertz (THz) frequency band (0.1-10 THz) is envisaged as the best potential wireless spectrum, in addition to existing spectra, to meet the upcoming bandwidth demands by the users [3], [4].
Although the THz technology has the potential to solve the capacity limitations of the current systems, there are many challenges about using this technology [3]. One of the major challenges for communications at this frequency band is the extremely high propagation loss due to high spreading loss and high atmospheric absorption, which can exceed 100 dB for a short 10 m link [5]. To address this issue, different solutions including adaptive design, reflectarrays, ultra massive multiple input multiple output (UM-MIMO), and hypersurfaces have been examined by [2] and a joint scenario is proposed which can solve path-loss and fading problems in none-line-of-sight (NLoS) communications.
In addition to high path loss, THz propagation exhibits other characteristics. While object surfaces can be treated as smooth at low frequencies, the same assumption is not valid for THz frequencies. At these frequencies, the dimensions of building material surfaces are comparable to wavelength (about 1 mm), therefore the scattering is more evident than in low-frequency bands. On the one hand, scattering can reduce the power received by the specular reflection component and make the path loss even worse for THz-bands. On the other hand, due to diffuse scattering, the wireless channel modeling in terms of multipath fading is rather complicated in comparison to the lower frequency scenarios [6]. To address this challenge, a number of research efforts in the literature modeled the THz wireless channel both empirically and theoretically.
The importance of measurements to model the THz channel behavior has not been neglected in the literature. In [7], [8], [9], [10], [11], [12] measurements are performed to characterize and model the unique behaviors of the THz channel. In [7] the ray-tracing technique in addition to measurements is leveraged to investigate the cluster behavior of the channel. The multipath channel data of a large-sized indoor environment at 28 GHz and 140 GHz is analyzed in [8] to investigate path-loss, delay spread, and angular spread. In [9] a line-of-sight (LoS) measurement is conducted for the ranges up to 5.5 m in 140-220 GHz frequencies to calculate path loss exponent and fading statistics, while in [10] the authors take a further step to model an outdoor channel with a longer range up to 35 m. To model rough surface scattering in THz frequencies, [11] studies the measured characteristics of building materials, along with modified reflection coefficients for vertical and horizontal polarization, considering scattering of the surfaces. In [12] a 2 × 2 MIMO configuration is employed in THz channel measurement and modeling at 300 GHz. Along with these measurement efforts, stochastic modeling of the THz channel is studied in [13], [14], [15], [16]. In [13] different stochastic parameters of a channel including the number of propagation paths, delay distribution, delay-pathloss correlation, and pathloss-angle correlation are predicted by the ray-tracing approach for three different types of transmitter and LoS/NLoS scenarios, although no information is provided about channel scattering characteristics. In [14] a hybrid channel model is presented, which includes combined spatial, temporal, and frequency domains. The proposed channel is able to model channel frequency dispersion, employing antenna arrays in the transmitter or receiver part, and includes THz scattering in the ray-tracing simulator. Based on this ray-tracing simulator, the stochastic characteristics of the channel parameters including amplitude, polarization, times of arrival, phases, frequency dispersion, and angles of arrival/departure are discussed and relative statistical distributions are proposed. In [15] a statistical model for the THz scattering channel is presented and applied for the MIMO system, and the channel matrix is presented in the spatial and angular domain. Finally, a discussion has been made for channel capacity and spectral efficiency in [15]. Moreover in [16], an analytical method based on the integral equation method (IEM) is presented for computing the scattering pattern of metallic scattering surface and comparison with measurement data. The impact of surface parameters, such as rms roughness height and correlation length, is also studied on the bi-static diffuse scattering pattern.
One of the most vital parts of channel modeling in THz frequencies is the modeling of rough-surface scattering. The scattered power estimated by directive scattering (DS) model [17] and measurement data of a 142 GHz channel for a short-range 1.5 m scenario is presented in [1], which shows satisfactory results. Although the specular reflection term in the DS model is a function of surface roughness, the diffuse scattering term does not contain the rough surface parameters (for example rms roughness height, correlation length, and refractive index) in the formulation. While in the original DS model, the diffuse scattering is modeled by a single lobe pattern, authors in [18] used a modified DS including an extra back-scatter lobe which creates more reliable results compared to measurement results. Furthermore, a comprehensive comparison between the DS model and the radar cross section (RCS) model is done for multiple sub-THz frequencies and different materials. Comparison of the DS model with full-wave simulations using FEKO software is also reviewed in [19], and the difference due to micro-structures is modeled using an additive random electric field component in the DS model's formulation.
In a different approach, [20], [21], [22], [23] used Beckmann-Kirchhoff (B-K) formulations to model THz channel scattering characteristics. While the B-K model is verified as an efficient approach to describe scattering parameters of randomly-distributed rough surfaces, including coherent (specular reflection) and non-coherent (diffuse scattering) components, there are some issues that need to be addressed in the THz channel modeling. Firstly, the original B-K formulation [24] only applies to random rough surfaces illuminated by uniform plane waves with unitary incident electric field amplitude. To address this issue, it is proposed that dividing the surface into a smaller hypothetical grid of small tiles, and summing the scattered power of all tiles, is an effective approach [23]. However, it needs to be carefully applied to the channel model, because the scattered power in the original B-K formulations is dependent on the surface area. Secondly, the model is valid when the dimensions of the surface are large enough compared to the wavelength and surface correlation length [24]. Equally important, the model assumes one polarization (horizontal or vertical) for the incident field, while in realistic channel modeling, we deal with combined polarization. Lastly, the transmitter-tosurface segment of the channel is not in the model and needs to be considered for a transmitter-to-receiver channel. While different research efforts tried to apply the B-K formula to model the THz multipath channel, still there is a gap to achieve physically justified channel impulse response (CIR).
In [20] and [21], the authors account the loss due to molecular absorption in the THz frequency band and present a model for the CIR by separating the LoS, diffracted, and scattered components. While the proposed model is comprehensive, the rough surface considered is small to accommodate the plane wave incidence assumption. Direct extension of the model to a larger surface divided into small tiles would lead to the tile size selection issue. Selecting a smaller tile size increases the received power from the tile due to the dependence of the scattering coefficient of each tile to its area, leading to a larger total received power that is not physical. A similar problem for the specular reflection component is reported in [23] and the authors proposed a different treatment with the specular tile.
In [22] the scattering model based on B-K formulation is combined with the ray-tracing simulator, and the results are used to examine the effect of the surface roughness on the THz channel behavior in terms of the spectral efficiency of the channel. In this paper, although the specular direction is treated separately, non-specular directions include a coherent term in the scattered power's equation. Moreover, the complete channel transfer function between the transmitter and receiver is not discussed.
The other parameter that should be considered in THz channel modeling is the effect of antenna directivity on channel characteristics such as path loss, root mean square (rms) delay spread and fading which are discussed thoroughly in [25].
Also, it is worth mentioning that we should be careful in selecting a model for rough surface scattering because a non-Gaussian random rough surface can have a completely different scattering behavior than a Gaussian surface [26], especially for the surfaces with high rms roughness heights.
In this paper, we focus on the multi-ray channel modeling of the THz wireless channel by taking into account the surface roughness as a Gaussian random variable with zero mean. We propose a physically meaningful use of the Beckmann-Kirchhoff scattering model. Similar to the previous works, we divide a reflection surface into small virtual tiles to account for the multipath scattering components generated by the entire surface roughness. The main advantages of our model are summarized as follows: 1) The diffuse scattering components are independent of the selection of the tile size for rough surfaces.
2) The specular reflection component is independent of tile size.
3) The total propagation paths between transmitter and receiver have been modeled. 4) Both horizontal and vertical polarizations are considered in received power and CIR. 5) The concept of effective scattering area is proposed. These features enable us to apply appropriate tile sizes for different objects and frequencies in channel simulations. The reason of this advantage in our proposed CIR is the independence of the total power of specular and diffuse components to virtual tile size selection of the scattering surface, which is not achieved completely in previous works. Given a scattering surface, regardless of the selection of the virtual tile size, in this paper the total power of specular reflection and diffuse scattering components in CIR as a result of the scattering surface remains unchanged, which matches the behavior of the physical world. Furthermore, the concept of the effective scattering area presented in this paper is a practical tool for reducing the number of accounted paths and simulation time for THz scattering channels. As a result of the abovementioned features, the resulting CIR gives a quantitative description of how indoor materials and geometry impact the channel.
The remainder of this paper is organized as follows. In Section II, the derivation of the mathematical expressions for the received power and CIR for the THz wireless channel is presented based on the Beckmann-Kirchhoff approximation for random rough surface scattering modeling. In Section III, different simulated environments are considered and various simulation results on the CIR, received power, the impact of tile size, and the impact of surface materials are presented and discussed. Finally, the conclusion is presented in Section IV.

II. MODELING OF THE SCATTERING MULTIPATH CHANNEL
In order to model a multipath CIR including scattering components of all rough surfaces in an environment, the first step is to divide each rough surface into small rectangular tiles. A reflection surface may have a large unevenly illuminated area by the incident wave from a transmitter, and hence it is hypothetically considered to be formed by smaller tiles, each of which is approximated by a plane wave incidence and one of its scattered rays reaches the receiver. The second step is to compute the geometrical parameters of wave incidence for each tile to model the multipath nature of the THz channel due to rough surface scattering. As depicted in Fig. 5, these parameters for the specific tile with index i include the angle of departure (θ tx i , ϕ tx i ) at the transmitter antenna, the distance of the transmitter to the tile (r t i ), the angle of incidence (θ 1,i ), the elevation and azimuth angles of scattering (θ 2,i , θ 3,i ), the distance of the tile to the receiver (r r i ), and the angle of arrival (θ rx i , ϕ rx i ) at the receiving antenna. The number of tiles is N, and hence 1 ≤ i ≤ N.

A. SCATTERING MODEL OF ONE ROUGH SURFACE
Denote P t as the total transmitted power at the transmitter. Due to the unequal distancing of the tiles to the transmitter, each tile is illuminated by a different field intensity, and hence we need to compute the electromagnetic field intensity at each tile separately. As a starting point, the Poynting vector which is a measure of the carried electromagnetic (EM) power is given by [27] where E is the electric field and H is the magnetic field vectors, neglecting the time-variant term of e jωt . In the farfield region of the antenna, these field vectors are in phase and make an orthogonal set with the direction of wave (k i ), and the ratio of electric field amplitude to magnetic field amplitude is equal to the intrinsic impedance of the free space, Z 0 = 120π . The direction of the Poynting vector iŝ k i and its time-average value is given by Because the antenna radiation field expression depends on the antenna type, we need to consider a general case including the transmitter and receiver antenna patterns. In this paper, we use the power radiation pattern of antennas in derivations. To get the wave intensity on the surface, we can write power density at the position of tile i as [27] wherek i is a unitary vector in the direction of transmitter to the i th tile and has a simple mathematical relation with the angle of departure at the Tx antenna, G t is the transmitter antenna gain, F t (θ Tx i , ϕ Tx i ) is the normalized power radiation pattern of the antenna in the direction of the i th tile, r t i is the distance between the transmitter and the i th tile, and k a is the frequency-dependent medium absorption coefficient due to molecular absorption in the THz frequency range [5].
Using (2), the magnitude of the incident electric field at the location of tile i can be obtained by where S in i is the magnitude of the incident Poynting vector given by (3), at the location of tile i. The incident fields at each of the tiles, in addition to having different amplitudes, also have different directions or polarizations. If we consider this polarization vector asû E , which is determined by the geometry and the type of transmitter antenna, it can be decomposed into vertical, i.e., transverse electric (TE) or s-polarized, and horizontal, i.e., transverse magnetic (TM) or p-polarized components. The direction of these two components in each tile is determined bŷ wheren is the normal vector of the tile. We define the ratio of the vertical component to the total electric field as where E in i = |E in i | is the amplitude of the incident electric field to tile i, and E in i,v is vertical or TE-polarized component of the incident electric field to tile i.
Rough surface scattering is derived using Helmholtz integral in [24], and formulated as the Beckmann-Kirchhoff model. Given that the formulation in [24] is based on having a unitary incident electric field (|E in i | = 1), the scattered power density of tile i with the non-unitary incident electric field at the receiver is written as where is the reflected field of a smooth perfect electric conductor (PEC) tile with the same size and location as the original tile which is illuminated by a plane wave having a unitary electric field and is equal to [24] where l x and l y are the dimensions of the tile, r r i is the distance between the i th tile and the receiver and k is the wave number which is equal to 2π/λ. The third term, ρρ * i is the average power scattering coefficient of a normally distributed random dielectric rough surface and is equal to where R Fres i is the Fresnel field reflection coefficient of the smooth surface, which is different for the vertical and horizontal polarizations and given by [24] R Fres Here, n the is complex refractive index of the scattering surface, which can be a frequency-dependent function of the intrinsic impedance of the scattering object [28]. In (9), the term ρρ * ∞,i is the average power reflection coefficient of a PEC tile with the normally-distributed rough surface at the same dimension of the original tile of i. Based on Beckmann-Kirchhoff [24], ρρ * ∞,i is given by where l corr is the correlation length of the rough surface, defined as the length that the autocorrelation function of the random rough surface, which is a Gaussian function, drops to 1/e [24]. Other parameters are given by [24] and In (12a), σ h is the standard deviation of the surface height function, or rms roughness, which is the main factor that determines the scattering behaviors of the surface. On the other hand, the correlation length of the rough surface shows how sharp the roughness is, i.e., a higher correlation length corresponds to a smoother surface and smaller values for correlation length show sharper variation in the surface height function. We would like to make a few remarks on (11). The first component e −g i ρ 2 0i is the specular component, which should only exist for the specular reflection direction (θ 1,i = θ 2,i , θ 3,i = 0). For the other scattering angles, this component will be zero in the physical world. The second component is called diffuse scattering term which is due to the random roughness of the surface. The extent of the surface roughness determines the relative value of specular and diffuse components. Surfaces with higher rms roughness heights have higher values of the diffuse part compared to the specular part.
Decomposing the incident electric field into v and h polarizations and using the definition for the κ i , the incident power density for these polarizations are In order to derive the received power, we know that the received power at Rx is given by [27] where G r and F r are the gain and normalized power pattern of the receiver antenna. The gain is fixed, however, the value of the normalized pattern is different for each tile. Combining (3), (6), (7), (8) and substituting the derived expression for S in i into (14), the received power scattered from tile i after a simple mathematical manipulation is obtained as The first impression of this equation is that the received power is dependent on the inverse squared product distance in the other literature formulations, such as CIR of [20]. In the next section, we will show that this new formulation is highly consistent with the conservation of the scattered power. One important thing about this formulation is adding the power received by the different tiles. As mentioned earlier, the Beckmann-Kirchhoff formulation is for a single definite rough surface, illuminated by a plane wave. However, in our problem, which is THz modeling of the indoor environment, different parts of the scattering objects, such as the walls and ceilings, are illuminated by nonuniform amplitudes of the electric field. Thus, the scattering surface needs to be divided into smaller tiles. If the dimension of the tile size (l x l y ) is small enough compared to the Tx distance (r t i ), the impinging wave to this tile can be considered as a plane wave.
As seen in (16) and (11), the diffuse scattered power of each tile is dependent on the area of the tile (l x l y ). Adding up diffuse scattered power of all tiles, given in (17), gives the total received diffuse scattering power that is independent of the chosen tile dimension, written as where |R Fres is the combined Fresnel power reflection coefficient of vertical and horizontal polarizations for the i th tile, and ρρ * ∞,i is given in (11) without the ρ 2 0i term. This is because the specular direction component is calculated separately as follows. Knowing that the field reflection coefficient for the Gaussian random rough PEC surface with rms height of σ h is e −2k 2 σ 2 h cos 2 θ sp [29], where θ sp is the specular angle of incidence, the specular power component received at the receiver is given by  where |R Fres is combined Fresnel power reflection coefficient of vertical and horizontal polarization for specular reflection. In (18) the θ 1 , θ 2 and θ 3 , as illustrated in Fig. 1, are incident elevation angle, scattering elevation angle, and difference of azimuth angles for the incident and scattering directions, respectively. δ(θ 1 − θ 2 )δ(θ 3 ) term is 1 when the incident angle is equal to the scattering angle (θ 1 = θ 2 = θ sp , θ 3 = 0). In other words, this term only exists in the specular direction. The other delta function δ pq , where p and q are the scattered and incident field polarization, shows that the depolarization effect is neglected for the specular direction that is consistent with the literature [24]. Due to the fact that the specular reflection component is calculated using (18), the diffuse scattering component for each tile given by (17), should not contain the coherent component. Finally, the total scattered power received in the Rx is a summation of the (17) and (18).

B. CHANNEL IMPULSE RESPONSE IN THE PRESENCE OF MULTIPLE REFLECTION OBJECTS
Knowing that the CIR is the electric field response of the channel, the received electric field amplitude is related to its power density as The phase information of the received electric field is not known in the literature. Given the random nature of the rough surface, we assume the phase distribution of the diffuse scattering component is uniform over [0, 2π ] [24]. The channel frequency response for the environment with N obj scattering objects is formulated in the frequency domain as (20). In (20), A is a constant and is equal to 2Z 0 P t G t × 4π λ 2 . It is worth mentioning that to be consistent with the literature we can consider A equal to √ P t G t G r , but with this selection, the amplitude of CIR components will be equal to the root square of power components. In (20), d is the LoS distance, τ LoS = d/c is the LoS path delay, r t sp,l + r r sp,l is the total specular path distance between the transmitter and receiver through the specular reflection of the l th scattering object, τ sp,l is the specular path delay for the l th object, r t i,l and r r i,l are the distances of the transmitter and receiver to the i th tile of the l th object, and τ i,l sc = (r t i,l + r r i,l )/c is the time delay for the scattering path containing the i th tile of the l th object. To consider antenna radiation patterns in the CIR, the pattern factors are included for all paths between transmitter and receiver. This is seen as F LoS t , F LoS r for LoS path, F sp,l t , F sp,l r for specular path of l th scattering object, and finally F l,i t , F l,i r for i th tile of l th scattering object. In (20), shown at the bottom of the page, i,l is random phase variable with uniform distribution over [0, 2π ] for the i th path of the l th object, and N l is the number of tiles assumed for the l th scattering object. The term e j∠R Fres sp,l is the phase factor due to vector summation of the reflected vertical and horizontal electric fields and is a mathematical function of the R Fres sp,v and R Fres sp,h angles and vertical and horizontal fields amplitudes. Based on (20), the time-domain multi-ray CIR including delay for frequency f m is written as (21), shown at the bottom of the page.
In (21) it is assumed that the phase change due to the complex Fresnel reflection coefficient has a negligible effect on the total path delay of the specular and scattering components, because when the phase is translated to time delay as ∠R Fres sp,l /2π f , it can be considered negligible in the THz frequency range.

C. SELECTION OF TILE SIZE
Tile size is a parameter to be chosen in the channel modeling. In general, for valid and reasonable tile sizes, the total power of the resulting CIR should not change much with different tile sizes, as the scattering phenomenon is only dependent on the physical rough surface. The selection of tile size has multi-facet effects on the CIR modeling. The value of the tile size determines the time resolution of the multipath components, i.e., to have a higher time resolution the tile size should be smaller. On the other side, the Kirchhoff scattering approximation and the derived CIR (eq. (20)) have validity ranges in terms of tile sizes. Based on [24] the basic criteria for Kirchhoff approximation is min l x , l y l corr .
If (22) is met, the tile size will be much larger than the wavelength, which according to [24] means that the edge effect has no significant impact on the validity of Beckmann-Kirchhoff scattering formulations. On the other hand, our assumption in deriving (20) is that the incident wave on each tile is a uniform plane wave, i.e., the electric field amplitude and phase on each tile are uniform. The main requirement of this condition is that each tile should be located in the far-field region of the transmitter/receiver antennas, which means that [27] where D Tx and D Rx are the dimensions of the transmitter and receiver antennas, respectively. As shown in Fig. 2, we can derive the binomial approximation for the R, which is the distance of the tile edge to the transmitter/receiver, as follows h A c  where L t is the maximum dimension of the tile, and the r is the minimum distance from the center of the tile to the transmitter/receiver location. In (24) if we can neglect the third term, the wave impinging on the tile can be considered a uniform plane wave. To achieve this, the third term should be smaller than λ/16, which gives 22.5 degrees of phase change for an incident wave on the tile. So the condition of the tile size will be Combining (23) and (25), uniform plane wave criteria for the tile can be summarized as follows As it is seen in (26), the upper limit for the tile size is dependent on the environment size and transmitter/receiver positions, because r t i and r r i are determined by environment size. Considering (22) and (26), a rule of thumb tile size for an example normal-sized indoor environment (min(r t i , r r i ) = 5 m) at 300 GHz is 10l corr ≤ l x , l y ≤ 5 cm. However, according to (26) in higher THz frequency bands smaller tiles should be chosen to avoid uniform plane wave approximation error on its surface. Using this criterion for simulation purposes we select two different sets of tile sizes: (1) small tile (l x = l y = 2.5 cm) which satisfies the lower bound for all materials, and (2) large tile (l x = l y = 5 cm). Later in the paper, we will use these tile sizes to study the effect of tile size on CIR.

A. SCENARIO
To evaluate the performance of the proposed formula in (20), we examine a wireless propagation environment consisting of multiple scattering objects. The surface of each scattering object is divided into N = N x × N y rectangular tiles, where N x and N y are the numbers of the horizontal and vertical tiles on the rough surface. Fig. 3 shows the configuration of the transmitter, receiver, and scattering objects in the simulation environment. The solid lines represent the specular reflections and the dashed lines are used for diffuse scattering. The transmitter is located at x = 5 m, y = 5 m, and the  receiver can move along the line which has an angle of π/3 with the horizontal axis. For the first part of the simulation, just the main scattering objects, i.e., walls, will be considered. Later, the secondary scattering objects are included in the simulation. After dividing each surface into small tiles, the received power is calculated using the ray-tracing technique and path gain values given by (20). It is noteworthy to mention that in the simulations the number and size of the tiles will be chosen according to the scattering object's distance to the transmitter/receiver and constructing materials. So for the secondary scattering objects located closer to the transmitter/receiver, we will use smaller tile sizes.
One of the important issues about the modeling of THz wireless channel in an indoor environment is the selection of the scattering parameters, i.e., the correlation length and rms roughness height, for the building materials. Our simulation is based on the values reported by [11] for 300 GHz, which are shown in Table 1. Plaster 3 is an exemplary plaster having moderate roughness compared to Plaster 1 and Plaster 2 and we will use it as an exemplary moderate rough surface. Using these values for the correlation length satisfies the Kirchhoff approximation criteria [24] l corr λ in the THz frequency range, especially in higher frequency bands, such as the 700 GHz band.

B. RESULTS
In the first step, before checking the results of (21) on the CIR, we add some discussions about the B-K's original model about the relative specular and diffuse scattering values of a single tile as shown in Fig. 4, which is made of Plaster 3 with moderate rms roughness height value. The scattering coefficient of a plaster tile with σ h = 0.1 mm and l corr = 1.5 mm is shown in  The incident plane wave in this case has an incidence angle of θ in = π/4, although the other incidence angles result in similar observations. The solid curve is the total scattering coefficient given by (11), while the dashed lines are the diffuse scattering part, excluding the specular term e −g ρ 0 . At the first glance, it is seen that the selection of a larger tile size results in a narrower beam in the specular direction, and the non-specular directions will have a smaller effect on the received power. However, the scattered power density given in (7) has a hidden (l x l y ) 2 which makes the diffuse scattering component unchanged for a specific scattering surface area. In other words, the diffuse scattered power of a 2l x × 2l y tile is equal to the sum of the diffuse scattered power of four l x × l y tiles. On the other hand, the peak which represents the specular reflection's coefficient is identical for different tile sizes, which makes the specularly reflected power density different for different tile choices given the same total scattering surface area, which is not physical. To avoid this effect, in Section II we separate the specular reflection term from (11) to make the final model physically justifiable. Also, the 3D scattering coefficient patterns (11) of a tile with size l0l corr and 50l corr and rms roughness height of σ h =0.3 mm are shown in Fig. 6, indicating that a tile with larger size has narrower specular component. Here, we have selected a high value for σ h for better visualization of the diffuse component. It should be noted that tiles are only a computational tool to model the  multipath components of the surface, and therefore its size should not affect the total received power. The only thing that can affect the received power is the physical characteristics of the rough surface including rms roughness, correlation length, and refractive index.
To have a comparison between the B-K scattering model and previous measurement data, in Fig. 7 we present B-K  (17) and (18) in this paper vs simulated and measured power reported in [30].
scattering results with the data reported in [16] for a synthetically manufactured metallic surface with σ h = 0.77 mm, and l corr = 19 mm. It can be seen that the width of the scattering pattern is similar, although the B-K model can not model the micro-structure scattering of the surface. Also at scattering angles that are away from the specular reflection, the B-K model's predictions are significantly different from the measured data, which is already stated in [24]. This behavior in diffuse scattering also can be seen in the original DS model, where the scattering pattern lacks micro fluctuations. Another thing to notice in Fig. 7 is the specular peak at B-K results, which can be because of the inaccurate selection of the surface size in our simulations, or the low resolution of angular measurements in [16].
To ensure the results of our model, we tested the formula (17) and (18) against the received powers presented in [30] for different rough woods. The receiver and transmitter are located at an identical point that has a perpendicular distance equal to 7.2 m from the scattering object. The frequency is 300 GHz and the gain of receiver/transmitter horn antennas is 25.45 dBi. To use (17) we divided the surface into 51 × 51 square tiles with dimensions of 5 cm × 5 cm. The results are presented in Table 2, which verifies that our formulations are on par with [30] simulation and measurement data.
Another aspect is to investigate the effective area of a scattering surface, which contributes to most of the diffuse scattered power. We define this effective scattering area, A e as a potion of the total surface area that the total diffused power is 90 percent of the actual diffused power. Using a simple scenario shown in Fig. 8, considering isotropic antennas in the transmitter and receiver, and using exemplary plaster parameters (σ h = 0.1 mm and l corr = 1.5 mm), A e consists of 31 × 31 tiles, which is equal to 31 2 (0.05) 2 ≈ 2.4 m 2 . It is clear that rough surfaces with higher rms roughness or smaller correlation length have a larger effective scattering area. This area is illustrated in Fig. 8(a) as the area inside the dashed line. It should be noted that the center point in the figure is the specular reflection point. Knowing this area, especially in the simulation of complex environments with multiple rough surfaces, we only need to take into account the tiles inside the effective area, because the contribution of the other tiles in the received scattered power is less than 10 percent. In our indoor environment with walls over 150 m 2 , the effective scattering area is under one percent of the total surface area, so the simulation can be over 100 times faster. Simulation results for different positions of the transmitter and receiver verified that the effective scattering area includes 200-400 tiles. It should be noted that according to Fig. 8(b), for the longer distances between the transmitter/receiver and the rough surface, the effective scattering area is larger. If we denote A e1 for the wall distance of d 1 , the effective scattering area for the wall distance d 2 is Thus the number of tiles needed for the calculation of the scattered power is (d 2 /d 1 ) 2 times higher than the base scenario in Fig. 8(a). Based on this observation, for the scenario shown in Fig. 3, adaptive tiling should be taken into account for different walls.
In the simulation procedure, the main investigation is to compare the analytical formulas (20) and (21) proposed in this paper with what is available in the literature. Since the channel model in [20] is concerned with only small rough surfaces as single tiles, we apply it to the rough surfaces characterized by multiple tiles and name it the modified [20] channel model. For simplicity, only one scattering object is considered in the simulation. This scattering object is an exemplary plaster surface, with an rms roughness σ h = 0.1 mm, correlation length of l corr = 1.5 mm and refracting index of n = 1.92 − j0.059. The scattering object is divided into N x N y tiles, each of which is like that in Fig. 4. The received power at the receiver is composed of 2 + N x N y components. The first two components are LoS and specular reflection components and the remaining N x N y are diffuse scattering components. The center frequency for the simulation is 0.3 THz and the bandwidth considered is 10 GHz. To consider the frequency-dependent response of the channel, N fft = 256 different frequencies are analyzed. The attenuation of the wave due to molecular absorption is calculated separately at each frequency and the time-domain results for the CIR are obtained using inverse fast Fourier transform (IFFT). Fig. 9 shows the CIR obtained by the formulations of this paper and that in [20] for the size of a single tile equal to l x = l y = 5 cm, and N x = N y = 51. The transmitter and receiver antennas are assumed isotropic antennas because with these antennas the multipath components due to surface scattering can be maximally seen. However, in reality, we probably need to use the directional antennas to compensate the high path-loss associated with the THz band [31]. The LoS component which is the first impulse in Fig. 9a is almost identical for the two formulas, while the NLoS components, based on (20) with considering A = √ P t G t G r , is slightly higher than the modified [20]. The reason for this difference is due to consideration of E ref , so that the magnitude of our CIR for non-specular or diffuse scattering components are (kl x l y cos θ 1,i /π ) × (r t i + r r i )/r t i r r i ≈ 1.4 times bigger than the [20] components. In Fig. 9(b) the comparison between the two formulations is extended to different materials, including Plaster 1, Plaster 3, and wallpaper. As expected, Plaster 1 which is smoother than other materials results in higher specular reflection and lower diffuse scattering components. Also, as shown in Fig. 9b, the difference between the proposed CIR and modified model based on [20] exists for all materials.
To study the impact of changing the tile size on the overall CIR including four scattering walls, simulation is performed using the half dimensions for the tiles, i.e., l x = l y = 2.5 cm. It should be noted that to consider almost the same surface area over the scattering object, the number of tiles in the second setting is doubled to N x = N y = 101. Note that we keep the number of tiles on odd numbers to have a central tile as a specular tile. The measure of interest is the total energy for NLoS components which can be computed using where h[i] is the sampled values of the CIR for NLoS components. Table 3 presents the total energy of the NLoS components using two different tile sizes. While the difference of the energy for the two settings in the proposed CIR is only 12 percent, the CIR of modified [20] shows above 12 times bigger energy value for the smaller tile size setting, which means the formula in [20] cannot be directly extended to multiple tile models. This result indicates that the CIR proposed in this paper has good consistency with the physical world, where we expect to have more or less constant received power (or CIR energy) with changing the tile size. Indeed, the choice of the tile size should only change the amplitude of multipath components, not the total received power in Rx. The smaller the tile size, the larger number of closely spaced multipath components in time with smaller magnitude. Furthermore, the CIR of the environment with four scattering objects for different surface materials, i.e., Plaster 1, Plaster 2, and Wallpaper in Table 1, shows similar results to  Table 3. In Table 4 the total energy and corresponding rms delay spread of the THz channel in 0.3 THz is shown for these materials. Two different tile size is used for the channel modeling: small tile (l x = l y = 2.5 cm) and large tile (l x = l y = 5 cm). For the large tile simulation, the number of tiles is set N x = N y = 51 for wallpaper and N x = N y = 71 for plasters. To keep the scattering area unchanged, for the small tile case, the numbers are doubled. As seen in Table 4,  TABLE 3. Total energy and rms delay spread calculated for two different tile size, using our presented formulation and the one for modified [20] for exemplary plaster with σh = 0.1 mm and lcorr = 1.5 mm. for all of the assumed materials, our model has comparable energy for both small and large tiles scenarios. The energy of modified [20] is much higher in small tile simulations.
According to Table 4, the environment with a Plaster 2 surface has the highest CIR energy or accordingly scattered power. On the other hand, the Wallpaper, due to a smaller refractive index, has the lowest scattered energy. The rms delay spread is also the highest for Plaster 2, which means smaller coherence bandwidth. On the other hand, although Wallpaper has a higher rms roughness height than Plaster 1, it has slightly lower energy and rms delay spread in comparison to Plaster 1. In terms of coherence bandwidth, the Wallpaper has the highest quantity, and in the LoS transmission scenario can have better characteristics. While in the NLoS scenarios, Plaster 1 because of having higher specular reflection components, can make a better physical channel for THz transmission.
To have a better view of this discussion, the CIR for an environment with Wallpaper surfaces is shown in Fig. 10 It is seen that in contrast to the modified [20] model, the proposed model shows comparable CIR for both small and large tiles and total energy of CIR for NLoS components just changes about 10 percent.
The CIR of the environment, comprising four main reflecting surfaces, and four secondary scattering objects, is shown in Fig. 11. The CIR components are labeled with the numbers, and the last component which has a delay equal to about 20 nanoseconds has not been shown in the figure. The first component is the line of sight path gain, and components 2, 3, 4, and 7 correspond to secondary scattering objects with Plaster 1 rough surfaces. The other components, i.e., 5, 6, and 8, are scattered responses of the main Wallpaper walls. Because of having a smaller distance to the transmitter/receiver, the tile sizes for the secondary objects are selected as half of the main scattering objects, so the corresponding components in CIR have narrower widths. Also, because of the smaller rms roughness height of Plaster 1, the secondary objects have stronger specular reflections, so we see higher spikes for 2, 3, 4, and 7. The rms delay spread of this channel is equal to 8.841 ns which is about 2.48 ns greater than the rms delay spread of the environment without secondary objects. The higher rms delay spread means that the coherence bandwidth for the environment, including secondary objects, is about 44 MHz smaller than that of an empty environment. As mentioned earlier, the CIR given by (21) is capable of operating in tile dimensions given by (26). Therefore, we can apply it to different environment dimensions, with various scattering objects and configurations. As a case in point, the environment shown in Fig. 12(a) is tested using the proposed CIR, and the associated CIR is shown in Fig. 12(b). This environment is an exemplary 6 m× 8 m × 3 m conference room, which has Wallpaper walls, a rough ceiling with σ h = 0.15 mm and l corr = 1.7 mm, a smooth glass (with a refractive index of n = 1.45) meeting table in the middle of the room, and three secondary scattering objects (S 1 , S 2 , and S 3 ) made of Plaster 1. As with the previous setup, the selected tile size for the walls (5 cm) is twice the tile size for the secondary scattering objects, which is selected equal to 2.5 cm. In this case, the reflection from the glass table appears close to LoS component, and the other eight scattered components follow after this strong reflection. This environment's resulting rms delay spread is equal to 4.333 ns, which corresponds to the channel's coherence bandwidth equal to B = 1/σ τ = 230 MHz.
Another important parameter that we can calculate for the presented THz channel is the angular spread as defined in [32]. Similar to the rms delay spread concept, the angular spread is the second central moment of the angular power profile. The angular spread can be calculated for a defined reference plane. In our case, this reference plane is the azimuth plane (θ = 0). The angular power profile associated with the environment of Fig. 12(a), depicted in Fig. 13(a), which has an angular spread equal to σ AoA,φ = 37 degrees. We also examine the angular spread dependence on distance. Fig. 13(b) shows the value of the calculated angular spread for the different distances between the transmitter and receiver. The position of the transmitter is fixed as shown in Fig. 12(a), and the receiver is moved along the dashed path from d = 0.5 m to d = 5 m. As can be seen, the angular spread is minimum for d = 0.5 m which has the highest LoS reception value, and increases for higher distances. For lower distances, the effect of scatterers is more evident and so the angular spread is higher.

IV. CONCLUSION
In this paper, we have proposed a complete model including rough surface scattering for the THz wireless channel. This model uses the Beckmann-Kirchhoff scattering model for non-line-of-sight transmission and incorporates two major components for each scattering object: a single specular reflection due to the average contribution of the whole surface and multiple diffuse scattering components due to the random roughness of the surface. The main advantage of our model is that the THz CIR of an environment scales in magnitude and time with the tile size selection for the rough surface and maintains the total power almost constant regardless of the tile size. It also models the entire propagation path between transmitter and receiver, which includes antennas pattern and polarization. A definition for effective scattering area has been developed, as an effective rule-of-thumb tool in computing diffuse scattering components of each surface. In addition to rough surface scattering, wave attenuation by molecular absorption is also included in the THz frequency range which has a significant effect on LoS and NLoS path loss. Based on the proposed model, implementing the ray-tracing approach to model a wireless propagation environment is appropriate because the selection of different tile sizes for the scattering surfaces has a negligible effect on the total power received in the receiver. A valid range of tile size choices has also been given.