{'fraction_non_alphanumeric': 0.03172835349270895, 'fraction_numerical': 0.011422207257375224, 'mean_word_length': 4.80339696377574, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Combined Stochastic and Physical Framework for Modeling Indoor 5G Millimeter Wave Propagation

Indoor coverage is a major challenge for 5G millimeter waves (mmWaves). In this paper, we address this problem through a novel theoretical framework that combines stochastic indoor environment modeling with advanced physical propagation simulation. This approach is particularly adapted to investigate indoor-to-indoor 5G mmWave propagation. Its system implementation, so-called iGeoStat, generates parameterized typical environments that account for the indoor spatial variations, then simulates radio propagation based on the physical interaction between electromagnetic waves and material properties. This framework is not dedicated to a particular environment, material, frequency or use case and aims to statistically understand the influence of indoor environment parameters on mmWave propagation properties, especially coverage and path loss. Its implementation raises numerous computational challenges that we solve by formulating an adapted link budget and designing new memory optimization algorithms. The first simulation results for two major 5G applications are validated with measurement data and show the efficiency of iGeoStat to simulate multiple diffusion in realistic environments, within a reasonable amount of time and memory resources. Generated output maps confirm that diffusion has a critical impact on indoor mmWave propagation and that proper physical modeling is of the utmost importance to generate relevant propagation models.


I. INTRODUCTION
Three use case families are defined for 5G mmWave applications. eMBB (enhanced Mobile Broadband) will deliver 10x higher data rates than 4G. mMTC (massive Machine Type Communications) will require high network capacity due to the huge number of connected devices.
URLLC (Ultra-Reliable Low Latency Communications) will support applications with stringent latency and reliability constraints.
Although mmWaves are a promising solution to the spectrum scarcity problem, their propagation properties, namely their very short wavelength, can severely impact transmission. These impairments become even stronger in indoor-to-indoor communication scenarios due to numerous constraints imposed by the environment, leading to a new major challenge: indoor coverage. This is where indoor-specific studies on mmWave propagation modeling become essential for planning indoor 5G networks. This subject is, however, very complex to investigate since indoor environments vary a lot in terms of spatial configuration. Even when considering a particular geometry, its interior characteristics (furniture positions, doors, etc.) and materials (wood, concrete, etc.) are very diverse. Moreover, indoor material surfaces are generally rough; this roughness is comparable to the wavelength of mmWaves and thus, electromagnetic (EM) diffusion cannot be ignored and must be studied.
Existing investigations rely on measurement campaigns or system simulations. Conducting measurements is a heavy process and requires a lot of resources. The few studies [1], [2] based on this method carry out measurements in a specific environment type and setting; collected data are accurate for that particular scenario only and do not account for the propagation diversity inherent to an environment's topology, morphology and material properties. Since it is impractical to conduct measurements in various environment configurations, generated empirical models cannot be used in other indoor mmWave applications, nor can they provide information regarding the influence of indoor characteristics on propagation mechanisms, especially diffusion.
System simulations offer a flexible alternative to predict propagation behavior. Existing indoor mmWave simulators allow users to study propagation in various environment settings. Each of these simulators uses different implementation techniques of propagation mechanisms; however, most of them suffer from major drawbacks when it comes to diffusion which is the most complex mechanism but certainly the one that should be mostly investigated. In [3] for example, the diffusion lobe which physically describes the angular distribution of EM energy is oversimplified and set by the user, choosing between a Lambert or directive model, and whether or not to include backscattering, setting also the back lobe. Some simulators as [4] utilize empirical models, and others as [5] do not consider diffusion due to its implementation challenges. Hence, these simulators cannot predict the actual indoor mmWave behavior, nor can generated models provide conclusions regarding the impact of indoor geometry and materials on mmWave propagation.
In this paper, we present a novel modeling framework for 5G indoor mmWave propagation that combines stochastic indoor environment generation with advanced physical propagation simulation. Its system implementation, so-called iGeoStat, utilizes the theory of stochastic geometry to generate parameterized typical environments that account for the statistical indoor variability, and employs the complex but comprehensive physical model of He [6] to simulate radio propagation based on the physical interaction between EM waves and indoor materials.
Practical implementation of this framework raises highly challenging computational tasks that we solve by formulating an adapted link budget and designing new memory management and optimization algorithms. This makes iGeoStat the first to simulate multiple diffusion in realistic environments, allowing us to study the actual indoor behavior of 5G mmWave propagation.
The proposed framework aims to statistically understand the influence of the environment's geometry and materials' parameters on mmWave propagation properties, especially coverage and path loss. Therefore, generated models are not defined by the geometry or material itself, but rather by the parameters that characterize them. This framework is not dedicated to a particular environment, material, frequency or use case, making it very efficient for indoor 5G mmWave propagation modeling.
The contributions of this paper are summarized as follows: 1) We introduce a new stochastic modeling approach that generates parameterized and realistic 3D environments, taking into account the indoor spatial variability.
2) We present the first implementation of the physical model of He for radio propagation simulation based on the physical interaction between EM waves and indoor materials.
3) We formulate an adapted link budget that enables multiple diffusion simulation and power tracking at any point in the environment. 4) We design advanced algorithms for memory optimization, multiple diffusion management and simulation acceleration.

