Radiative Transfer: Application to Indoor Propagation

1. Introduction

In recent years, with the increasing interest in indoor wireless communications systems, the development of appropriate tools for modeling the propagation within an indoor environment is becoming of utmost importance. A versatile technique for studying propagation in such a complex scenario is ray tracing [1-2]. Through this approach, a number of paths, stemming from the transmitter, are traced along their way to the receiver, accounting for reflection over the obstacles within the scenario. Other mechanisms of interaction between the wave and the environment, such as diffraction, can be accommodated in ray tracing procedures by appropriate generalization of the basic theory [2].

While this method is purely deterministic, in actual environments with many randomly placed scatterers of size comparable to the wavelength (Mie scattering), statistical characterization of the multipath channel [3-6] may be the only viable approach in order to have an accurate model of the propagation [7-8]. Statistical modeling built on iteration of ray tracing results suffers of its inner computational intensity. As for the radio channel design what it is of interest are the fluctuations about a mean value of the received power, one can use a simpler and efficient method to take into account variations due to randomly placed obstacles in the propagation environment. Radiative transfer theory seems to be appropriate, as it deals with the wave propagation within a random medium characterized by randomly placed scatterers. Based on a phenomenological description of the transfer of energy, the basic equation (referred to as radiative transfer equation, RTE) simply states the conservation of energy in terms of the specific intensity $$I(r,\widehat{s})$$ , i. e. the power per unit area and per unit solid angle propagating along $$\widehat{s}$$ , and which is a function of position r in the random medium [9]. The development of the theory is heuristic and does not entail any information about the phase of the wave. Chandrasekhar first studied the RTE within the context of astrophysics [10]. Later, the same formulation was employed successfully for many other applications, such as propagation modeling in the atmosphere or in forested environments, heat transfer through insulating materials, neutron scattering, and power absorption within biological tissues [11-14].

In this chapter, the author reports the RTE results for the evaluation of the power fluctuations in an indoor environment described as a homogeneous medium filled with scatterers arbitrarily placed. The radiative transfer outcome is compared with the ray tracing predictions to assess its limits of applicability.

For the sake of simplicity and without introducing inessential complications in our analysis, the investigation is limited to a propagation environment that can be modeled as a layered parallel plane medium. This modeling is a first approximation of an indoor environment including people, benches or, for instance, a row of chairs in an auditorium, whose location is not fixed. The random medium layer accounts for the average condition out of many possible spatial configurations, given or assumed the number density of scatterers of a succession of layers where each layer is modeled as random medium containing randomly placed scatterers. Moreover, we consider infinite-length circular cylinder as scatterers so as to simplify the solution of the RTE. It should be clear that while the ray tracing approach can be in principle used for any geometry and provide information about the phase of the wave, RT is in practice only applicable to simple geometries and can only yield information about the second order statistics of the wave. The scatterer’s number density is chosen so as interference and interaction between scatterers could be neglected. RT reliability under this condition has been thoroughly investigated [15-17]; here, it is shown how it can be a useful and simple tool for indoor propagation analysis regarding the spatial correlation as defined in Section 5.

The chapter is organized as follows. A review of the radiative transfer is presented in Section 2, together with the definition of main quantities. Section 3 and 4 are devoted to the description of the RTE and of the numerical techniques for its solution. Iterative procedure and its limits of applicability are discussed. Section 5 reports the numerical results in two study cases in a 2D geometry with focus on the impact of the system parameters on the specific intensity and spatial correlation. The results of the comparison between the radiative transfer results and the predictions obtained through ray tracing are also reported. Section 6 is devoted to comments and conclusions.

2. Radiative Transfer Theory: physical background

Two basic theories have been developed in order to approach the study of wave propagation within a random medium characterized by randomly placed scatterers. The first is the analytical theory, where taking into account the scattering and absorption characteristic of the particles solves the Maxwell equations. This approach is mathematically and physically rigorous since in principle the effects of the mechanisms involved in multiple scattering, diffraction and interference can be appropriately modeled. However, in practice, various approximations have to be made in order to obtain feasible solutions (see [9] or [10] for an overview). Recent developments [18] adopt the random medium as paradigm to describe the propagation channel. Stochastic Green’s functions are computed to obtain the channel transfer matrix T in MIMO applications. On the other hand, the transport theory is based on a phenomenological description of the transfer of energy. The basic equation simply states the conservation of energy expressed in terms of the specific intensity and it is equivalent to the Boltzmann equation used in kinetic theory of gases [9]. The development of the theory is heuristic and it does not entail any information about the phase of the wave.

In this chapter, the basic concepts and the quantities of interest when dealing with radiative transfer theory are reviewed, while the reader is referred to literature for a deeper insight. As the problem is defined in a two-dimensional domain, the classic theoretical formulation is adapted to this framework. Accordingly, the chosen reference system is in cylindrical coordinates, as shown in Figure 1. The analysis deals with monochromatic signals with frequency f and the phasor notation is used.

Figure 1.

Coordinate system.

To study the propagation of a wave in presence of randomly distributed particles, the main results related to scattering and absorption of a single particle in vacuum are reviewed. Let us consider an elliptically polarized incident plane wave E _i (ρ) in the direction t _i described by the azimuth angle ϕ _i:

$${\widehat{t}}_i=cos{\varphi }_i\widehat{x}+sin{\varphi }_i\widehat{y}$$ E1

that is:

$$E_i=E_i^0{e}^{j{k}_i{}_{\cdot }\rho }=\left[E_{v,i}^{0}v+E_{h,i}^{0}h\right]e^{j{k}_i\cdot \rho }$$ E2

