We define radio maps as path integrals, similar to how measurements are defined in computer graphics [23, 24, Section 3.7]. Consider a point source $𝐬$. The measurement surface, denoted by $ℳ$, is divided into $N_M$ cells $M_i$, for $i=1,\mathrm{\dots },N_M$, each of area $\left|M_i\right|$. Ignoring diffraction for now, let $𝒜$ denote the union of all scatterer surfaces in the scene. We define

as the set of all sequences of $\mathrm{\ell }$ vertices representing interactions with scatterers surfaces—that is, all paths of length $\mathrm{\ell }$ that are not diffracted,

A path from $𝐬$ to a point $𝐭\in ℳ$ can be written as $(𝐬,p,𝐭)$, though for brevity we refer simply to the path $p$. Note that ${𝒜}^0$ corresponds to the LoS path from $𝐬$ to $𝐭$.

Focusing on paths of a fixed length $\mathrm{\ell }$, we define the radio map value for the $i$-th cell as follows:

where $d𝐭$ is an infinitesimal surface element on the measurement surface $M_i$, and ${\mu }^{(\mathrm{\ell })}(p)$ is the area-product measure on the scattering surface, i.e.,

$𝒱\left(p\right)$ is the indicator function which equals to $1$ if the path is not occluded and $0$ otherwise, and $e(p)$ is the channel gain for the path, defined as

Here, $\lambda$ is the wavelength, ${𝐂}_{\text{T}}$ is the $2$-dimensional complex vector of the transmit antenna pattern evaluated at the departure angles $({\theta }_{\text{T}},{\phi }_{\text{T}})$ of $p$, and $𝐓$ is a $2\times 2$ complex matrix modeling the effects of scatterers and free-space path loss along $p$. Unlike in (4), the receive antenna pattern is not applied; instead, the squared norm of the electric field is used. This corresponds to assuming a receiver equipped with an isotropic antenna whose polarization is matched to the one of the incident wave. Note that since $e(p)$ is the squared norm of the electric field amplitude, the radio map represents a non-coherent sum of the path gains.

In (41), $w(p)$ restricts which paths are considered:

where ${\widehat{𝐤}}_i^{(i)}$ is the incident wave direction at the $i$-th interaction, ${\widehat{𝐤}}_s^{(i)}$ is the scattered wave direction at the $i$-th interaction, and $w^{(i)}$ depends on the types of interactions considered for the $i$-th interaction. For instance, if only specular reflections and refractions are allowed at the $i$-th interaction,

where $\widehat{𝐧}$ is the normal to the scatterer’s surface. Although $w(p)$ could in principle be absorbed into $e(p)$, as usually done in the computer graphics literature, we keep it separate for clarity.

The overall radio map value for the $i$-th cell is given by summing over all possible path depths $\mathrm{\ell }=0,1,\mathrm{\dots },\mathrm{\infty }$:

While the definition (46) is comprehensive, it is not practical for computation as-is. Indeed, directly integrating over the surfaces of the measurement cells would require connecting points on the measurement surface to the source, for example by placing one or more targets within each cell and using the path solver to compute the connecting paths. However, because accurately computing fine-grained radio maps typically demands a large number of target points, this approach is computationally prohibitive. Instead, we use Monte Carlo integration with a single SBR loop to efficiently estimate the value of each cell. However, SBR is not applicable to paths that include diffraction: as explained in Section 3.1, the probability of a ray intersecting an edge is zero. To address this, Sionna RT computes the radio map in two steps, as illustrated in Figure 26. First, the radio map is computed for all interaction types except diffraction using a single SBR loop; wedges near the source are detected and stored during this step. Second, only diffracted paths are computed. Currently, Sionna RT supports only first-order diffraction for radio map computation, i.e., paths connecting the source to the measurement surface through a single diffracting wedge. The two radio maps are then summed to obtain the final result. The following describes the computation for each step.

When diffraction is ignored, the radio map can be efficiently computed using a single SBR loop. This is achieved by reformulating the radio map definition (46) so that the integration is performed over the directions of rays originating from points that spawn ray tubes intersecting the measurement surface, rather than directly over the measurement surface itself. This is enabled by applying the change of variables $d𝐭=\frac{r^2}{\cos \theta }d\omega$, where $r$ is the length of the ray tube on the measurement surface, $\theta$ is the angle between the incident ray and the normal of the measurement surface, and $d\omega$ is the infinitesimal solid angle subtended by the ray tube (see Section 2.4). Thus, $d𝐭$ represents the intersection area of the ray tube on the measurement surface. With this change of variables, the radio map definition (46) can be rewritten as:

where $\mathrm{\Omega }$ is the set of all directions of rays originating from points that spawn ray tubes intersecting the measurement surface, typically a hemisphere or a sphere.

