Experimental Investigation of Terahertz Scattering: A Study of Non-Gaussianity and Lateral Roughness Influence

The scattering phenomenon caused by rough surfaces has a dominant role in shaping the reflected field at terahertz (THz) frequencies, both in specular and non-specular directions. Most surfaces in nature are randomly rough, and the surface height obeys a certain statistical distribution. A Gaussian probability density function (PDF) for height distribution is often considered, and the correlation length is assumed to be longer than the wavelength. However, a clear understanding of how changing these assumptions affect the angular distribution of the scattered field is still lacking. In the first part of this work we investigate via microscopic measurements the statistical distribution of realistic indoor materials, and its deviation from the assumed normal distribution. After that, the influence of non-Gaussianity on the specular reflection in the low THz region is shown analytically. In the second part, a measurement campaign of diffuse scattering, caused by structured statistically-controlled surfaces, is reported. The correlation length assumption has been proven experimentally and via full-wave simulation to affect the diffuse scattering by rough samples, when the other statistical parameters are kept without changes.


I. INTRODUCTION
Global mobile data traffic is expected to exceed 77 exabytes per month by 2022 [1], which is more than double the data traffic in 2019. In order to fulfill such extreme data rate demands, the carrier frequencies of wireless links need to be expanded into higher ranges beyond 100 GHz. Due to the abundant availability of bandwidth, terahertz (THz) wave technologies are proposed as one of the key fields of wireless communications for 5G and beyond. Although THz applications have already reached maturity in some fields of science, like astrophysics and imaging, THz frequency range is one of the least investigated regions of the electromagnetic spectrum for wireless communications. Very high propagation losses The associate editor coordinating the review of this manuscript and approving it for publication was Santi C. Pavone .
confine communication applications at THz frequencies and restrict them to short-range links in indoor environments [2]. Nevertheless, advances with respect to the generation of extremely high power THz waves [3] give a good indication that communications with THz signals are becoming more and more feasible. This is also reflected in the establishment of a new IEEE standard (IEEE 802.15.3 d standard) [4]. Here, a portion from the low THz band (252 to 325 GHz) was assigned to upcoming wireless communication links.
One significant difference of THz propagation compared to the currently allocated frequencies is the scattering phenomenon that appears due to the roughness of indoor surfaces. Material roughness determines the degree of scattering that occurs, with several parameters contributing to a specific distribution of the scattered rays. The non-specular scattering in the indoor environment causes a major multi-path contribution and thus has a significant impact on propagation characteristics, primarily in the absence of a direct propagation link [5]. Empirical methods have been suggested to estimate the distribution of the scattered field. For example, in [6], the effective roughness (E-R) model was developed. Here, the field distribution in the angular domain is approximated by assuming a Lambertian scattering pattern and developed later to a directive pattern. However, such models characterize roughness by only measuring or estimating one single parameter, the surface heights' standard deviation.
Moreover, given the lack of a closed-form solution for the scattering problem, the effect of scattering is also included in ray-tracing models by employing analytical approximations. The most utilized method is the physical optics (PO), also known as the Kirchhoff approximation (KA). Scattering from surfaces obeying the KA was comprehensively investigated [5], [7]- [9]. In [5], the Beckmann-Kirchhoff (B-K) model, which is an advanced form of the KA, was implemented in these scattering models. The authors showed that the diffuse scattered power has a substantial impact on the non-line-ofsight (NLOS) propagation paths. A comparison between the E-R method and the B-K model for THz propagation has been described in [9].
When applying the Kirchhoff model, it has been assumed in all the previous studies that the surface heights follow a Gaussian distribution. The reason behind this assumption is apparent, and thus significantly reducing the complexity of the scattering problem, as the slope distribution function can be obtained directly. In reality, no surface is normally distributed all in all, and minimal attempt is made to characterize the roughness distribution of real surfaces and their deviation from normality. In an original paper from Beckmann [10], the scattered field from a normally distributed surface was compared to an exponentially-distributed one possessing the same standard deviations of heights. The scattering was significantly different, even for the specular direction, except for a slightly rough surface. The deviation has also been confirmed analytically in the optical domain, even for moderately rough surfaces [11]. To what extent does the surface height distribution of real materials fit a Gaussian distribution? How does the randomness assumption affect scattering behavior? These two questions are addressed in the first part of this work.
Another mandatory assumption for applying the B-K model is that the surface is locally smooth, i.e., the correlation length of the surface is greater than the wavelength of the electromagnetic wave [5]. The surface correlation length is a spatial roughness parameter and generally defined as the distance at which the surface autocorrelation function falls to 1/e of its maximal value. The assumption mentioned above may not be valid in real materials. For example, in [12], the correlation length of the two examined materials was 0.18 mm and 0.29 mm. If we consider a THz wave with a central frequency of 300 GHz, the correlation length in both cases is less than a third of the wavelength. Until today, experimental investigations about scattering, especially in the THz range, from rough surfaces, are very scarce. Reflection scattering measurements in the millimeter-wave range in [13] show that material surfaces exhibiting a roughness standard deviation equal to or greater than 0.3 mm induce a distinct influence on the roughness dependent power reduction of the reflected ray. In [14], [15], measurements of diffuse scattering at THz were reported. Here, sandpapers with different grit sizes were chosen to represent different levels of roughness. The authors in [7] proved experimentally that non-specular NLOS paths play a valuable role in enabling THz wireless links.
It is worth mentioning that since a rough surface is generated from a random process, the field scattered will be a random process. Therefore, measurement of the scattered field from one realization of a random surface is not enough to draw a conclusion about the scattering behavior, and therefore better estimated by an analytical model. In case we aim to investigate the lateral correlation effect, the statistics condition will be fulfilled if the sample size is much larger than the correlation length.
The second part of this work is dedicated to the study of correlation length assumption and its effect on the diffuse scattering profile. This is carried out via bistatic full-wave simulations and scattering measurements on constructed rough materials. The statistical parameters of the samples are monitored and realized by high accuracy photolithographic techniques.