being (E⁰_v,i, E⁰_h,i) the field amplitudes for the vertical (z) and the horizontal (unit vector h lying on the xy plane) components, respectively. k is the wave number vector:

$$k_i=k_0{t}_i=\omega \sqrt{{\mu }_0{\epsilon }_0}{t}_i$$ E3

and (μ₀, ε₀) represent the vacuum dielectric constant and its permeability. As the propagation takes place in the xy plane, the observed directions are denoted by a unit vector $$\widehat{t}$$ or by the corresponding azimuth angle ϕ according to (1). The particle is characterized by a complex dielectric relative constant:

$${\epsilon }_r=\frac{{\epsilon }_{}}{{\epsilon }_0}={\epsilon }_r^{\text{'}}+j{\epsilon }_r^{\text{'}\text{'}}$$ E4

and it is assumed for simplicity to be homogeneous with:

$${\epsilon }_r^{\text{'}\text{'}}=-\frac{\sigma }{\omega \text{}{\epsilon }_0}$$ E5

where σ is the conductivity [Sm^-1].

In far field, the scattered field behaves like a cylindrical wave:

$$E_s=\frac{1}{\sqrt{r}}F\left({\varphi }_i,\varphi \right)E_i^0{e}^{j\left(k\rho -\frac{\pi }{4}\right)}$$ E6

where F(ϕ_i,ϕ) is the 2×2 scattering matrix accounting for the amplitude, the phase and the polarization of the scattered wave in direction ϕ when illuminated by a plane wave propagating in direction ϕ_ι. For a 2D problem the polarizations are independent and the scattering matrix is diagonal:

$$F({\varphi }_i,\varphi )=\left[\begin{array}{cc}{F}_{vv}({\varphi }_i,\varphi )& 0\\ 0& F_{hh}({\varphi }_i,\varphi )\end{array}\right]$$ E7

Therefore, for each polarization p (p = v, h) we can write:

$$E_{p,s}=\frac{1}{\sqrt{r}}{F}_{pp}\left({\varphi }_i,\varphi \right)E_{p,i}^0{e}^{j\left(k\rho -\frac{\pi }{4}\right)}$$ E8

In the following, where not stated otherwise the subscript p is dropped for simplicity of notation.

2.1. Scattering and absorption cross widths

The power dP_s scattered along direction ϕ within the differential width dl subtended by the differential angle dϕ,dl = r dϕ is:

$$d{P}_s=S_{s}dl$$

where S_s [Wm^-1] is the amplitude of the Poynting vector of the scattered wave, $$S_s={\left|E_S\right|}^2/2\eta$$ and $$\eta =\sqrt{{\mu }_0/{\epsilon }_0}$$ is the vacuum electromagnetic impedance. It follows from (8) that:

$$d{P}_s={\left|F({\varphi }_i,\varphi )\right|}^2{S}_{i}d\varphi$$ E9

and integrating over all angles:

$$P_s={\displaystyle \underset{0}{\overset{2\pi }{\int }}d}{P}_s={\displaystyle \underset{0}{\overset{2\pi }{\int }}{\left|F({\varphi }_i,\varphi )\right|}^2{S}_{i}d\varphi =S_i}{\sigma }_s$$

where the scattering cross width σ_s is defined as:

$${\sigma }_s=\frac{P_s}{S_i}={\displaystyle \underset{0}{\overset{2\pi }{\int }}{\left|F({\varphi }_i,\varphi )\right|}^{2}d\varphi }$$ E10

This quantity represents the equivalent width that would produce the amount P_s of scattered power if illuminated by a wave with power density S_i.

The geometric cross width σ_g [m] of a particle is its geometric width projected onto a plane that is perpendicular to the direction of the incident wave t_i. The relationship between the geometric and scattering widths can be investigated in two regimes. If the size of the object D (maximum distance between two points inside the object) is much smaller than the wavelength λ, it follows that:

$$\frac{{\sigma }_s}{{\sigma }_g}<<1(D<<\lambda )$$

According to Rayleigh scattering theory [19], it means meaning that in this regime the power scattered by the particle is much smaller than the product of geometric cross width and the amplitude of the Poynting vector. Besides, in the high frequency regime D>>λ:

$$\frac{{\sigma }_s}{{\sigma }_g}\to 1$$

which is known as geometric optics limit.

Similarly to (10), the absorption cross width can be defined as the ratio between the absorbed power P_a and the incident Poynting vector S_i. From Ohm’s law [20]:

$$P_a=\frac{1}{2}\omega {\displaystyle \underset{A}{\int }d\rho }{\epsilon }_r^{\text{'}\text{'}}{\left|E_{int}(\rho )\right|}^2$$ E11

where E_int(ρ) is the internal field within the particle and A is the particle’s area.

2.2. Extinction cross width and albedo

The extinction cross-section σ_ext of a particle is defined as:

$${\sigma }_{ext}={\sigma }_s+{\sigma }_a$$ E12

and it represents the total power loss from the incident wave due to scattering and absorption. The fraction of scattering over the extinction cross width is defined as albedo α:

$$\alpha =\frac{{\sigma }_s}{{\sigma }_{ext}}\le 1$$ E13

The computations of the extinction cross width (12) can be carried out from the knowledge of the scattering matrix, as stated by the forward scattering theorem (also known as optical theorem):

$${\sigma }_{ext}=-4\sqrt{\frac{\pi }{2k_0}}Re\left(F({\varphi }_i,\varphi )\right)$$

This result can be proved either undertaking the explicit computation of (11) [15] or computing the received power over a given width and relating this quantity to the geometric dimensions [19]. Besides, it can be shown that for the high frequency regime [9]:

$$\frac{{\sigma }_a}{{\sigma }_g}\to 1$$ E14

$$\frac{{\sigma }_{ext}}{{\sigma }_g}\to 2$$ E15

Equation (15) is also known as extinction paradox.

2.3. Example: scattering from a circular cylinder

Figure 2 shows the scattering, the absorption and the extinction cross widths are computed for a circular cylinder of radius a [21]. In figure 2 these cross width values are normalized over the geometric cross width σ_g=2a and they are plotted versus the radius a normalized over the wavelength λ. The particle is characterized by $${\epsilon }_r\text{'}=4$$ and σ=10⁻³. One could notice that as the size of the particle increases the absorption and scattering cross widths tend to σ_g width; whereas the extinction cross width tends to 2σ_g (extinction paradox, (15)). Also, Figure 2 distinguishes three different regions: Rayleigh scattering (a<< λ), Mie scattering (a ≈λ) and optical region (a>> λ).

Figure 2.

Scattering, absorption and extinction cross widths normalized over the geometric cross width σ_g =2a versus a/λ ( $${\epsilon }_r\text{'}=4$$ and σ= 10⁻³).

2.4. Specific intensity

Transport theory deals with the propagation of energy in a medium containing randomly placed particles. For a given point ρ and a given direction specified by vector $$\widehat{t}=cos\varphi \widehat{x}+sin\varphi \widehat{y}$$ (or equivalently by the azimuth angle ϕ), the power flux density within a unit frequency band centered at frequency f within a unit angle is defined as specific intensity and denoted by I(ρ,ϕ) [W m⁻¹rad⁻¹Hz^-1]. Hence, the amount of power dP flowing along direction ϕ within an angle dϕ through an elementary width dl with normal that forms an angle Δϕ with ϕ (see Figure 3) in a frequency interval (f, f+df) is:

$$dP=I(\rho ,\varphi )cos\Delta \varphi dld\varphi df$$ E16

The specific intensity I(ρ,ϕ) as it appears in (16) could be related either to the power emitted from a surface or to the power received by the unit width. As far as the single particle of § 2.3 is concerned, the specific intensity carried by the incident plane wave (2) is:

$$I_i(\rho ,\varphi )=S_i\delta (\varphi -\varphi )$$ E17

in which δ(.) is the Dirac Delta function, whereas, recalling (9), the specific intensity for the scattered wave can be written as:

$$I_s(\rho ,\varphi )={\left|F({\varphi }_i,\varphi )\right|}^2{S}_i={\left|F({\varphi }_i,\varphi )\right|}^2{I}_i(\rho ,{\varphi }_i)$$ E18

Hence, the square modulus of the scattering function F(ϕ_i,ϕ) relates the incident specific intensity with the scattered specific intensity. In a random medium, the specific intensity is computed as the ensemble average of the power per unit angle, frequency and area over the distribution of the random scatterers.

Figure 3.

Geometry used for the definition of specific intensity: the amount of power dP flows along direction $$\widehat{\varphi }$$ within an angle dϕ through an elementary width dl whose normal $$\widehat{n}$$ forms an angle Δϕ with the direction of propagation $$\widehat{\varphi }$$ .

3. Radiative transfer equation

The radiative transfer equation is an integro-differential equation that governs the propagation of specific intensity within a random medium. Let us assume that the random medium is a made of uniform slabs in the x-direction (see Figure 4). Therefore, specific intensity is a function of position in space only through y. In the following the derivation of the radiative transfer equation is briefly reminded. The reader is referred to [15, 17] for further details.

Figure 4.

Geometry of the propagation environment.

3.1. Reduced intensity

Let us consider an area with side 1 in the x-direction and dy in the y-direction The area dy contains Ndy particles, where N is the number of particle per unit area. The particles are uniformly distributed in space. Each particle absorbs the power P_a = σ_aI and scatters the power P_s = σ_sI, so that the decrease of specific intensity due to area dy is (recall (12)):

$$sin\varphi \cdot dI(y,\varphi )=-dyN{\sigma }_{t}I(y,\varphi )=m-dy\cdot k_{ext}I(y,\varphi )$$ E19

where k_est = Nσ_ext [m^-1] is the extinction coefficient. Notice that if there are different kinds of particles (say m), each with possibly different orientation, density N_j and extinction cross width σ_ext,j [9]:

$$k_{ext}={\displaystyle \sum _{j=1}^m{N}_j<{\sigma }_{ext,j}}>$$ E20

where < > represents the ensemble average over the distribution of particles orientations. Equation (19) defines the so-called reduced intensity since it only takes into account the extinction of the incident wave [9].

3.2. Independent scattering and limits of the transfer theory

Scattering of waves impinging on the area dy from all directions ϕ_i increases the intensity along direction ϕ according to (18). As the scatterers are assumed to be independent, the specific intensities due to different particles can be added. This is strictly true under some conditions that are discussed in the following as this point entails the limits of applicability of the radiative transfer theory [15]. As the Maxwell equations are linear, the total scattered field E can be written as the sum of the E_j fields scattered from each particle:

$$E={\displaystyle \sum _j{E}_j}$$ E21

The specific intensity is proportional to the square modulus of the electric field averaged over the distribution of scatterers:

$$<{\left|E\right|}^2>={\displaystyle \sum _j<{\left|E_j\right|}^2>}+{\displaystyle \sum _j{\displaystyle \sum _{l\ne j}<E_j}}{E}_l^{*}>$$

Now, let $$E_j=\left|E_j\right|e^{j{\alpha }_j}$$ , we get that:

$$<E_j{E}_l^{*}>=<\left|E_j\right|\left|E_l\right|e^{j({\alpha }_j-{\alpha }_l)}$$

The phase difference (α_j− α_l) depends on the distance between the particles d_jl through the product k_0djl. If the distribution of particle separation is not much smaller than the wavelength, i.e. the standard deviation (S.D.) of d_jl satisfies:

$$S.D.(d_{jl})\ge \frac{\lambda }{4}$$ E22

then (α_j− α_l) is approximately uniformly distributed within [0, 2π] so that:

$$<E_j{E}_l^{*}>=0$$

From the discussion above, the assumption underlying the radiative transfer equation of independent scattering limits the applicability of the transfer theory to cases where the distance between particles is large enough (see (22)) so as to make negligible the near far interactions between particles. Experimental studies confirm this conclusion: for the radiative transfer to be applicable the spacing between scatterers must be larger than λ/3 and 0.4D where λ is the wavelength of the host medium and D the diameter of the scatterers [10]. Theoretical studies on the relationship between the radiative transfer approach and the wave approach using Maxwell’s equations can be found in [22-23]. Moreover, the limits of the transfer theory are investigated for a two-dimensional problem similar to the one considered here in [16] through comparison with the wave approach.

Assuming independent scattering, the increase on the specific intensity along direction ϕ due to scattering within the area dy is:

$$dy{\displaystyle \underset{0}{\overset{2\pi }{\int }}N{\left|F({\varphi }_i,\varphi )\right|}^{2}I(\rho },{\varphi }_i)d{\varphi }_i=dy{\displaystyle \underset{0}{\overset{2\pi }{\int }}p({\varphi }_i\varphi )I(\rho },{\varphi }_i)d{\varphi }_i$$ E23

where $$p({\varphi }_i\varphi )=NF({\varphi }_i,\varphi )$$ is the phase function. Notice that if the random medium contains particles of different kinds, the overall phase function is defined as [9]:

$$p({\varphi }_i,\varphi )={\displaystyle \sum _{j=1}^m{N}_j<}{p}_j({\varphi }_i,\varphi )>$$ E24

3.3. Progressive and regressive intensity

The radiative transfer equation is obtained by combining (19) and (23):

$$sin\varphi \frac{\partial I(y,\varphi )}{\partial y}=-k_{ext}I(y,\varphi )+{\displaystyle \underset{0}{\overset{2\pi }{\int }}p(y,{\varphi }_i,\varphi )}I(y,{\varphi }_i)d{\varphi }_i$$ E25

In this formulation, the extinction coefficient and the phase function are considered function of the position y. A first step toward the solution of scalar radiative transfer problem is converting equation (25) into two coupled integro-differential equations by introducing the progressive intensity I⁺, that corresponds to propagating directions 0 < ϕ < π, and the regressive intensity I^-, that accounts for the propagating directions π < ϕ < 2π. Τhe two specific intensities are defined as:

$$\begin{array}{l}{I}^{+}(y,\varphi )=I(y,\varphi )\\ I^{-}(y,\varphi )=I(y,\varphi +\pi )\end{array}$$

where the azimuth direction ϕ ranges within 0 < ϕ< π. The scalar radiative transfer equation (25) can be equivalently stated as:

$$\begin{array}{c}sin\varphi \frac{\partial I^{+}(y,\varphi )}{\partial y}=-k_{ext}{I}^{+}(y,\varphi )+{\displaystyle \underset{0}{\overset{\pi }{\int }}p(y,{\varphi }_i,\varphi )}{I}^{+}(y,{\varphi }_i)d{\varphi }_i+\\ +{\displaystyle \underset{0}{\overset{\pi }{\int }}p(y,{\varphi }_i+\pi ,\varphi )}{I}^{-}(y,{\varphi }_i)d\varphi \\ -sin\varphi \frac{\partial I^{-}(y,\varphi )}{\partial y}=-k_{ext}{I}^{-}(y,\varphi )+{\displaystyle \underset{0}{\overset{\pi }{\int }}p(y,{\varphi }_i,\varphi +\pi )}{I}^{+}(y,{\varphi }_i)d{\varphi }_i+\\ +{\displaystyle \underset{0}{\overset{\pi }{\int }}p(y,{\varphi }_i+\pi ,\varphi +\pi )}{I}^{-}(y,{\varphi }_i)d\varphi \end{array}$$ E26

This equivalent formulation of (25) makes it easier to set the boundary conditions as explained in §3.3.1 and §3.3.2. For a uniform distribution of the particles over the random medium and circular cylindrical particles, the phase matrix becomes a function ϕ-ϕ_i only and the extinction matrix becomes independent on ϕ [17]. Therefore, in this case the scalar radiative transfer equations (26a-26b)can be written as:

$$\begin{array}{l}sin\varphi \frac{\partial I^{+}(y,\varphi )}{\partial y}=-k_{ext}{I}^{+}(y,\varphi )+{\displaystyle \underset{0}{\overset{\pi }{\int }}{p}^{+}(y,\varphi -{\varphi }_i)}{I}^{+}(y,{\varphi }_i)d{\varphi }_i+\\ +{\displaystyle \underset{0}{\overset{\pi }{\int }}{p}^{-}(y,\varphi -{\varphi }_i)}{I}^{-}(y,{\varphi }_i)d\varphi \end{array}$$ E27

$$\begin{array}{l}-sin\varphi \frac{\partial I^{-}(y,\varphi )}{\partial y}=-k_{ext}{I}^{-}(y,\varphi )+{\displaystyle \underset{0}{\overset{\pi }{\int }}{p}^{-}(y,\varphi -{\varphi }_i)}{I}^{+}(y,{\varphi }_i)d{\varphi }_i+\\ +{\displaystyle \underset{0}{\overset{\pi }{\int }}{p}^{+}(y,\varphi -{\varphi }_i)}{I}^{-}(y,{\varphi }_i)d\varphi \end{array}$$ E28

where $$p^{+}(y,\varphi -{\varphi }_i)=p(y,{\varphi }_i,\varphi )$$ and $$p^{-}(y,\varphi -{\varphi }_i)=p(y,{\varphi }_i,\varphi +\pi )=p(y,\varphi +{\pi }_i,\varphi )$$ , based on the 2π periodicity of both p⁺ and p^- [17].

3.3.1. Boundary conditions on the specific intensity

The radiative transfer equation has to be solved by imposing appropriate boundary conditions. Here, the boundary conditions that the specific intensity must satisfy on a plane boundary between two media with indices of refraction n₁ and n₂ are considered. The expression for the transmission T and reflection coefficient R for the two polarizations are known in literature [20]. The reflected specific intensity I_r is related to the incident specific intensity I_i as:

$$I_r={\left|R\right|}^2{I}_i$$

As far as the transmitted specific intensity is concerned, we can write the conservation of power on a segment dl belonging to the boundary (recall (16)):

$$I_{i}cos{\varphi }_{1}d{\varphi }_{1}dl=I_{r}cos{\varphi }_{1}d{\varphi }_{1}dl+I_{t}cos{\varphi }_{2}d{\varphi }_{2}dl$$ E29

As by the Snell’s law $$n_{1}sin{\varphi }_1=n_{2}dsin{\varphi }_2$$ , it ends:

$$\frac{I_t}{I_i}=\frac{n_2}{n_1}\left(1-{\left|R\right|}^2\right)$$

4. Solution of the radiative transfer equation through numerical quadrature

The radiative transfer equation is an integro-differential equation whose solution in analytical form is very difficult, if not impossible. However, efficient numerical solutions can be devised. A comprehensive treatment of the main techniques can be found in [15, 24].

The case under study concerns a random medium where relevant scattering occurs. An approximate solution can be obtained by computing the integrals in (26) by numerical quadrature as firstly proposed in [15]. The continuum of propagation directions ϕ is discretized into a set of n directions ϕ_iI = (i-1) Δϕ, where Δϕ = π/n and i = 1,..., n, and the corresponding n×1 vectors gathering the progressive and regressive intensities:

$$I^{+}(r)=\left[\begin{array}{c}{I}^{+}(h,{\varphi }_1)\\ I^{+}(h,{\varphi }_2)\\ ⋮\\ I^{+}(h,{\varphi }_n)\end{array}\right]$$ E30

$$I^{-}(r)=\left[\begin{array}{c}{I}^{-}(h,{\varphi }_1)\\ I^{-}(h,{\varphi }_2)\\ ⋮\\ I^{-}(h,{\varphi }_n)\end{array}\right]$$ E31

The two integro-differential equations (27) can be approximate as follows:

$$\frac{d{I}^{+}(y)}{dy}=-K(y)I^{+}(y)+P^{+}(y)I^{+}(y)+P^{-}(y)I^{-}(y)$$ E32

$$\frac{d{I}^{-}(y)}{dy}=-K(y)I^{-}(y)-P^{-}(y)I^{+}(y)-P^{+}(y)I^{-}(y)$$ E33

where then ×n matrices P⁺ and P^- are defined as:

$${\left[P^{\pm }(y)\right]}_{i,j}=\frac{1}{sin{\varphi }_i}{p}^{\pm }(y,{\varphi }_i-{\varphi }_j)\Delta \varphi$$

so that:

$$\frac{1}{sin{\varphi }_i}{\displaystyle \underset{0}{\overset{\pi }{\int }}{p}^{\pm }(y,\varphi -\varphi \text{'})I^{\pm }(y,\varphi \text{'})d{\varphi }^{\text{'}}\approx {\displaystyle \sum _{j=1}^n{\left[P^{\pm }(y)\right]}_{i,j}}{I}^{\pm }(y,{\varphi }_j)}$$ E34

The n×n K(y) matrix is defined as:

$$\begin{array}{c}{\left[K(y)\right]}_{i,j}=\frac{1}{sin{\varphi }_i}{k}_{ext}(y)i=j\\ =0i\ne j\end{array}$$ E35

Now, defining the 2n×1 vector $$I(y)=\left[I^{+}{(y)}^T{I}^{-}{(y)}^T\right]$$ , one obtains the system of first order linear equation:

$$\frac{dI(y)}{dy}=G(y)I(y)$$ E36

where

$$G(y)=\left[\begin{array}{cc}-K(y)+P^{+}(y)& P^{-}(y)\\ -P^{-}(y)& K(y)-P^{+}(y)\end{array}\right]$$

4.1. Discrete ordinate Eigen analysis

If the space contains slabs of homogeneous random (and non-random) media, the linear differential system (36) within each slab (say the l ^th) can be written as

$$\frac{dI(y)}{dy}=G_{l}I(y)$$ E37

being matrix G _l a constant. This follows from the fact that, the extinction coefficient and the phase function within each slab are constant. In this case, the system (37) can be solved by the discrete ordinate analysis method as follows. Let $$I(y)=\beta exp(-\lambda g)$$ be a tentative solution where β is a 2n×1, one get by substitution in (36) - dropping l for simplicity:

$$(\lambda U+G)\beta =0$$

As a result, (β,λ) represent pair of eigenvectors/eigenvalues of the matrix –G. Notice that given the symmetry relations:

$$\begin{array}{l}{\left[G_l\right]}_{(n+i))n+i)}=-{\left[G_l\right]}_{i,i}\\ {\left[G_l\right]}_{(n+i))n+j)}=-{\left[G_l\right]}_{iji}i,j=\mathrm{1,}\cdots ,n\end{array}$$

by ordering the eigenvalues in increasing order ( $${\lambda }_i\ge {\lambda }_{i-1}$$ ) it follows:

$${\lambda }_{(n+i)}=-{\lambda }_{i}i=\mathrm{1,}\cdots ,n$$

and the corresponding eigenvectors satisfy the condition:

$${\left[{\beta }_{(n+i)}\right]}_{(n+)}={\left[{\beta }_i\right]}_{j}i=\mathrm{1,}\cdots ,n$$

The solution of (36) can be written as the linear combination

$$I(y)={\displaystyle \sum _{i=1}^{2n}{\beta }_i{e}^{(-{\lambda }_{i}y)c_i}}=BD(y)c$$ E38

where $$B=\left[{\beta }_1\cdots {\beta }_{2n}\right]$$ , $$c={\left[c_1\cdots c_{2n}\right]}^T$$ , and D(y) is a diagonal matrix with elements $$\left[e^{(-{\lambda }_{1}y)}\cdots e^{(-{\lambda }_{2n}y)}\right]$$ .

The vector of unknown constant c can be computed from the knowledge of I(y) for one value of y, say $$\overline{y}$$ , as:

$$c={\left(BD(\overline{y})\right)}^{-1}I(\overline{y})$$ E39

4.1.1. Setting the boundary conditions for a single slab

Given geometry depicted in figure 4, with only one slab of random medium (ranging within 0 ≤ y ≤d), the semi-infinite medium along the positive direction of y is a dielectric one with given dielectric constant. The boundary conditions read:

$$I^{+}(0)={\overline{I}}_0^{+}$$ E40

$$I^{-}(d)=A_R{\overline{I}}^{+}(d)$$ E41

where $${\overline{I}}_0^{+}$$ is the incident specific intensity and A_R is the n×n reflection matrix relative to the interface at y = d between the random medium and the last slab. This can be easily computed from the results in §2.4. The value of the unknown constant vector (39) for the specific intensity within the random medium is obtained through the computation of the specific intensity $$I(\overline{y})$$ . To do this, as first step one need to calculate from (38) and (39):

$$I(d)=\underset{T}{\underbrace{BD(d)B{}^{-1}}}I(\overline{y})I(0)$$ E42

Then, using (40) and (41), and partitioning the matrix T as:

$$T=\left[\begin{array}{cc}{T}_{11}& T_{12}\\ T_{21}& T_{22}\end{array}\right]$$ E43

in which T_ij is n×n, one gets:

$$I{(0)}^{-}={\left(A_R{T}_{12}-T_{22}\right)}^{-1}\left(T_{21}-A_R{T}_{11}\right){\overline{I}}_0^{+}$$

The substitution of this finding in (39) is used for the computation of the vector c.

4.1.2. Setting the boundary conditions for a multi-layer medium

In the multi-layered geometry of figure 4, each slab has depth dl and the specific intensity is (38). The unknown is the constant vector c _l and the subscript l runs over the different slabs, l = 1,..., L. Notice that slab can contain either a random medium characterized by its matrix G_l or vacuum (matrix G_l=0 for vacuum). The semi-infinite medium along the positive direction of y is a general dielectric with given dielectric constant. The boundary conditions (40) and (41) still hold with $$d={\displaystyle \sum _{l=1}^L{d}_l}$$ .

The specific intensity I_l,2 - at the interface between the l and the (l + 1)^th slab - is related to the specific intensity I_l,1 at the interface between the (l – 1) and l ^th slab as:

$$I_{l\mathrm{,2}}=\underset{T_l}{\underbrace{B_{l}D(d_l)B{}_l{}^{-1}}}I(\overline{y})I_{l\mathrm{,1}}$$ E44

Therefore, similarly to (42) we can write:

Proceeding as in the previous Section, the specific intensity I(0) is now calculated and used in (39) to obtain the constant vector c₁ for the first slab. Then, equation (44) can be iteratively applied so as to compute the specific intensity at the beginning of each slab. The latter allows through (39) the computation of c_l for each slab.

5. Numerical results

In this section two examples illustrate how radiative transfer theory could be employed to study the beam broadening and the corresponding spatial correlation for an indoor environment within the context of a communication system.

The main assumption is that the propagation environment can be modeled as a layered parallel plane medium. This situation is a first approximation of an open space office made of a succession of tables where each table is modeled as random medium containing randomly placed scatterers. These are modeled as circular cylinder so as to simplify the solution of the radiative transfer equation.

5.1. Case study I: one table

In this example, the geometry under study is depicted in Figure 5. The whole xy plane is characterized by the vacuum dielectric constant ε₀. One slab ranging within 1 ≤ y ≤ 1.8 contains N [m⁻²] uniform randomly distributed circular cylinders with radius a = 6 cm. Where not stated otherwise, frequency of operation is f = 5.2 GHz, according to standard wireless local area networks (WLAN) such as IEEE802.11x and Hyperlan/2.

Figure 5.

Geometry of the propagation environment (case study I).

5.1.1. Beam broadening

Figure 6 shows the specific intensity I⁺(y,ϕ) for N=10, ε^’_r =4, σ=10⁻³ Sm^-1 and vertical polarization. Specifically, the image represents the beam broadening of a plane wave propagating in the direction ϕ_i = 90° with specific intensity $$I_0^{+}$$ = 1 [W m⁻¹rad⁻¹Hz].

The specific intensities are assumed to be normalized with respect to $$I_0^{+}$$ and thus are shown in dB. Right after the entrance into the table region the progressive beam I⁺(y, ϕ) broadens since the energy is scattered in all directions by the cylinders. As expected the regressive beam I⁻(y, ϕ) is zero on the right side of the table and almost uniform in the yϕ plane (not shown).

Figure 6.

Specific intensity I⁺(y,ϕ) versus y (vertical polarization, ε^’_r = 4, σ = 10⁻³ Sm^-1, N = 10).

In order to get a quantitative insight into the beam broadening discussed above, Figure 7 shows the specific intensity I⁺(y, ϕ) for y = 1.9 m (on the right of the table) for different values of the density N=2, 10 (vertical polarization). The reduced or line of sight (LOS) contribution as shown in the box decreases for increasing object densities. However, a larger density of scatterers entails a more relevant contribution of the diffuse energy (i.e., I⁺(y, ϕ) for ϕ ≠ 90°).

Figure 7.

Specific intensity I⁺(y, ϕ) for y = 1.9 m (on the right of the table) for different values of the density N = 2, 10 (vertical polarization, ε^’_r = 4, σ=10⁻³ Sm^-1).

Figure 8.

Specific intensity I⁺(y, ϕ) for y = 1.9 m (on the right of the table) for different values of the frequency f = 5.2, 52 GHz (vertical polarization, N = 6, ε^’_r = 4, σ = 10⁻³ Sm^-1).

The effect of an increase in the carrier frequency, envisioned for next generation wireless LAN, is shown in Figure 8. A larger carrier frequency (f = 52 GHz) yields a less consistent contribution of the LOS direction (see box) and a greater number of interference fringes. Employing the horizontal polarization instead of the vertical one yields qualitatively similar results as shown in Figure 9 for N=6. Therefore, in the following only the vertical polarization is considered.

Figure 9.

Specific intensity I⁺(y, ϕ) for y = 1.9 m (on the right of the table) for horizontal and vertical polarization frequency (f = 5.2 GHz, N = 6, ε^’_r = 4, σ = 10⁻³ Sm^-1).

5.1.2. Spatial correlation

In a communication link, it is of great interest to assess the degree of correlation between the signals received by different antennas as a function of their inter-spacing Δ [m]. In fact, wireless links can capitalize on the uncorrelation between the received signals to increase the degree of diversity of the system that in turn rules the asymptotic performance of the link in terms of probability of error [25].

The correlation r_y(Δ) of the signals received by antennas separated by Δ can be expressed in terms of the probability density function p_y(ϕ) of the direction of arrivals of the waves impinging on the receivers. This function can be obtained by interpreting (after appropriate scaling) the power received over a certain direction ϕ as a measure of the probability that a signal is received through direction ϕ:

$$p_y(\varphi )=\frac{I^{+}(y,\varphi )}{{\displaystyle \int I^{+}(y,\varphi )d\varphi }}$$ E45

Then, recalling that the ratio between signals received at two points separated by Δ in the direction x is $$exp(-j2\pi /\lambda cos\varphi \Delta )$$ :

$$r_y(\Delta )={\displaystyle \underset{0}{\overset{\pi }{\int }}{p}_y(\varphi )e^{\left(-\frac{j2\pi }{\lambda }cos\varphi \Delta \right)}}d\varphi$$ E46

The correlation (46) is evaluated for the example at hand under the same assumptions as in Figure 7.

Figure 10.

Correlation r_y(Δ) versus Δ/ λ for y = 1.9 m (on the right of the table) for different values of the density N = 2, 10 (vertical polarization, ε^’_r = 4, σ = 10⁻³ Sm^-1).

The results are shown in figure 10. The correlation decreases with increasing object density. Thus, from the perspective of communication system performance, increasing the object density is beneficial in terms of degree of diversity. For instance, for Δ = λ the correlation decreases from around 0.9 to around 0.7.

5.2. Case study II: two tables with interface

In this example the environment has a more complicated geometry, where two tables are followed by a semi-infinite dielectric slab (dielectric constant ε_r,I) as shown in figure 11. The simulation parameters follow the setting described for Case study I, in particular the radius of the scatterers is a=6 cm and the carrier frequency is f=5.2 GHz. The dielectric constant of the slab is chosen as ε_r,I =3.

Figure 11.

Geometry of the propagation environment (case study II).

5.2.1. Beam broadening

Figure 12 shows the specific intensity I⁺(y, ϕ) for N = 6, ε^’_r = 4, σ = 10⁻³ Sm^-1 and ϕ_i=90° and vertical polarization. Again, the specific intensity is normalized with respect to $$I_0^{+}$$ and thus are shown in dB. As discussed for case study I, right after the entrance into the table regions the progressive beam I⁺(y, ϕ) broadens since the energy is scattered in all directions by the cylinders. Moreover, the refraction over the interface at y=4.6 m focuses the beam within the semi-infinite dielectric, hence reducing its angular spread. The refraction and the reflection of the reduced (or LOS) component are apparent in Figures 12 and 13, respectively. The incidence angle is chosen equal to 45°. Refracted and reflected angles satisfy the Snell’s law.

Figure 12.

Specific intensity I⁺(y, ϕ) versus y (vertical polarization, ε^’_r = 4, σ = 10⁻³ Sm^-1, N = 6, ϕ_i=90°).

Figure 13.

Specific intensity I⁺(y, ϕ) versus y (vertical polarization, ε^’_r = 4, σ = 10⁻³ Sm^-1, N = 6, ϕ_i=45°).

Figure 14.

Specific intensity I^-(y, ϕ) versus y (vertical polarization, ε^’_r = 4, σ = 10⁻³ Sm^-1, N = 6, ϕ_i=45°).

5.2.2. Spatial correlation

The spatial correlation r_y(Δ) is shown for the same parameters of the previous example in Figure 15. It is computed for y = 1.9 m (just on the right of the first table) and y = 3.7 m (just on the right of the second table). As expected, the increased scattering contribution due to the random scatterers of the second table decreases the spatial correlation. In particular, for inter-antenna spacing around Δ = 5λ the correlation is decreased from 0.75 to around 0.61.

Figure 15.

Correlation r_y(Δ) versus Δ/ λ for y = 1.9 m (on the right of the first table) and y = 3.7 m (on the right of the second table) for N = 6, vertical polarization, ε^’_r = 4, σ=10⁻³ Sm^-1, ϕ_i=90°.

5.3. Limits and validity of radiative transfer predictions

The validation of numerical results presented in §5.1 and §5.2 relies on a discussion, as there is no empirical evidence (measurements) to compare with. The two chosen cases are not realistic but they represent a possible benchmark to be employed against the numerical results that issue from analytical methods. Among these, the T-matrix approach [26, 27] seems to be the best suited to compute scattering from a random distribution of cylinders and to compare the radiated fields on a realization-by-realization basis with ray tracing. At this stage of development, the only comparison that the author is able to provide to assess the validity of a radiative transfer approach is against a ray tracing technique based upon the beam tracing method [26]. The numerical code used for the simulations was developed at the Politecnico di Milano [29] and it has already been used for different purposes related to the indoor propagation [30].

To make the beam tracing procedure suitable for the study of electromagnetic propagation, it is necessary to include the reflection coefficient associated with the interactions of the path with the environment. In particular, the signal received from each path has to be scaled by the product of the reflection coefficients corresponding to the bounces each path goes through when propagating from the transmitter to the receiver. To compare the outcome of the ray tracing simulation with the results of radiative transfer theory some approximations were made. Case I was reproduced in such a way to replicate the same geometrical and radio electrical conditions. The reader is referred to [31] for further details about the numerical computation and the comparison results, while here the author recalls the general approach and the basic information. The transmitted plane wave was approximated by building a linear antenna array of N_T elements with inter-element spacing Δ_T lying parallel to the y axis. To fulfill the requirement of far field regime, the transmitting array was placed at a range distance $$R>{\left(N{\Delta }_T\right)}^2/\lambda$$ [9].

The circular scatterers were approximated by polygons of N_P sides. The quantities of interest were averaged over N_I realizations of the random medium. The signal is received at the desired points by a linear antenna array of N_R elements with inter-element spacing Δ_R lying parallel to the y axis. Figure 16 illustrates the geometry together with the transmitting and the receiving antenna arrays.

Figure 16.

Geometry of the experiment studied by means of ray tracing.

The specific numerical values for the ray-tracing algorithm were selected after careful empirical investigation to yield negligible approximation errors. The transmitting array has N_T= 100 elements and inter-element distance Δ_T=λ/2, while the receiving array has N_R= 80 elements with inter-element spacing Δ_R = λ/4, thus resulting in an angular resolution of approximately Δϕ=6°. The key parameter that discriminates the reliability of the radiative transfer approach was identified as the fraction of area that is effectively occupied by the scatterers. Let $$A_{ext}=\pi {({\sigma }_{ext}/2)}^2$$ be the effective area occupied by each particle. The fraction of effective area that is occupied by the scatterers is $$\eta =N{A}_{ext}$$ . The numerical investigation proved that for scatterers that yield a relatively low value of η the radiative transfer provides a solution that closely matches the second order statistics (i.e. power fluctuations) as given by the ray tracing results.

Figure 17 shows the progressive specific intensity I ⁺ (y,ϕ) computed according to radiative theory compared with $$\widehat{I}(y,\varphi )$$ from the ray tracing procedure. To compare the two approaches, the effect of limited angular resolution is accounted for by processing the outcome of the radiative transfer I⁺(y,ϕ) and the number of averaging iterations for ray tracing is N_I = 100 [31]. By increasing the density N, and consequently the fraction η of effective area occupied by the scatterers, the difference between the prediction of radiative transfer and those from ray tracing increases.

Figure 18 shows the comparison between the spatial correlation r _y (Δ) computed by the radiative transfer theory and the same quantity $${\widehat{r}}_y(\Delta )$$ from ray tracing versus Δ/λ for y = 1.9 m (on the right of the table) and two values of the density N=6,10.

Figure 17.

Specific intensity I⁺(y,ϕ) computed according to the radiative theory with and without limited resolution and $$\widehat{I}(y,\varphi )$$ from the ray tracing procedure for ε^’_r = 4, σ = 10⁻³ Sm^-1, y = 1.9 m (on the right of the random medium), N = 10 m⁻² and vertical polarization.

Figure 18.

Spatial correlation r_y(Δ) computed by the radiative transfer theory and $${\widehat{r}}_y(\Delta )$$ from ray tracing versus Δ/λ for y = 1.9 (on the right of the random medium) and different values of the density N = 6, 10 m⁻² (vertical polarization and ε^’ _r = 4, σ = 10⁻³ Sm^-1).

6. Conclusions

In this work, the use of radiative transfer theory to study the propagation in an indoor environment was reported. A particular focus was made on two specific aspects related to the channel performances, such as the beam broadening and the spatial correlation, as they are parameters of interest when deploying a network. The pretty good agreement, although under well-controlled conditions, encourages proceeding along two distinct directions. From one side, an improvement in the modeling and in the sketching of actual environments is mandatory if one wants to use the radiative transfer technique as a possible approach to analyze indoor channel performances in practical scenarios. From the other side, validation should be performed either against measurements or - at least - against appropriate analytical solution without limitations or approximations in the frequency domain (or wavelength scale) of interest.

Acknowledgments

The author would like to gratefully acknowledge the work of Osvaldo Simeone (now Assistant Professor at the New Jersey Institute of Technology), who contributed with ideas, numerical developments and simulations to this research during his doctoral studies in Milano at the Politecnico, Faculty of Telecommunication Engineering.