Expression (51) can be estimated as follows. Paths are generated by launching $N_S$ rays from the source $𝐬$ and tracing their interactions with the scene geometry, thus performing a Monte Carlo integration over the surfaces of the scatterers $𝒜$. Intuitively, when a ray intersects a scatterer surface, it samples a surface element $\mathrm{\Delta }A$ that corresponds to the ray tube intersection area, given by $\mathrm{\Delta }A=\frac{r^2}{\cos \theta }\mathrm{\Delta }\omega$, where $r$ denotes the ray tube length, $\theta$ is the angle between the incident ray and the surface normal, and $\mathrm{\Delta }\omega$ is the solid angle subtended by the ray tube. The maximum path depth is set to , and since a single path may intersect the measurement surface at multiple depths, all path depths $\mathrm{\ell }=0,1,\mathrm{\dots },L$ are considered. Interaction types (specular reflection, refraction, or diffuse reflection) are sampled according to (20).

This process leads to the following Monte Carlo estimator for the radio map (51):

where $p_n^{(\mathrm{\ell })}$ denotes the $n$-th path at depth $\mathrm{\ell }$; $𝕀\left(p_n^{(\mathrm{\ell })}\right)$ is an indicator function equal to $1$ if $p_n^{(\mathrm{\ell })}$ intersects the measurement surface at depth $\mathrm{\ell }$ and $0$ otherwise; and $\mathrm{Pr}\left({𝝌}_n^{(\mathrm{\ell })}\right)$ is the probability of sampling the specific sequence of interaction types ${𝝌}_n^{(\mathrm{\ell })}$ corresponding to $p_n^{(\mathrm{\ell })}$. ${\beta }_n^{(\mathrm{\ell })}$ is a weighting factor that incorporates the intersection area of the ray tubes on the scatterer and measurement surfaces, and thus absorbs the factor $\frac{r^2}{\cos \theta }$,

where, for the $n$-th path at depth $\mathrm{\ell }$, $n_{\text{rt}}$ is the number of ray tubes that constitute the path, $r_i$ is the length of the $i$-th ray tube, ${\theta }_i$ is the angle between the incident ray tube and the surface normal, and $\mathrm{\Delta }{\omega }_i$ is the solid angle subtended by the $i$-th ray tube. In the current propagation model of Sionna RT, only the source $𝐬$ and diffuse reflection points act as point sources that spawn new ray tubes, whereas specular reflection and refraction points merely deviate existing ray tubes rather than spawning new ones. This is illustrated in Figure 27. The source emits $N_S$ ray tubes, each subtending a solid angle of $\frac{4\pi }{N_S}$, while each diffuse reflection point spawns a single ray tube, uniformly sampled over the hemisphere, with a solid angle of $2\pi$. Thus, $n_{\text{rt}}$ equals the number of diffuse reflection points along the path plus one.

The Sionna RT radio map solver currently supports only first-order diffraction, i.e., paths connecting the source to the measurement surface via a single diffracting wedge. For clarity, we consider diffraction on a single edge $ℰ$, defined by an origin $𝐨$, an edge vector $\widehat{𝐞}$, and length $L_{ℰ}$. Sionna RT processes all such wedges found near the source (see Figure 26) in parallel.

For a single diffraction wedge, the radio map definition (46) is rewritten by integrating along the edge $ℰ$ instead of the scattering surface $𝒜$:

Here, the inner integral is over the diffraction point position $x$ along the edge $ℰ$, i.e., $𝐯=𝐨+x\widehat{𝐞}$, while the outer integral is over the measurement cell $M_i$. The function $w(p)$ ensures the scattered wave lies on the Keller cone:

where ${\beta }_0^{\prime }$ and ${\beta }_0$ are the angles between the edge and the incident and diffracted wave, respectively (see Figure 5(b)).

Direct integration over the measurement surface is inefficient, as it would require connecting surface points to the source (for example, by placing targets within each cell and running the path solver). To overcome this, consider an orthogonal basis $({\widehat{𝐭}}_0,{\widehat{𝐧}}_0,\widehat{𝐞})$, where ${\widehat{𝐧}}_0$ is the normal to a wedge face and ${\widehat{𝐭}}_0\times {\widehat{𝐧}}_0=\widehat{𝐞}$, as illustrated in Figure 28. For a given angle of incidence ${\beta }_0^{\prime }$, any point on $M_i$ that can be reached by a diffracted ray can be written as:

Here, ${𝐯}^{\prime }$ is a point on the edge $ℰ$ parametrized by $x^{\prime }$, ${\widehat{𝐤}}_s$ is the direction of the diffracted ray (with $\varphi$ the Keller cone azimuth), and $\gamma >0$ is the distance for which the ray from ${𝐯}^{\prime }$ in direction ${\widehat{𝐤}}_s$ intersects the measurement cell $M_i$ (see Figure 28). Since a ray and a plane intersect at most once, $𝐭$ is uniquely determined by $x^{\prime }$ and $\varphi$, and thus $𝐭(x^{\prime },\varphi )$ forms a reparametrization of the subset of $M_i$ that can be reached by a diffracted ray when the angle of incidence is ${\beta }_0^{\prime }$.

With this reparametrization, (54) becomes:

where $n\pi$ is the wedge exterior angle, and ${\parallel \frac{\delta 𝐭}{\delta x}\times \frac{\delta 𝐭}{\delta \varphi }\parallel }_2$ accounts for the change of variables [26].

The condition in (55) is only satisfied when $x=x^{\prime }$ for a given $\varphi$, leading to:

To estimate (58), $N_S^{\prime }$ samples are drawn on the wedge: $x$ uniformly in $[0,L_{ℰ}]$ and $\varphi$ uniformly in $[0,n\pi ]$. This yields the Monte Carlo estimator:

where $p_n$ is the $n$-th path, $𝕀\left(p_n\right)$ is 1 if $p_n$ intersects the measurement surface and is not occluded by other scatterers and 0 otherwise, and the weighting factor ${\parallel \frac{\delta 𝐭}{\delta x}\times \frac{\delta 𝐭}{\delta \varphi }\parallel }_2$ is computed as described in Annex D.

The proposed definition of radio maps involves a non-coherent summation of path gains (43). As a result, it does not capture fast fading effects, which result from the constructive and destructive interference of multiple paths. Figure 29 demonstrates this by comparing three scenarios: a radio map generated using the radio map solver, a radio map created by non-coherently summing the path coefficients (4) calculated with the path solver (Section 3) by placing a single target at the center of each measurement cell, and a radio map produced by coherently summing the path coefficients (4) calculated with the path solver, also by placing a single target at the center of each measurement cell. The radio map obtained using the radio map solver aligns with the non-coherent summation of path coefficients calculated using the path solver, as anticipated. However, it does not account for fast fading effects, which requires the coherent summation of path coefficients.

The Sionna RT radio map solver supports antenna arrays of any size at the transmitter and allows for user-defined precoding vectors. However, it only supports synthetic arrays. In this setup, paths are generated from a single source located at the center of the array, and the channel response for each individual antenna is calculated by applying the appropriate phase shifts synthetically. This approximation holds true as long as the array size is small compared to the distance from the transmitter to the measurement plane and scatterers.

Using synthetic arrays enables efficient scaling with respect to the number of antennas. To see this, let $N_A^{\text{tx}}$ denote the number of antennas at the transmitter, and let ${𝐮}_P\in {ℂ}^{N_A^{\text{tx}}}$ be the precoding vector. For a given path connecting the source to the measurement plane, let $𝐄\in {ℂ}^2$ represent the electric field observed at the measurement plane. By employing a synthetic array, the corresponding electric fields for each antenna element are obtained by applying the array response vector (36), denoted as ${𝐮}_A^{\text{tx}}$, leading to

Here, ${𝐄}_A$ is a matrix of size $2\times N_A^{\text{tx}}$, where each column represents the electric fields radiated by each transmit antenna. The contribution of this path to the radio map is given by

This expression allows the path’s contribution to the radio map to be represented as

where $\alpha ={\left({𝐮}_A^{\text{tx}}\right)}^{𝖳}{𝐮}_P$ is a scalar, which depends on the departure direction of the path from the transmitter and thus varies for paths originating from the source. By precomputing $\alpha$ for each path, the effect of the transmitter array can be efficiently simulated by multiplying the electric field of each path by $\alpha$ whenever the path intersects the measurement surface. This operation does not depend on the size of the array.

Figure 30(b) illustrates the computation time needed to generate a radio map for a single source as the number of samples $N_S$ and the number of transmit antennas $N_A^{\text{tx}}$ are varied. The setup, depicted in Figure 30(a), employs the “simple street canyon” scene built into Sionna RT. The red ball indicates the source, positioned 20 meters above the ground. The radio map in Figure 30(a) is calculated using ${10}^8$ samples and a linear array of 512 transmit antennas. Notably, the radio map solver can compute radio maps with ${10}^9$ samples on a single NVIDIA RTX 4090 GPU. The computation time increases with both the number of samples and antennas, exceeding one second only for the most demanding computations. It is important to note that while the computation time increases with the number of antennas due to the precomputation of the array weighting factor $\alpha$, it does not increase the memory usage of the radio map solver.