A. ROUGHNESS CHARACTERIZATION
A random rough surface is given by the function h = f (x, y), where h is the surface height at a point (x, y). It is commonly described in statistical terms using two distribution functions, the heights' probability distribution, and the autocorrelation function. The heights' probability density function describes the surface height deviation from a certain mean level and the autocorrelation function describes the x from each other. (1) The shape of the probability density function provides meaningful information about the nature of roughness, and the shape is described mathematically by an infinite number of higher moments, and approximated by the first four moments of the random process h in case of surface characterization. The first moment is the average heights, and it is calculated ash The second central moment is one of the most vital parameters to characterize roughness. The root mean square height is calculated by The third central moment of the heights' distribution determines the skewness, it represents the asymmetric spread of VOLUME 8, 2020 heights' distribution.
The fourth central moment of the heights' distribution determines the kurtosis, and it represents the spikiness of the statistical distribution, In 2002, a working group (TC213/N499) was created aiming to develop international standards for surface characterization. The standard ISO 25178 was then released in 2005 [16]. The first-order moment is the average surface height. Assuming that the random height is an ergodic process, estimates of the moments can be calculated by averaging with respect to the area of a sample A. The second, third and fourth heights' parameters are defined as follows: σ is the only parameter of roughness necessary for calculating the Rayleigh parameter, which is one of the most applied criteria in determining if the surface is considered to be rough in scattering models. The Rayleigh parameter is given as where i is the angle of incidence of the electromagnetic wave and λ is the wavelength.
• Estimated skewness ( S k ) The skewness value within a given limited area is estimated by Surfaces with high peaks and shallow valleys have positive values of skewness and vise versa. Symmetrical distribution functions, including Gaussian ones, have zero skewness.
• Estimated Kurtosis ( K u ) Similarly, the fourth moment is estimated as When K u > 3 the statistical distribution has a distinct peak, and the higher K u is the sharper the peak becomes.
Values of K u < 3 indicate a rather flat distribution.