5)
We generate the first physical-based output maps (power, SINR, coverage, path loss and delay spread) for 5G Fixed Wireless Access at 60 GHz and industry 4.0 at 26 GHz.
The rest of the paper is organized as follows. Section II presents the theory behind the stochastic modeling of indoor environments. Section III presents the physical model of radio propagation simulation. The system implementation of our framework is presented in Section IV, along with the implementation challenges in Section V. In Section VI, we evaluate the performance of iGeoStat for two major 5G applications, and validate our simulation results with existing measurement data. Finally, conclusions and perspectives are drawn in Section VII.

II. STOCHASTIC GEOMETRY
Stochastic geometry is a powerful tool used to generate mathematical models based on spatial probabilities [7]. It provides statistical information about the random configurations we wish to analyze. This tool is used in many disciplines like computer vision [8], wireless [9] and fixed [10] networks, where analytic calculation and simulation can be done with few input parameters.
Stochastic geometry incorporates the theory of random tessellations [11], defined as the random division of space into convex non-overlapping polygonal regions; i.e. partitions. This section sketches the theory behind two random tessellations implemented in iGeoStat to generate typical parameterized indoor environments.

A. Poisson Line Tessellation (PLT)
A PLT is a set of random lines driven by an underlying Poisson Point Process, where each point represents the origin of a line drawn according to an angular distribution in [0, π[. The associated PLT in the plane is a set of random lines delimiting polygonal cells (Fig. 1), characterized by its intensity L A , i.e. the mean total edge length per unit area. Anisotropy can be described by a probability distribution R in the space of directions and is measured via the anisotropy parameter ξ (Eq. 1), where | sin (u, v)| is the area of a parallelogram drawn by unit vectors u and v [12]: The isotropic case corresponds to uniform R and ξ = 2/π. The simulator uses an anisotropy coefficient α derived from ξ and defined between 0 (isotropic) and 1 (anisotropic).

B. STable by ITeration tessellation (STIT)
The STIT tessellation introduced in [13] is indexed with time. The initial window has a random exponentially distributed lifetime t after which it is divided into two new cells by a random line.
These new cells are each attributed independently a lifetime that drives their division by a new random line. The parameter of the time probability distribution is chosen here inversely proportional to the cell perimeter so that larger cells tend to die sooner. The whole process is stopped at an arbitrary time that can be transformed into the mean total edge length L A depending on the anisotropy. The result (Fig. 1) is a STIT formed of polygons that do not depend on the initial window and is parameterized by L A . Anisotropy ξ can be also considered here. We note that the interior of a typical STIT (L A , ξ) cell has the same distribution as a PLT (L A , ξ) one.

C. Morphological features
The morphology of a stationary planar tessellation is described by the mean values of the number of nodes, cells, etc. per unit area. Mean value formulae [12] exist for PLT and STIT topologies (Table I), allowing us to fit the tessellation parameters to real data as in [14].

III. EM PROPAGATION MODELING
This section presents the propagation model implemented in iGeoStat to simulate reflection and diffusion. Reflection occurs when a wave is incident on a smooth surface; the outgoing 'specular' direction is given by Snell-Descartes law and the fraction of reflected power is provided by Fresnel's coefficients. Diffusion, on the other hand, is more complex and occurs when a wave is incident on a rough surface. This mechanism is crucial for indoor mmWave propagation and must not be neglected. Its behavior should not be user-defined but solely based on the physical interaction between EM waves and material properties.

A. Diffusion from rough surfaces
Surface roughness is defined here with respect to the incident wavelength. Typical indoor materials (carpets, furniture, etc.) have rough surfaces with regards to mmWaves. Fig. 2 illustrates the transition from specular reflection to diffusion, induced by increasing surface roughness. Fig. 2: Impact of increasing surface roughness (left to right) on the diffusion lobe (from [15]).
The problem of diffusion from rough surfaces is addressed according to three different approaches. (I) Empirical ones reproduce measured data without any physical interpretation [16]- [19]. (II) Geometrical ones use ray optics assuming that irregularities are much larger than the wavelength [20]- [22]. (III) Physical ones apply EM theory; they are very complex but more general [6], [15]. We note that physical diffusion models are used in computer graphics and video games [23], but have never been implemented to simulate radio propagation.
In the context of (III), we chose to implement the physical model of He [6] since it encompasses all the other models, supports different types of surface roughness and includes all major EM propagation mechanisms: reflection, diffusion, surface diffraction, interference, polarization, masking and shadowing. Although developed for visible light, its physical approach allows it to be used for any frequency.
The geometry of the diffusion problem is illustrated in (Fig. 3). An EM wave propagating at a wavelength λ is incident on a rough interface S separating the air and a medium of complex refractive index n c (λ). (0, x, y, z) are the natural Cartesian coordinates defined by the surface's average plane (z = 0) and the incident plane wave vector k i . We are interested in finding the diffused intensity at any observation point Q in the far field (θ d ∈ [0, π 2 ]; φ d ∈ [0, 2π]) where the diffused wave can be assumed planar with wave vector k d . The illuminated area is L x × L y such that L x , L y λ. The polarization directions ì s and ì p, associated to a planar wave ì k are respectively in the z = 0 and ( ì k, ì z) planes.

B. He's Model
He's model solves the diffusion problem using the Kirchhoff Approximation integral (KA) [15] (which only holds for very large surface curvatures compared to the wavelength), coupled with a random process modeling of the surface irregularities, since they cannot be accurately described. The model assumes that the height distribution z = ζ(x, y) is spatially isotropic and Gaussian of zero mean value and σ 0 standard deviation. The dimensionless parameter σ 0 /λ defines the surface roughness. Height values taken at two points on the surface are separated by an autocorrelation distance τ (Fig. 4). Using this method, He's model provides the statistical average of the diffused intensity I d (Q) over the realizations of the random surface S, through the Bidirectional Reflectance Distribution Function (BRDF) [24], [25].

C. BRDF
The BRDF, denoted ρ, represents the spatial distribution of I d (Q) and is defined as the ratio of total diffused intensity in the direction (θ d , φ d ) to incident intensity on the surface element in the direction θ i .
ρ is the sum of ρ sr (specular), ρ dd (directional diffuse) and ρ ud (uniform diffuse) (Fig. 5): The first two components result from first-order surface reflections on the tangential planes: ρ sr is due to reflection from the mean surface and ρ dd is due to diffraction scattering from surface irregularities. The third component ρ ud results from multiple surface and subsurface reflections.
We note that due to transmission losses in the lower half-space, the integration of ρ over the upper-half space should verify He's model provides the ρ expression for a polarized incident fields (pol = ì s or ì p): S denotes a shadowing and masking function, N 1 is the specular contribution due to reflection from the mean surface, and N 2 is the directional diffuse one due to diffraction scattering from surface irregularities. F s and F p describe cross-polarization effects due to the orientation of the ζ height field and are functions of the Fresnel coefficients for the angle bisector of ì k i and ì k d .
The missing third component ρ ud is generally approximated by a function a(λ), estimated from experiments [26]. In our case, it is deduced approximately from an estimation of the transmitted component that ensures energy conservation (Section V-B).
For non-specular directions the field is always incoherent [15] and it is impossible to reconstruct the polarization vector of the diffused field. Hence the unique way to consider a series of multiple diffusion is to assume non-polarized incident field and use power link budget (Section IV-B). The non-polarized BRDF is defined as the statistical average of ρ pol considering an incident polarization angle uniformly distributed in [0, 2π]. Noting F nonp =< F s + F p >: Detailed analytic expressions of S, N 1 , N 2 , F s,p and F nonp are found in [6]. Their numerical computation is far from being straightforward and was performed as a first step in Mathematica [27] for validation purposes by comparing the results to figures from papers [6], [15], [28], [29].
The input parameters to compute ρ are presented in Table II. The wavelength λ defines the reference length scale. The material's complex refractive index n c (λ) required for the Fresnel coefficient is not usually known a priori; measurement values can be found in [30], [31] for some typical indoor materials, measured at mmWave bands. Fig. 6 illustrates spherical 3D plots of ρ nonp that correspond to parameters (Table II) inspired by He's example of roughened aluminum. These plots allow us to visualize how the incident energy diffuses in space as a function of (θ d , φ d ). The effect of an increase in surface roughness is illustrated by the transition from a near specular reflection to a large diffusion lobe. As roughness increases, the contribution of ρ sr diminishes with respect to ρ ud . The case of a smooth surface (σ 0 = 0) is naturally included in ρ nonp and does not require specific treatment.   (Table II). The incoming direction is materialized by the thin blue pencil.
We note that the incidence angle θ i has a strong impact on the directivity of ρ nonp , tending to specularity as θ i increases. For the sake of simplicity, ρ nonp + a(λ) is hereafter denoted by ρ.

IV. SYSTEM IMPLEMENTATION
The system implementation of our framework, so-called iGeoStat, is written in C++ and consists of three main modules: the first one (Section IV-A) generates parameterized indoor environments based on the random tessellations discussed in Section II, the second one (Section IV-B) simulates radio propagation according to the mechanisms presented in Section III. From this, the third one generates various output maps based on the measurement plane (Section IV-C).
Note that all parameters below are user-defined in the iGeoStat configuration file.

A. Generating parameterized indoor environments
This module uses stochastic geometry and random tessellations (Section II) to generate various typical 3D indoor environments like apartments, commercial centers, open-space offices, etc., using a minimum number of input parameters.
Generating a random tessellation: the simulator starts by defining the boundaries of the area on which to build the environment. The reference system is the orthogonal 2D Cartesian coordinate system (Fig. 7) plus a vertical coordinate h. The 'initial cell' represents the common floor, centered at the origin. Its polygonal shape is defined by the radius of its circumscribed circle and the number of its sides. This cell is then divided into sub-cells according to a PLT or STIT topology, choosing as reference the areaĀ or perimeterP of a typical cell (Table I).
To generate realistic floor plans with rectangular partitions, we set the anisotropy coefficient to The tessellation can be then rotated by a tilt angle with respect to the x-axis ( Fig. 7(c)).   Table III.
All the cells, including the initial one, are then stored in a vector of cells as a circular linked list of edges. Each edge is a 2D coordinates vector of its 'start' and 'end' vertices, where the 'start' vertex of the first edge is the 'end' vertex of the last one.
Drawing a 2D floor plan: the cells have to be modified to represent the layout of the different rooms or sub-spaces in the indoor environment, separated by corridors. The simulator proceeds by iterating through the vector of cells, creating inside each one of them a polygon of equal number of parallel sides, distant by a half corridor width w c (minus sampling technique [11]).
Segment lines delimiting the cells are then removed to visualize the spacing between adjacent polygons, now separated by the width of a corridor.
In order to represent doors or entrances of the indoor rooms, a single opening is added to a randomly selected edge of center m and length l. The user chooses among two options 'pcent' or 'rand'. With 'pcent' and parameter w 0 (%), a door of width w 0 l centered at m is drawn. With 'rand', a door is drawn with a uniformly random w 0 l around m. These modifications transform a random tessellation from an arrangement of polygons to a typical 2D floor plan (Fig. 8).     A custom mode for generating indoor environments is also implemented in the simulator.
In contrary to the random mode, it allows users to generate a personalized environment based on given 3D coordinates data. The above mentioned fitting methods [14] can determine the tessellation parameters that would best fit a given deterministic environment.

B. Radio propagation simulation
This module simulates radio propagation in the generated 3D indoor environment, where diffusion is simulated according to He's physical model (Section III-B).
Propagation methodology: to consider diffusion, it is fundamental to keep track of the (θ d , φ d ) angular distribution of ρ (Section III-C) at each impact point, which is strongly dependent on θ i . Our implementation is based on the discretization of diffusion directions, with steps ∆θ d and is attributed a solid angle and diffused power, enabling straightforward power tracking through an adapted link budget. Diffusion is thus monitored at each impact in a physical way; this allows us to reduce the set of (θ d , φ d ) directions by adapting discretization steps or adding thresholds on ρ and power values (Section V-C).
Antenna parameters: the antenna is placed in the reference 3D Cartesian coordinate system of the environment at (x a , y a , h a ). Its position can be user-defined with respect to the geometry of a deterministic environment, or randomly set in a stochastic one. This approach aims to derive statistics on various propagation indicators related to the global parameters of the environment.
By assumption, the antenna emits uniformly a total power P 0 at wavelength λ in a solid angle Ω a centered in the spherical direction (θ a 0 in azimuth, φ a 0 in elevation). Using dedicated 13 discretization steps ∆θ 0 and ∆φ 0 , Ω a is divided into N 0 elementary solid angles d 2 Ω 0 , each carrying in its own direction (θ 0 , φ 0 ) a power P 1 = P 0 d 2 Ω 0 Ω a to be included in the link budget.
BRDF: the BRDF from He's model is implemented separately from the main core of iGeoStat using the mathematical library PARI [32]. To avoid computing ρ for the same input parameters over and over again, numerical simulations of ρ are executed for a whole set of input parameters and stored in databases using SQLite [33]. Hence, at each impact point, the incoming direction and surface normal define the incidence angle ψ j (Fig. 10) and the ρ value of each diffusion Each BRDF database corresponds to computing ρ for a single material (σ 0 , τ, n c ), wavelength centered at S 1 (x 1 , y 1 , z 1 ). P 1 (W) denotes the power received by d 2 A 1 from some source in the far field direction ψ 1 with respect to its normal (Fig. 10). The power diffused in the solid angle where ρ 1 (ψ 1 , θ 1 , φ 1 ) is the non-polarized BRDF computed for the surface parameters, the corresponding angles, and wavelength λ. The direction (θ 1 , φ 1 ) from S 1 defines a point S 2 (x 2 , y 2 , z 2 ) (assumed to exist) at distance r 1 on another planar surface, as well as angle ψ 2 with respect to the normal at S 2 . The relationship between the illuminated surface d 2 A 2 and the solid angle d 2 Ω 1 allows us to write P 2 as a function of surface element and distance as P 2 = P 1 ρ 1 d 2 A 2 cos ψ 2 /r 2 1 . The power received on an infinitesimal surface after j diffusion is then given by: Gas attenuation is not relevant here, but transmission and absorption losses in material are taken into account by ρ.  Measurement points are also pointers that contain the same information as a regular impact point, only that these information are evaluated and stored at the measurement plane and used to generate various output maps such as the impacts map, received power, Signal to Interferenceplus-Noise Ratio (SINR), antenna coverage and delay spread. colors, but the output is very precise as each grid element covers a small environment area.

V. IMPLEMENTATION CHALLENGES
We present in this section the main challenges of implementing the data structure that manages diffusion, as well as the challenges of computing, storing and retrieving BRDF values and preventing the explosion of diffusion directions.

A. Tracking impact points
As discussed in Section IV-B, when diffusion occurs, multiple propagation directions are launched after every impact. Each direction corresponds to a potential source of diffusion or new impact point. The main challenge is to implement a data structure that keeps track of the impact points in order to accurately compute their corresponding propagation information (received power, traveled distance, etc.) when undergoing multiple diffusion.
Choosing the optimal data structure: several structures were investigated: the linear ones (arrays, linked lists, etc.) and the hierarchical ones (trees, graphs, etc.). However, the dependency of an impact point on its previous one induces a parent −→ child relationship, which mostly correspond to a tree structure. BFS is a corecursive level-order algorithm, i.e. it uses self-produced data bit-by-bit as they become available to produce further bits of data. It is implemented using a queue; starting from the root, tree nodes are explored by visiting horizontally all the current level adjacent nodes before moving to the next level. Referring to Fig. 12  Various external open-source libraries were investigated like PARI, SageMath [34], Boost [35], etc. We selected PARI since it is written in C, contains all the required functions, and supports multi-precision computations.

Memory management and optimization: PARI initializes a stack in the memory where
BRDF values are computed. The objective behind our implementation is to minimize memory usage and allow numerical simulations regardless of the system's performance. However, without any optimization, computing ρ for a single incident angle consumes a huge amount of memory and eventually crashes when simulating for all θ i . This major issue was tackled by implementing a customized algorithm for manual memory management and optimization (Algorithm 2).
Starting from ∆θ d = ∆φ d = 1 • , their optimal values are automatically selected by decreasing them until the variation of (Eq. 6) is less than an error threshold. Computing ρ for a single incident angle, material and wavelength, takes 21 seconds for a 1 • discretization and up to 30 minutes for a 1/32 • discretization.

VI. FIRST SIMULATIONS PERFORMANCE RESULTS
We present in this section the simulation results for two major 5G indoor applications: Fixed Wireless Access (FWA) in an apartment at 60 GHz and industry 4.0 (i4.0) in a factory at 26 GHz (Fig. 13). Four scenarios are considered ( Table V)     Using measurements from [30] and [31], three BRDFs (Table VI) were generated for the following angles: θ i ∈ [0, π 2 ], θ d ∈ [0, π 2 ] and φ d ∈ [0, 2π] with ∆θ d = ∆φ d = 1 • . This corresponds to 2916000 ρ values per BRDF, computed and stored in database tables using the advanced hardware configuration in 40 minutes with 1215 inserts/second, and is done only once before launching the simulations in iGeoStat. The surface roughness parameters σ 0 and τ are chosen such that the diffusion profile corresponds to a moderately rough surface (Fig. 6(b)).
To test the full potential of our implementation, no thresholds are imposed on ρ nor on power, and all the propagation information are computed and stored for each impact.

A. FWA -reflection vs. single diffusion
In this use case, we first impose a reflection-only simulation to study whether this mechanism can illustrate or not the actual propagation behavior of 5G mmWaves in an indoor environment.
We then introduce single diffusion and compare its contribution on simulating indoor mmWave propagation to the reflection-only scenario.
Indoor environment and materials: using the custom mode in iGeoStat, we generate an apartment with a total area of 150 m 2 where all the walls go up to the ceiling of height h c = 3 m ( Fig. 13(a)). All the obstacles (walls, floor, and ceiling) are covered in "Plasterboard (A)" [30], configured with a moderately rough surface profile (σ 0 /λ = 0.57).  (Table VI).
Antenna and measurement plane: An antenna is placed on the main corridor's left wall at h a = 2.8 m (Fig. 13(a)). It emits P 0 = 20 mW in a large angular aperture Ω a = π sr, Simulation output: the power map of the reflection scenario ( Fig. 14(a)) shows multiple coverage holes in the environment. For convenience, received power values are expressed in dBm, but do not have any physical meaning since such low values cannot be measured. Filtering at -200 dBm ( Fig. 14(b)) shows a 20% environment coverage only. On the other hand, using the same antenna configuration, the power map of single diffusion (Fig. 14(c)) shows a fully covered environment, with typical received power varying between -208 and -17 dBm.
These power maps highlight the crucial contribution of diffusion on indoor mmWave propagation and its advantage over reflection when simulating the interaction between EM waves and materials, especially if properly modeled and implemented. This shows that a classical reflectiononly scenario cannot reflect the actual indoor behavior of 5G mmWaves even with a huge number of reflections. Hence, existing classical ray tracing techniques that mainly employ reflection are not adapted to study indoor 5G mmWave propagation. Additional output maps can be generated in iGeoStat for the single diffusion scenario, based on the propagation information stored by the impacts. Fig. 15(a) is a delay spread map that shows the average arrival time difference of the various multipath trajectories in each grid element.
Values vary between 0.5 ns close to the antenna and 44 ns in its far region. Fig. 15(b) is a path loss variability map that allows to visualize how power decreases in the environment, showing an analogy to Fig. 14(c). Values vary between 29 and 215 dB as we move away from the antenna.

B. i4.0 -single vs double diffusion
Previous scenarios showed that diffusion must be taken into account when simulating indoor mmWave propagation. However, some cases might exist where single diffusion is insufficient to predict mmWaves behavior. Therefore, we first test this mechanism in a different environment geometry and antenna configuration. We then introduce double diffusion and compare its contribution on simulating indoor mmWave propagation to the single diffusion scenario.
Indoor environment and materials: using the custom mode in iGeoStat, we generate a factory with a total area of 1000 m 2 ( Fig. 13(b)). The ceiling is at h c = 5 m, however the walls' height is fixed at h wa = 3 m, allowing radio waves to propagate in the whole environment.
Antenna and measurement plane: two antennas are placed on opposite walls of the environment at h a = 4.8 m (Fig. 13(b)). Both antennas emit P 0 = 200 mW in a narrow angular aperture Introducing double diffusion has a major effect on coverage as shown in the output maps of both antennas ( Fig. 16(b)). They illustrate a fully covered environment by each antenna separately with a very slight increase in received power that vary between -190 and 18 dBm.
These results confirm that in some cases, simulating single diffusion is not sufficient to study indoor mmWave propagation; multiple diffusion simulation would be required to properly illustrate the indoor behavior of 5G mmWaves. These results also show that the main power is concentrated in the first diffusion level, which gives an idea regarding the adequate level of multiple diffusion simulation for a given scenario.

C. Results validation
As a first results validation, we generate the LOS and NLOS path loss PL(dB) curves with respect to the transmitter-receiver distance r(m) for both use cases and compare them to the ones obtained from existing measurements at 60 GHz [36] and 26 GHz [37], [38]. We note that these measurements are conducted in different environments and antenna configurations from ours. Hence, we are more interested in comparing the trend of the path loss curves rather than their exact values, as more simulations are needed to characterize path loss for a given scenario.
60 GHz use case: the power map in Fig. 14(c) are used to generate the path loss curves by fitting 10000 positions based on linear regression, without any filtering. For the LOS scenario ( Fig. 18(a)), the plotted curve PL = 3.2 r + 30 is based on 1648 positions (red dots). The NLOS scenario ( Fig. 18(b)) has a higher and steeper path loss curve PL = 5.5 r + 41 and is based on 8352 positions (blue dots). The resulting path loss curves from measurements in [36] are based on fitting 21 combined positions for LOS and NLOS.
26 GHz use case: the power maps in Fig. 16 For the LOS scenario (Fig. 18(c)), the plotted curve PL = 1.2 r + 55 is based on 6272 positions (green dots). The NLOS scenario (Fig. 18(d)) has a higher path loss curve PL = 0.5 r + 80 with a more gentle slope and is based on 13728 positions (purple dots). The resulting path loss curves from measurement in [37] and [38] are based on 10 and 264 positions respectively for LOS only as they do not consider NLOS.  [37] and [38] (Fig. 18(c)) confirms that diffusion is strongly impacted by the environment settings which leads to various indoor mmWave propagation behavior. Hence, one measurement campaign is also insufficient to characterize path loss for a given scenario.

D. Performance results
These simulation results are very promising; they highlight the ability of our system implementation to simulate multiple diffusion in a reasonable amount of resources and to generate various output maps that illustrate the actual indoor behavior of 5G mmWave propagation.
We note that when simulating any propagation mechanism (reflection, single, double diffusion, etc.), the total RAM usage is equal to the amount used by a single antenna trajectory, independently of the number of antennas and trajectories N 0 . This is due to the efficiency of Algorithm 1 where the same memory resources are recycled and used over and over again. This allows us to anticipate the amount of required RAM prior to the simulation.
The performance results are presented in Table VII. Although FWA scenarios are simulated on a mid-range hardware configuration with N 0 = 16200 trajectories, simulation times are very reasonable and RAM usage is extremely low considering the huge amount of registered impacts, containing all the propagation information. i4.0 scenarios were simulated on a more advanced hardware configuration with N 0 = 16 trajectories; for the single diffusion case, simulation time is very fast and memory usage is extremely low and is as expected, equal to the FWA scenario, despite a huge difference in number of impacts. As stated above, this is due to the efficiency of Algorithm 1. For the double diffusion scenario, simulation time and RAM usage are both relatively high compared to the previous scenarios. As expected, this mechanism is extremely complex due to the explosion of diffusion directions and consequently the number of impacts. Section V-C, like filtering diffusion directions whose power is below thermal noise, reduces the number of impacts. These techniques will also allow us to go beyond double diffusion if needed in a reasonable amount of resources.
Generating such output maps as in iGeoStat would have been impossible with classical ray tracing techniques that mainly simulate reflection. Including diffusion would create many complicated challenges like power evaluation from multiple rays and selecting/verifying the appropriate number of rays that ensures the representation of the whole diffusion lobe. This would correspond to computing the BRDF at each impact, which requires huge computational resources. In our system implementation, the BRDF is computed only once prior to a simulation, and is stored in a database which can be used for other simulations.

VII. CONCLUSION AND FUTURE WORK
This paper introduced a novel modeling framework for indoor 5G mmWave propagation, combining stochastic environment modeling with physical propagation simulation. Its system implementation, iGeoStat, generates parameterized typical environments that account for the statistical indoor variability, then simulates radio propagation based on the physical interaction between EM waves and indoor materials. Challenging computational tasks were solved by formulating an adapted link budget and designing new memory management and optimization algorithms, making iGeoStat the first to simulate multiple diffusion in realistic environments.
First simulations were carried for FWA at 60 GHz and i4.0 at 26 GHz to investigate the contribution of reflection, single and double diffusion on indoor mmWave propagation simulation, where diffusion is the most complex but certainly the most essential mechanism. Illustrated results confirmed that a reflection-only simulation cannot describe indoor mmWaves behavior. They also showed that received power is mostly concentrated in single diffusion; however, simulating this mechanism is not always sufficient to study indoor mmWave coverage, and may require double diffusion. Such conclusions cannot be obtained with classical ray tracing techniques.
Simulation results were successfully validated by comparing the generated LOS and NLOS path loss curves to the ones obtained from measurements at 60 and 26 GHz. The performance results showed the efficiency of iGeoStat to simulate multiple diffusion and generate various output maps (power, SINR, coverage, path loss and delay spread) in a reasonable amount of time and memory. Our implementation enables us to determine the impact of a given multiple diffusion level and to predict the required time and RAM for any propagation mechanism.
Optimization techniques presented in this paper can further reduce simulation time and memory usage, allowing us to go beyond double diffusion if needed. iGeoStat's flexibility will enable us to implement additional modules (indoor furniture, antenna placement optimization, throughput analysis, etc.) to solve indoor planning and dimensioning 5G mmWaves challenges. Extensive statistical studies on simulations results will be carried for various sets of parameterized stochastic environments (geometry and materials) to propose the first analytical formulae of indoor propagation models that are not dedicated to a particular environment, frequency or use case.