B. METHODS AND MATERIALS
For standard measurements of roughness, two main methods can be applied, linear or aerial surface roughness measurements. A linear roughness measurement system investigates roughness in two dimensions (x,y) in the form of single lines. In contrast, an areal three dimensional (3D) roughness measurement is conducted over an area of the surface. In this work, the latter method is used in the form of a high accuracy non-contact laser confocal microscopy system. With confocal microscopy, it is possible to investigate the three dimensional surface structures at sub-micron resolutions. A sequence of digital images is acquired for each sample. Then, the output stack of photos is post-processed with the surface metrology software MountainsMap [17]. The reconstructed 3D model is a digital image representing the surface, where each pixel has a value that refers to the surface height. Based on the reconstructed 3D model, the four surface heights parameters, as defined by the ISO 25178 standard and described in section II. A, are then calculated.
In this study, real materials, which are commonly found in an indoor environment, are investigated. Table 1 lists the 15 investigated material, including their measured areas and root mean square heights σ . The maximum patch area that can be measured by the confocal microscopy is 7 mm x 7 mm. For each material, five patches of 7 mm x 7 mm were measured at statistically independent positions, S k and K u parameters are then averaged over the repeated measurements. Only in two cases (L1, W1), the resumed area size is 3 mm x 3 mm.

C. MEASUREMENTS RESULTS
Representative surface maps and histograms of four samples (B1, A2, W2, L2) are obtained from the aforementioned measurements and post-processing procedure. The results are displayed in figure 1. As can be seen from the histograms, some real surfaces do not follow Gaussian heights' distribution. For example, comparing the surface maps and histograms of sample A2 and W2 indicates that two different materials have completely different surface structures, even if their height standard deviations are comparable. The last column of Table 1 presents the calculated standard deviation of the surface heights, ranging between 0.005 mm for the smoothest sample and 0.145 mm for the roughest one.
In figure 2, the calculated skewness and kurtosis values are represented as blue dots for all the samples. The blue bars in both x and y directions represent the error caused by averaging over the repeated measurements. The red dot in this figure has skewness S k = 0 and kurtosis K u = 3, which represents the Gaussian distribution. The green dot represents a Gamma distribution of the fourth degree (κ = 4), which has skewness value of S k = 1 and kurtosis K u = 4.5. Note that the exponential distribution, which is most compared to the Gaussian in older studies, has very extreme values (S k = 2, K u = 9), and is not close to any of the measured real materials, so we did not include it in this figure.
Most of the investigated materials tend to have a skewness value between −0.5 and 0.5, and the kurtosis values are closer to 2 than 3.
Four samples (B1, W5, WP2, WP3) have more extreme values of S k and K u and their deviations from a normal distribution is obvious, WP2 and WP3 could be better represented by a Gamma distribution. In the following section, we will check analytically the effect of the deviation from the assumed Gaussian distribution due to the scattering, within our frequency range of interest.

D. ANALYTICAL DESCRIPTION OF THE SCATTERING
In this section, the general theory of scattering by random rough surfaces, described in [18], is adopted. The surface is assumed homogeneous and continuous, with an average surface height equal to zero. The surface height is a random field. The specular reflection factor depends on k, and is given as in [18] S spec (k ⊥ ) = where I tot is the total integrated intensity of the reflected and scattered fields, and I spec is the intensity of the specular direction. The specular reflection factor has a direct relation to the surface heights' probability density function.

E. THE GAUSSIAN CASE
For a random surface with a surface height that follows the Gaussian distribution, the PDF of the surface height is given By employing in relation (11), the specular reflection factor for the Gaussian case is

F. THE NON-GAUSSIAN CASE
The specular reflection factor can be calculated if f h (h) is known as mentioned in equation (11). A non-Gaussian surface refers to a very general case as any surface which does not follow the normality assumption is non-Gaussian. Here, we will give two examples of non-Gaussian surfaces and calculate the expected specular scattered intensity. The first is the exponential distribution (S k = 2, K u = 9), one of the most compared to the Gaussian in previous works. The second is the Gamma distribution of the fourth degree (S k = 1, K u = 4.5), which fits well to some of the measured samples, as mentioned in the previous section. An exponential distribution has the following form [18] while the Gamma distribution is modeled as in [18] Integrating into equation (11), the specular reflection factor from an exponentially distributed surface is given by The specular reflection factor for Gamma distributed surface is then calculated as By evaluating equations (13), (16) and (17) numerically, the intensity of the scattered field in the specular direction is calculated for the three statistically different surfaces. Figure 3 shows the different curves of specular reflection factor with respect to increasing the heights' standard deviation σ at a normal incident angle with a carrier frequency of 300 GHz. From figure 3, we notice that the assumption of a Gaussian height distribution underestimates the reflected signal in the specular direction for large values of σ . Already for σ = 0.2 mm the specular reflection component in case of the Gamma distribution is by 4 orders of magnitude larger than assuming a Gaussian distribution. Moreover, the exponential distribution shows a similar specular intensity like the Gaussian one for low values of σ , but exhibits a much higher specular intensity for very rough surfaces where σ reaches the wavelength of the incident wave. Nevertheless, if we take a look back at the measured materials, the roughest material has a standard deviation of heights equal to 0.0145 mm. At a 300 GHz carrier frequency, the deviation in the specular intensity, caused by changing the distribution, is small and can be neglected. In other words, the Gaussian distribution is still valid even when the real statistics of the material differ from Gaussian. This is true for low THz frequencies, which attracts our interest in this study. Note that the important factor here is the ratio between σ and the wavelength λ.
Changing the frequency to a higher level has the same effect as increasing σ , and the surface distribution assumption is expected to have a noticeable deviation in the scattered intensity.

III. EXPERIMENTAL MEASUREMENTS ON THz SCATTERING BY ROUGH CONTROLLED SAMPLES
The aim of this section is to investigate the THz wave scattering by 3D printed samples with defined surface statistical parameters. The resulting measurements can then be used to develop a better understanding of how the correlation length assumptions affect the accuracy of THz propagation models.

A. SAMPLE PREPARATION
The numerical generation of random surfaces with specific roughness parameters is well studied in the literature. Here, we adopted a spectral method similar to the one described in [19]. A random rough surface with the desired properties is generated starting with generating an uncorrelated Gaussian random rough surface distribution using a Gaussian random number generator, with an average of zero and different standard deviations σ . To achieve a correlation of surface points, the distribution is then convolved with correlation function having an exact value of the correlation length. This can be implemented numerically by applying a Fast Fourier Transform (FFT). The samples are then exported as standard tessellation language (STL) three dimensional (3D) models. The models are printed following Stereolithography (SLA) 3D printing method. SLA provides a high vertical resolution of 25 microns. The material used is standard Resin. The steps for obtaining rough samples with controlled statistics are depicted in figure 4.
The dimensions of the discussed surfaces are 100 mm x 100 mm x (5 ± σ ) mm (length x width x height). Finally, a layer of pure copper (99.9% copper) is applied to prevent the scattering inside the material, and the copper can be considered as perfectly conducting material at low THz frequencies. The layer thickness of the coating is not known accurately but is significantly thin, so it does not affect the surface statistics. A representation of a numerical 3D generated surface, the STL model and the realized 3D sample are shown in figure 4. Table 2 provides presumed and actual statistical parameters of the printed samples, the S0 sample is a smooth surface having the same size and coated with copper.

B. MEASUREMENT SETUP
The measurement setup, shown in figure 5, was used to calculate the bistatic angular distribution of the THz signal scattered by the rough samples. The Frequency-domain terahertz platform (i.e., THz-FDS), with 40 MHz frequency resolution, is used to generate the THz incident signal. The transmitter is a continuous-wave THz transmitter, consists of photodiode coupled with integrated antenna and a silicon lens.
The THz signal emitted by the transmitter is collimated with a 2'' PTFE lens with a back focal length of 63.2 mm. The lens is placed 100 mm from the center of the sample and the angle of incidence is fixed at θ i = 22 • . The lens assures that the collimated THz beam becomes nearly a plane wave. The THz receiver is made of an antenna-integrated photo-conductor with a silicon lens and is placed over a motorized rotation stage. Moreover, a pinhole with a diameter of circa 1 mm is placed in front of the lens, the pinhole eliminates angular dependence of the radiation pattern on the receiver and rises the spatial resolution. The receive antenna is placed 100 mm away from the rotation center, and by rotating the rotation stage, the receiver antenna rotates on a circular track around the sample, in the range between 60 • and 0 • with a step of 1 • . The observation angle is considered relative to the angle under which the specular reflection is estimated. VOLUME 8, 2020 C. RESULTS AND DISCUSSIONS Figure 6 and 7 present examples of the scattering by rough samples in the form of bistatic radar cross section (RCS) full-wave simulation for a specific frequency (300 GHz) with a normal incident angle. The blue curve in both figures represents the normalized RCS in the case of a smooth perfectly conductive surface, whereas rough surface RCS values are calculated from 10 realizations of each rough sample and then normalized according to the peak of the blue curve, which is the specular reflection. In Figure 6, the rough sample standard deviation is similar to the maximum value of a normal indoor real material, and the correlation length has the value of three times the wavelength l > λ, i.e., the values which are frequently assumed when applying KA. Compared to the smooth surface, a degradation of around 20 dB is observed in the specular reflection, and the strongest side lobe of the diffuse scatter in around 15 dB less than the specular reflection. The diffuse scattering effect is evident in the non specular angles, and it is around (20-30) dB higher compared to the smooth surface reflection in the range −50 • to 50 • around the specular. In figure 7, the rough sample has the same heights' standard deviation, but a shorter correlation length l < λ. In this case, the reflected specular power degradation is 10 dB, with a small amount of diffuse scattering outside the specular direction. It appears that although samples in figure 6 and 7 differ only in the lateral correlation length, and the rest of the statistical parameters are the same, the scattered field distribution is different. The sample with a correlation length shorter than the wavelength acted more like a smooth surface, and the reflection is more specular. The measured and simulated RCS of samples S0, S1, S3, and S4 are depicted in figure 8 for an incident angle of 22 • and frequency range 200 to 400 GHz. In the worst case, at 200 GHz, the ratio of the surface diameter to the wavelength is 66.7. A good agreement is found between the plane wave scattering simulation and the real measurements. The scattering angular distribution of the two samples where the correlation is longer than the wavelength (S2 and S3) shows no noticeable difference. Thus, S3 angular scattering is only shown. The reflection by smooth surface S0, shown in figure 8 (a) and (b) causes a clear specular narrow peak in the complete range of frequency, with a slight decrease of specularly reflected power close to 400 GHz. A comparison between measurements results in figure 8 (c) and (e), and also full-wave simulation in the same figure (d) and (f), indicates that a correlation length shorter than the wavelength makes a surface with fixed roughness height's standard deviation more similar to a smooth reflector. Thus, it increases the specular reflection at a given frequency. The previous effect of the correlation is dominant for slightly rough surfaces. In the case of highly rough materials, the specular reflection is completely destroyed, and the power is scattered randomly all over the angular field. Figure 8 (g) and (h) depict an example of how a very rough sample, S4 here, would destroy the specular reflection.
In order to keep the comparison fair, this study is limited to perfectly conducting samples. No parameters regarding the nature of the surface were included, such as its permittivity and conductivity. However, it would be interesting to calculate the RCS of rough samples considering also these two parameters. The metal coverage of the samples incurs a complete reflection of the waves, without losses inside the material. A further investigation of dielectric samples, with several values of permittivity and conductivity, is planned to follow this work. (a) S0 measurements, (b) S0 simulated,(c) S1 measurements, (d) S1 simulated,(e) S3 measurements, (f) S3 simulated, (g) S4 measurements, (h) S4 simulated.

IV. CONCLUSION AND FUTURE WORK
We report measurements-based investigation to validate two assumptions, which are commonly adopted when calculating the scattered field distribution by Kirchhoff analytical approximation. First, real statistical distributions of rough indoor materials are measured, and their deviations from the assumed Gaussian distribution are calculated. Thereafter, an analytical study of specular scattering with changing the Gaussian assumption is conducted. We conclude that the distribution assumption has a minor effect in the lower THz band, within the limits of real indoor materials. However, a Gaussian height distribution significantly underestimates the measured signal in the specular direction for large values of the surface heights' standard deviation compared to the wavelength. This indicates that the assumption has no effect only in lower THz band, but such an assumption should be avoided in case of higher frequencies. Statistically controlled rough samples are designed and fabricated, and measurements were also performed on the manufactured samples to inspect the effect of alterations in the rough sample correlation length. If the surface correlation length is larger than the wavelength, the diffuse scattering caused by the surface roughness increases drastically. It is observed that when such a condition is not fulfilled, the rough surface causes more reflection, taking the specularity roll-off to higher frequencies. We believe that the obtained results act as an advanced step in developing realistic models for angular distribution of the scattered fields in the THz band.