3 Path Solver

Refer to caption — Figure 7: The path solver constructs propagation paths that connect a source and a target in a scene.

A path solver aims to determine a set of paths that connect two endpoints in a scene, as illustrated in Figure 7. From this set of paths, a time- or frequency-domain CIR can be computed. The path solver leverages SBR to efficiently compute paths in a scene. SBR operates as follows: First, rays²² 2 The processing of the rays is parallelized for accelerated computation. are traced from the source into directions specified by a lattice on the unit sphere (see Section 3.1.1). Then, an SBR loop iteratively tests the intersection of the rays with the scene geometry. At every intersection, an interaction type is selected from specular reflection, diffuse reflection, and refraction, and a new ray is spawned according to the selected interaction type. The method for choosing interaction types is detailed in Section 3.1.2. The loop terminates when no more rays remain active, meaning they have either bounced out of the scene, reached a predefined maximum number of bounces, or satisfied another termination criterion. Notably, at each interaction of a ray with a surface in the scene, only one scattered ray is spawned, as spawning multiple rays would result in an exponential increase in the number of rays, quickly overwhelming the available computational resources.

When a diffuse reflection occurs during the SBR loop, the solver checks whether there is a line-of-sight (LoS) to the target, meaning that the ray segment from the interaction point to the target point is not occluded. If such a connection exists, the path is deemed valid and is returned by the solver. This process is known as next-event estimation. Subsequently, a new ray is spawned in a random direction within the hemisphere defined by the surface normal of the primitive involved. This new ray may lead to the discovery of additional paths. Next-event estimation is only possible for diffuse reflections that scatter energy into all directions. For specular reflection and refraction, only a single ray is spawned in the only possible direction. Consequently, the probability of finding paths connecting the source to the target that end with a specular reflection or refraction by SBR is effectively zero, as depicted in Figure 8. This is because hitting a specific point (the target) with a specular reflection or refraction requires perfect alignment of the incident ray, which has zero probability in continuous space. Therefore, paths ending with a specular suffix cannot be determined solely by SBR. Especially, specular chains cannot be computed by SBR. Since specular chains typically carry a significant amount of the transported energy, SBR alone is insufficient for computing CIRs.

To compute paths ending with a specular suffix, including specular chains, the path solver uses an image method-based algorithm. However, a brute-force implementation of the image method would test all possible combinations of primitives that constitute the scene and hence is prohibitive even for moderately large scenes. Therefore, the path solver uses SBR as an heuristic to find path candidates. A path candidate consists of a sequence of primitives and interaction types ending with a specular reflection or refraction. As the path vertices determined during the SBR do not form a path that connects to the target, the image method is used to adjust the vertex locations to establish the connection. Intuitively, SBR pre-selects path candidates in the vicinity of the transmitter and receiver for the image method to compute corresponding valid paths to be returned by the solver. These steps of the Sionna RT path solver are depicted in Figure 9. Notably, only path candidates ending with a specular suffix (including specular chains) require further processing by the image method algorithm. Finally, for each valid path identified, the path solver computes the corresponding channel coefficients and propagation delays that characterize the EM wave propagation along that path.

Alternative methods for identifying paths that include specular chains or specular suffixes may involve detection spheres centered at the targets to capture paths ending with a specular suffix, and/or extending the direction of specular reflection (or refraction) to a lobe around the specular (or refracted) direction. Such methods only approximate the path vertex locations and the directions of the scattered rays, yielding inaccurate path coefficients and consequently resulting in incorrect CIRs.

The remainder of this section follows the architecture of the path solver shown in Figure 9. First, Section 3.1 presents the SBR-based candidate generator. Section 3.2 then describes the image method-based algorithm that processes candidates to find valid paths with specular suffixes. Finally, Section 3.3 details how the solver computes the complex-valued channel coefficients and propagation delays, which characterize the radio channel.

3.1 Generating Candidates by Shooting and Bouncing of Rays (SBR)

The candidate generator utilizes the SBR-based method as illustrated in Figure 10. We first provide an overview of the SBR procedure, followed by a detailed explanations. We consider a scene with a single source located at $\mathbf{s}\in\mathbb{R}^{3}$ and $N_{T}$ targets at $\mathbf{t}_{k}\in\mathbb{R}^{3}$ , where $k=1,\dots,N_{T}$ . Although the path solver built into Sionna RT can handle multiple sources, we focus on a single source for the sake of clarity. The process is identical for all sources. The paths traced by the candidate generator during the SBR procedure are referred to as samples in this section to avoid confusion with the valid paths returned to the user and the candidate paths forwarded to the image method-based algorithm. A sample therefore consists of a sequence of rays connected end-to-end, traced by the candidate generator. As we will see, a single sample can result in multiple valid or candidate paths, as illustrated in Figure 11.

The SBR procedure begins by spawning a user-specified number of rays, denoted by $N_{S}$ , from the source $\mathbf{s}$ in directions $\widehat{\mathbf{k}}_{n}^{(0)}$ on a Fibonacci lattice, where $n=1,\dots,N_{S}$ . Each ray starts a sample that is traced by the candidate generator. Details on the implemented Fibonacci lattice can be found in Section 3.1.1. The SBR loop then iterates until either all samples bounce out of the scene or reach a user-specified maximum number of bounces, denoted as the maximum depth $L$ . We denote the depth of the samples by $\ell=0,\dots,L$ . The first step of the SBR loop involves testing the intersection of the samples with the scene geometry. If a sample does not intersect any object, it has bounced out of the scene, and is deactivated. When a sample $n$ with depth $\ell$ intersects the scene, we denote by $\mathbf{v}_{n}^{(\ell)}\in\mathbb{R}^{3}$ the intersection point, by $\widehat{\mathbf{n}}_{n}^{(\ell)}\in\mathbb{R}^{3}\leavevmode\nobreak\ \left(% \left\lVert\widehat{\mathbf{n}}_{n}^{(\ell)}\right\rVert_{2}=1\right)$ the surface normal at the intersection point oriented towards the half-space containing the incident field, by $o_{n}^{(\ell)}\in\mathbb{N}_{0}$ the index of the intersected object, and by $m_{n}^{(\ell)}\in\mathbb{N}_{0}$ the index of the specific intersected primitive within that object. An interaction type, denoted by $\chi_{n}^{(\ell)}$ , is then sampled among specular reflection, diffuse reflection, and refraction, which determines how the sample will continue to propagate. Details on the sampling of the interaction type are provided in Section 3.1.2. The rest of the loop depends on the sampled interaction type.

Diffuse Reflection

As previously indicated, diffuse reflection is assumed to scatter energy in all directions of the half-space containing the incident field, which enables next event estimation. Therefore, if a diffuse reflection is sampled (i.e. $\chi_{n}^{(\ell)}={\cal S}$ ), the LoS between the intersection point and each target $\mathbf{t}_{1},\cdots,\mathbf{t}_{N_{T}}$ is tested. For each target $\mathbf{t}_{k}$ with an unobstructed LoS to the intersection point, a path connecting the source to the target

p_{n}^{(\ell)}=\left(\mathbf{s},\left(\mathbf{v}_{n}^{(1)},o_{n}^{(1)},m_{n}^{% (1)},\chi_{n}^{(1)}\right),\cdots,\left(\mathbf{v}_{n}^{(\ell)},o_{n}^{(\ell)}% ,m_{n}^{(\ell)},\chi_{n}^{(\ell)}\right),\mathbf{t}_{k}\right)

(11)

is recorded as a valid path to be returned to the user. Note that the path here is enriched with the intersected primitives and the interaction types. Finally, a new ray is spawned to continue the sample path, with its direction selected uniformly at random from the hemisphere defined by $\widehat{\mathbf{n}}_{n}^{(\ell)}$ , i.e. $\widehat{\mathbf{k}}_{n}^{(\ell+1)}\sim{\cal U}_{2\pi}(\widehat{\mathbf{n}}_{n% }^{(\ell)})$ , where ${\cal U}_{2\pi}(\widehat{\mathbf{n}}_{n}^{(\ell)})$ represents the uniform distribution on the hemisphere.

Specular Reflection and Refraction

When a specular reflection or refraction is sampled (i.e. $\chi_{n}^{(\ell)}={\cal R}$ or $\chi_{n}^{(\ell)}={\cal T}$ ), the path (11) is a valid path for target $\mathbf{t}_{k}$ , where $k\in\left\{1,\dots,N_{T}\right\}$ , only if

\frac{\mathbf{t}_{k}-\mathbf{v}_{n}^{(\ell)}}{\left\lVert\mathbf{t}_{k}-% \mathbf{v}_{n}^{(\ell)}\right\rVert_{2}}=\widehat{\mathbf{k}}_{n}^{(\ell+1)},

(12)

where

\widehat{\mathbf{k}}_{n}^{(\ell+1)}=\begin{cases}\widehat{\mathbf{k}}_{n}^{(% \ell)}-2\left(\left(\widehat{\mathbf{k}}_{n}^{(\ell)}\right)^{\mathsf{T}}% \widehat{\mathbf{n}}_{n}^{(\ell)}\right)\widehat{\mathbf{n}}_{n}^{(\ell)}&% \text{if }\chi_{n}^{(\ell)}={\cal R}\\ \widehat{\mathbf{k}}_{n}^{(\ell)}&\text{if }\chi_{n}^{(\ell)}={\cal T}\end{cases}

(13)

and where $\mathbf{u}^{\mathsf{T}}$ denotes the transpose of $\mathbf{u}$ . Since targets are modeled as points, condition (12) effectively occurs with zero probability, and, consequently, the path (11) is not a valid path. However, the specular suffix $p_{n}^{(\ell_{d}:)}$ , where $\ell_{d}$ is defined as in (10), can be further processed to determine vertices $\left(\widetilde{\mathbf{v}}_{n}^{(\ell_{d}+1)},\cdots,\widetilde{\mathbf{v}}_% {n}^{(\ell)}\right)$ such that

\widetilde{p}_{n}^{(\ell)}=\Bigg{(}\mathbf{s},\left(\mathbf{v}_{n}^{(1)},o_{n}% ^{(1)},m_{n}^{(1)},\chi_{n}^{(1)}\right),\cdots,\\ \underbrace{\left(\widetilde{\mathbf{v}}_{n}^{(\ell_{d}+1)},o_{n}^{(\ell_{d}+1% )},m_{n}^{(\ell_{d}+1)},\chi_{n}^{(\ell_{d}+1)}\right),\cdots,\left(\widetilde% {\mathbf{v}}_{n}^{(\ell)},o_{n}^{(\ell)},m_{n}^{(\ell)},\chi_{n}^{(\ell)}% \right),\mathbf{t}_{k}}_{\text{Processed specular suffix}}\Bigg{)}

(14)

is a valid path, as illustrated in Figure 12. This additional processing is performed using the image method in the algorithm described in Section 3.2.

Importantly, for a sequence of primitives and interaction types

\left(\mathbf{v}_{n}^{(\ell_{d})},\left(o_{n}^{(\ell_{d}+1)},m_{n}^{(\ell_{d}+% 1)},\chi_{n}^{(\ell_{d}+1)}\right),\cdots,\left(o_{n}^{(\ell)},m_{n}^{(\ell)},% \chi_{n}^{(\ell)}\right),\mathbf{t}_{i}\right)

(15)

where $\mathbf{v}_{n}^{(0)}=\mathbf{s}$ , there exists at most a single set of vertices $\widetilde{\mathbf{v}}_{n}^{(\ell_{d}+1)},\cdots,\widetilde{\mathbf{v}}_{n}^{(% \ell)}$ such that $\widetilde{p}_{n}^{(\ell)}$ is a valid path. Consequently, if $\ell_{d}=0$ , i.e. $p_{n}^{(\ell)}$ is a specular chain, then further processing by the image method of multiple identical candidates (15) would result in multiple identical valid paths, which would lead to erroneous CIRs. Therefore, if the sample $p_{n}^{(\ell)}$ is a specular chain, a hashing-based de-duplication mechanism, detailed in Section 3.1.3, ensures that specular chain candidates are only evaluated once. Notably, the path hashing mechanism discards redundant candidates “on-the-fly”, significantly reducing memory and compute requirements. However, when $\ell_{d}>0$ , then $\mathbf{v}_{n}^{(\ell_{d})}$ corresponds to a diffuse reflection point that is the result of a randomly sampled ray direction (either at the source or the last previous diffuse reflection point). Since the probability of two different paths having identical randomly sampled diffuse reflection points is effectively zero, path uniqueness is guaranteed in these cases without requiring the de-duplication mechanism. As shown in Figure 10, a LoS occlusion test to the target is used as a heuristic to discard path candidates that have little chance of resulting in a valid path after processing by the image method-based algorithm. Finally, a new ray is traced into the direction (13) to continue the sample path.

The candidate generator depicted in Figure 10 thus yields a set of valid paths connecting the source to all targets that do not require further processing (except for the computation of the corresponding path coefficients and delays) and end with a diffuse reflection. Additionally, it produces a set of candidate paths ending with specular suffixes, which require further processing to be either discarded or refined into valid paths. The path solver constructs paths that may include any combination and number of specular reflections, diffuse reflections, and refractions.

3.1.1 Sampling Initial Ray Directions

The SBR procedure starts by spawning rays from the source with directions obtained from a spherical Fibonacci lattice with $N_{S}$ points, which is defined as follows:

\widehat{\mathbf{k}}_{n}^{(0)}=\begin{bmatrix}\sin{\theta_{n}}\cos{\phi_{n}}\\ \sin{\theta_{n}}\sin{\phi_{n}}\\ \cos{\theta_{n}}\end{bmatrix}

(16)

where

	$\displaystyle\theta_{n}$	$\displaystyle=$	$\displaystyle\arccos{\frac{2n}{N_{S}}}$
	$\displaystyle\phi_{n}$	$\displaystyle=$	$\displaystyle 2\pi\frac{n}{\varphi}$

and $-\left\lfloor\frac{N_{S}}{2}\right\rfloor\leq n\leq\left\lceil\frac{N_{S}}{2}% \right\rceil-1$ , $\varphi=\frac{1+\sqrt{5}}{2}$ is the golden ratio, and $\left\lfloor\cdot\right\rfloor$ and $\left\lceil\cdot\right\rceil$ the floor and ceiling function, respectively. Note that when defining the Fibonacci lattice, the index $n$ takes values in the range $-\left\lfloor\frac{N_{S}}{2}\right\rfloor,\dots,\left\lceil\frac{N_{S}}{2}% \right\rceil-1$ instead of the usual $1,\dots,N_{S}$ . Sampling rays from the Fibonacci lattice results in ray tubes sharing approximately the same opening angle $\frac{4\pi}{N_{S}}$ .

The Fibonacci lattice was selected for its simple construction and its benefits over other lattices for numerical integration [12, 13]. Figure 13 shows the spherical Fibonacci lattice with 100, 1000, and 10000 points. The Sionna RT path solver uses $10^{6}$ points as the default value for $N_{S}$ , and higher values will ensure that more valid paths are found.

3.1.2 Sampling Interaction Types

Within the SBR loop, when a ray intersects the scene geometry, the algorithm generates at most one scattered ray. This constraint is necessary because generating multiple rays at each intersection, such as both reflected and refracted rays, would lead to an exponential increase in computational complexity and memory usage with the depth of the path. Instead, a single interaction type is randomly chosen based on a distribution:

{\cal Q}=\left\{q\left({\cal R}\right),q\left({\cal S}\right),q\left({\cal T}% \right)\right\},\qquad\text{where }q\left({\cal R}\right)+q\left({\cal S}% \right)+q\left({\cal T}\right)=1\text{ and }q\left(\cdot\right)\geq 0.

(17)

Here, $q\left(\chi\right)$ represents the probability of selecting the interaction type $\chi\in{\cal I}$ . Importantly, the choice of distribution ${\cal Q}$ does not affect the correctness of the calculated paths, which are solely dependent on the propagation environment. However, the choice of distribution ${\cal Q}$ significantly impacts the sample efficiency of the ray tracer, which is the number of samples $N_{S}$ needed to compute paths that connect a source to a target and capture the majority of the transported energy. Sampling only one interaction type per intersection may appear limiting, but this is offset by the large number of samples generated from each source. In most scenarios, a high sample count allows for identifying all significant paths.

For each sample $n$ at depth $\ell$ , the interaction type is sampled independently from other interactions and such that the probability of sampling each event type is proportional to the squared amplitudes of the reflection and refraction coefficients:

$\displaystyle q\left({\cal R}\right)$	$\displaystyle=R^{2}\frac{\left\|r_{\perp}\right\|^{2}+\left\|r_{\parallel}\right\|% ^{2}}{\left\|r_{\perp}\right\|^{2}+\left\|r_{\parallel}\right\|^{2}+\left\|t_{\perp% }\right\|^{2}+\left\|t_{\parallel}\right\|^{2}}$	(18)
$\displaystyle q\left({\cal S}\right)$	$\displaystyle=S^{2}\frac{\left\|r_{\perp}\right\|^{2}+\left\|r_{\parallel}\right\|% ^{2}}{\left\|r_{\perp}\right\|^{2}+\left\|r_{\parallel}\right\|^{2}+\left\|t_{\perp% }\right\|^{2}+\left\|t_{\parallel}\right\|^{2}}$
$\displaystyle q\left({\cal T}\right)$	$\displaystyle=\frac{\left\|t_{\perp}\right\|^{2}+\left\|t_{\parallel}\right\|^{2}}% {\left\|r_{\perp}\right\|^{2}+\left\|r_{\parallel}\right\|^{2}+\left\|t_{\perp}% \right\|^{2}+\left\|t_{\parallel}\right\|^{2}}$

Here, $r_{\perp}$ , $r_{\parallel}$ , $t_{\perp}$ , and $t_{\parallel}$ are the Fresnel coefficients (116), $S$ is the scattering coefficient (see Section A.8), and $R$ is the reflection reduction factor (119). This approach ensures that the interaction types that scatter more energy are sampled more frequently–a technique known as importance sampling [14]. This is similar to techniques for importance sampling of light transport paths as used in computer graphics.

To illustrate the importance of optimized sampling strategies, we compare uniform sampling, where each interaction type is selected with equal probability, to importance sampling using the scene shown in Figure 14(a). The scene consists of a metallic box containing a glass screen, with a source and target positioned on opposite sides. Diffuse reflection, specular reflection, and refraction are enabled, and the maximum depth $L$ is set to $5$ . As shown in Figure 14(b), using importance sampling leads to better sample efficiency when compared to uniform sampling, as convergence requires approximately $100\times$ less samples to capture most of the received energy. This demonstrates the benefit of selecting interactions according to the material properties of the scatterers.

3.1.3 Hashing of Specular Chains

As mentioned in Section 3.1 and illustrated in Figure 10, to avoid storing multiple identical specular chain candidates, a path hashing method is used to discard duplicates of previously found candidates “on-the-fly”. Recall that a specular chain candidate is defined by a sequence of interactions with the scene, each characterized by the intersected object $o_{n}^{(\ell)}\in\mathbb{N}_{0}$ , the corresponding intersected primitive $m_{n}^{(\ell)}\in\mathbb{N}_{0}$ , and the interaction type $\chi_{n}^{(\ell)}\in{\cal I}$ , as in (15) with $\ell_{d}=0$ , i.e.

c^{(\ell)}_{n}=\left(\mathbf{s},\left(o_{n}^{(1)},m_{n}^{(1)},\chi_{n}^{(1)}% \right),\cdots,\left(o_{n}^{(\ell)},m_{n}^{(\ell)},\chi_{n}^{(\ell)}\right)% \right).

(19)

For specular chain candidates, we do not need to store the path vertices. If the candidate leads to a valid path, these vertices will be computed. Invalid candidates will be discarded.

\texttt{hash}_{n}^{(0)}\leftarrow 0

\ell\leftarrow 1

3: while <Stop condition> do

\triangleright

SBR loop

\cdots

\texttt{inter\_hash}_{n}^{(\ell)}\leftarrow\textsc{hash\_intersection}(o_{n}^{% (\ell)},m_{n}^{(\ell)},\chi_{n}^{(\ell)})

\texttt{hash}_{n}^{(\ell)}\leftarrow\textsc{hash\_update}(\texttt{hash}_{n}^{(% \ell-1)},\texttt{inter\_hash}_{n}^{(\ell)})

\cdots

\ell\leftarrow\ell+1

9: end while

10:

11: procedure hash_intersection(

o, m, s

)

12:

h_{1}\leftarrow\textsc{cantor\_pair}(o,m)

\triangleright

Pair object and primitive

13:

h_{2}\leftarrow\textsc{cantor\_pair}(h_{1},s)

\triangleright

Pair with interaction type

14: return

h_{2}

15: end procedure

16:

17: procedure hash_update(

\texttt{curr\_hash},\texttt{inter\_hash}

)

18: return

(1373\cdot\texttt{curr\_hash}+\texttt{inter\_hash})

\triangleright

Update rolling hash

19: end procedure

20:

21: procedure cantor_pair(

x, y

)

22:

z\leftarrow\frac{x^{2}+x+2xy+3y+y^{2}}{2}

23: return

z

24: end procedure

Algorithm 1 Hashing of specular chains

Algorithm 1 demonstrates how path hashing works during ray tracing to avoid duplicate specular chain candidates. It starts by initializing a hash value $\texttt{hash}_{n}^{(0)}$ to $0$ for each sample $n=1,\dots,N_{S}$ being traced. Then, as the sample interacts with objects in the scene, the algorithm computes a hash value $\texttt{inter\_hash}_{n}^{(\ell)}$ for each intersection based on the object that was hit $o_{n}^{(\ell)}$ , the specific primitive of that object $m_{n}^{(\ell)}$ , and the type of interaction $\chi_{n}^{(\ell)}$ . The overall path hash $\texttt{hash}_{n}^{(\ell)}$ gets updated by combining it with the hash of each new intersection, creating a unique identifier for the entire sample up to that point.

Computing the hash after a given intersection is implemented by combining the three components of an intersection into a single hash value through a two-step process using the Cantor pairing function. The Cantor pairing function is a mapping between pairs of natural numbers and single natural numbers [15]. This bijective function guarantees that distinct input pairs always yield distinct output values, making it suitable for combining multiple values into a single hash. Hashing of the interaction first pairs the object and primitive indices, then combines the resulting value with the interaction type.

The hash_update function updates the sample hash value $\texttt{hash}_{n}^{(\ell)}$ through a polynomial rolling hash function [16] with a prime base of $1373$ . For each new intersection, it updates the running hash by multiplying the current hash value by the prime number and adding the new intersection hash. When a specular chain candidate is considered for a target $k\in\left\{1,\dots,N_{T}\right\}$ , the hash value of the candidate is paired with the target index $k$ to form a new hash value

\texttt{hash}_{n,k}^{(\ell)}=\textsc{cantor\_pair}(\texttt{hash}_{n}^{(\ell)},% k).

(20)

This ensures that the specular chain candidate is unique for each target.

Practical Considerations

To track information about previously encountered specular chains, an array of integers of size $N_{H}$ is allocated in memory and initialized to zero, referred to as the hash array. When considering a specular chain candidate for target $k$ with corresponding hash value $\texttt{hash}_{n,k}^{(\ell)}$ , the corresponding index $i_{n,k}^{(\ell)}$ in the hash array is computed as

i_{n,k}^{(\ell)}=\texttt{hash}_{n,k}^{(\ell)}\mod N_{H}.

(21)

If the hash array entry at index $i_{n,k}^{(\ell)}$ is greater than zero, it indicates that this path has already been encountered and should be discarded. Otherwise, the entry is incremented by one and the candidate is retained for further processing by the image method. It is important to note that the index $i_{n,k}^{(\ell)}$ is not guaranteed to be unique for different specular chains. When two different specular chains share the same index value despite being different, the path solver will discard one of them. Such an event is referred to as a collision, and can result in the discarding of valid paths. However, when two candidates are identical, their indices values will also be identical.

Valid paths and candidates are stored in a buffer of size $N_{B}$ . This buffer holds, for each path, the sequence of interaction types, vertices, intersected objects, primitives, and more. When a valid path is found or a candidate is identified as unique, indicated by its corresponding entry $i_{n,k}^{(\ell)}$ (21) in the hash array being set to zero, the path data are stored in the buffer. It is important to note that the index $i_{n,k}^{(\ell)}$ is solely used for indexing the hash array and not for storing path data in the buffer. Instead, path data are stored contiguously in the next available slot of the path buffer, as shown in Figure 15. If the buffer becomes full, any newly found paths or candidates are discarded. Discarding paths rather than overwriting existing ones is based on the reason that paths with lower depth, which are added to the buffer first as the SBR loop advances from $\ell=0$ to $\ell=L$ , typically carry more energy and should not be replaced by higher depth paths when the buffer becomes full. Consequently, the buffer size $N_{B}$ determines the maximum number of paths and candidates that can be stored.

Due to the substantial amount of data stored for each path, setting $N_{B}$ to excessively high values can lead to memory issues. This necessitates separating the hash array size $N_{H}$ from the buffer size $N_{B}$ . This is because a small hash array size can increase the probability of different hash values mapping to the same array index through the modulo operation (21), leading to collisions. To reduce the risk of collisions, $N_{H}$ is set to

N_{H}=\max\left(N_{B},10^{6}\right)

(22)

ensuring a minimum array size of $10^{6}$ elements. In the interface of the Sionna RT built-in path solver, $N_{B}$ corresponds to the samples_per_src parameter, and $N_{H}$ corresponds to max_num_paths_per_src.

3.1.4 Potential Improvements

Several potential improvements can be made to enhance the SBR-based candidate generator. To begin with, optimizing the sampling of ray directions can be achieved by incorporating importance sampling based on the transmit antenna pattern for initial ray spawning and the scattering pattern for diffusely reflected rays. This approach is motivated by the substantial sample efficiency gains shown in Section 3.1.2 through the use of importance sampling. By applying optimized sampling strategies for ray directions, further improvements may be realized as illustrated in Figure 16.

Moreover, the current method for handling multiple targets involves iterating over all targets during each SBR loop iteration. This approach may be optimized by selecting only targets near intersection points or the source, which will enhance scalability when dealing with a large number of targets. Figure 17(b) illustrates this, showing the compute time of the candidate generator and the average number of candidates per target found as the number of targets varies. In this experiment, the built-in Sionna RT scene “simple street canyon” depicted in Figure 17(a), was used. Targets were sampled uniformly at random on a plane parallel to the ground at an elevation of 1.5 meters. The maximum depth $L$ was set to 5, diffuse reflections were disabled, and the number of samples $N_{S}$ was set to $10^{6}$ . As observed, the compute time increases with the number of targets due to the iteration over all targets in each SBR loop iteration. However, the number of candidates per target remains stable, although with variations due to hash collisions. Note that, as discussed in Section 3.2, most candidates do not result in valid paths.

Scaling with the number of sources presents its own set of challenges. When $N_{S}$ samples are generated per source, the total number of samples produced by the candidate generator becomes $N_{S}\cdot N_{O}$ , where $N_{O}$ denotes the number of sources. This results in increased computation time as the number of sources grows, as illustrated in Figure 17(c). However, the number of candidates per source remains stable. In this experiment, a single target is positioned at the center of the scene at an elevation of 1.5 meters, while sources are randomly sampled on a plane parallel to the ground at an elevation of 70 meters, above all buildings. Recall that variations in the number of candidates per source may occur due to collisions. For a comprehensive view, Figure 17(d) presents results where the total number of samples is kept constant at $10^{6}$ , resulting in $N_{S}=\left\lfloor\frac{10^{6}}{N_{O}}\right\rfloor$ samples per source. In this scenario, the computation time increases at a much slower rate with the number of sources, but the number of candidates found per source decreases rapidly as the number of sources increases.

3.2 Image Method-based Candidate Processing

The image method processes path candidates to determine valid paths that are specular chains, or end with specular suffixes. Consider a candidate path

c_{n}^{(\ell)}=\left(\mathbf{s},\left(o_{n}^{(1)},m_{n}^{(1)},\chi_{n}^{(1)}% \right),\cdots,\underbrace{\left(o_{n}^{(\ell_{d}+1)},m_{n}^{(\ell_{d}+1)},% \chi_{n}^{(\ell_{d}+1)}\right),\cdots,\left(o_{n}^{(\ell)},m_{n}^{(\ell)},\chi% _{n}^{(\ell)}\right),\mathbf{t}}_{\text{Specular suffix}}\right)

(23)

where $\ell$ denotes the path depth, $\ell_{d}$ is the specular suffix index defined as in (10), and $\mathbf{t}$ is the target of the path. Recall that the sequence of interactions $\boldsymbol{\chi}^{(\ell)}=\left(\chi_{n}^{(1)},\cdots,\chi_{n}^{(\ell)}\right)$ was determined by the SBR-based candidate generator. For path candidates that end with a specular suffix but are not specular chains (i.e. $\ell_{d}>0$ ), the image method processes only the specular suffix. In these cases, the vertex of the last diffuse reflection $\mathbf{v}_{n}^{(\ell_{d})}$ serves as the source point for the image method calculations. Since the processing algorithm is identical for both specular chains and specular suffixes, we treat them equivalently in the following discussion. For simplicity of notation, we will refer to both types as specular chains and assume $\ell_{d}=0$ throughout the remainder of this section.

Figure 18 illustrates the image method process, which consists of two main steps. First, the method iteratively computes images of the source by reflecting it across each plane at which a specular reflection occurs. Since refracted paths are modeled without angular deflection, incorporating refraction into the image method is straightforward: refraction events may simply be skipped when computing the source images, as they do not alter the ray’s direction. Then, the algorithm backtracks from the target to determine the actual path vertices by intersecting lines between successive image sources.

Formally, the image method iteratively computes images of the source $\widetilde{\mathbf{s}}_{n}^{(i)}$ for $i=0,\dots,\ell$ , starting from the original source $\mathbf{s}$ :

\widetilde{\mathbf{s}}_{n}^{(i)}=\begin{cases}\mathbf{s}&\text{if }i=0\\ \widetilde{\mathbf{s}}_{n}^{(i-1)}-2\left(\widetilde{\mathbf{s}}_{n}^{(i-1)}-% \mathbf{p}_{n}^{(i)}\right)^{\mathsf{T}}\widehat{\mathbf{n}}_{n}^{(i)}&\text{% if }\chi_{n}^{(i)}={\cal R}\\ \widetilde{\mathbf{s}}_{n}^{(i-1)}&\text{if }\chi_{n}^{(i)}={\cal T}\end{cases}

(24)

Here, $\mathbf{p}_{n}^{(i)}$ represents any point on the plane containing the primitive $\left(o_{n}^{(i)},m_{n}^{(i)}\right)$ , and $\widehat{\mathbf{n}}_{n}^{(i)}$ denotes the surface’s normal vector. For specular reflections ( $\chi_{n}^{(i)}={\cal R}$ ), the image is computed by reflecting the previous image across the surface plane, while for refractions ( $\chi_{n}^{(i)}={\cal T}$ ), the image position remains unchanged.

During backtracking, the algorithm traces rays iteratively from the target $\mathbf{t}$ back to the source $\mathbf{s}$ by intersecting lines between successive image points. Specifically, a sequence of $\ell$ rays $\left(\widetilde{\mathbf{t}}_{n}^{(i)},\widehat{\mathbf{d}}_{n}^{(i)}\right)$ , $i=\ell,\dots,0$ , is traced according to

$\displaystyle\widetilde{\mathbf{t}}_{n}^{(i)}$	$\displaystyle=\begin{cases}\mathbf{t}&\text{if }i=\ell\\ \widetilde{\mathbf{v}}_{n}^{(i+1)}&\text{if }1\leq i<\ell\end{cases}$	(25)
$\displaystyle\widehat{\mathbf{d}}_{n}^{(i)}$	$\displaystyle=\frac{\widetilde{\mathbf{s}}_{n}^{(i)}-\widetilde{\mathbf{t}}_{n% }^{(i)}}{\left\lVert\widetilde{\mathbf{s}}_{n}^{(i)}-\widetilde{\mathbf{t}}_{n% }^{(i)}\right\rVert_{2}}$
$\displaystyle\widetilde{\mathbf{v}}_{n}^{(i)}$	$\displaystyle=\texttt{intersection}\left(\left(\widetilde{\mathbf{t}}_{n}^{(i)% },\widehat{\mathbf{d}}_{n}^{(i)}\right),\left(o_{n}^{(i)},m_{n}^{(i)}\right)% \right),\text{only if }i>0$

where $\widetilde{\mathbf{v}}_{n}^{(i)}$ is the intersection point between the ray $\left(\widetilde{\mathbf{t}}_{n}^{(i)},\widehat{\mathbf{d}}_{n}^{(i)}\right)$ and the primitive $\left(o_{n}^{(i)},m_{n}^{(i)}\right)$ , yielding the path vertex $\widetilde{\mathbf{v}}_{n}^{(i)}$ as shown in Figure 18(b). The path is considered valid only if two conditions are met: (i) each ray successfully intersects its corresponding primitive $\left(o_{n}^{(i)},m_{n}^{(i)}\right)$ , and (ii) the line segment between $\widetilde{\mathbf{t}}_{n}^{(i)}$ and $\widetilde{\mathbf{v}}_{n}^{(i)}$ is unobstructed by any other scene geometry. If either condition fails, the path will be discarded. The case $i=0$ corresponds to the final segment of the path, which connects to the source $\mathbf{s}$ . If the path is valid, then the returned path will be (14).

Typically, most of the candidates are discarded by the image method. This is illustrated in Figure 19, which shows the number of candidates discarded as a function of the path depth for a single source and target. The considered scene is shown in Figure 19(a). The number of samples $N_{S}$ is $10^{6}$ and diffuse reflection is disabled. As seen in Figure 19(b), most of the candidates are discarded by the image method because only specular chains are generated by the candidate generator, and most of them do not result in valid paths.

3.3 Channel Coefficients, Delays, and Doppler Shifts Computation

The final step of the path solver computes complex-valued channel coefficients and real-valued delays for each valid path (see Figure 9). For clarity, we will focus our discussion on a single valid path with index $n$ and depth $\ell$ , though the same procedure applies to all valid paths computed by both the SBR-based candidate generator and image method. We denote the path vertices as $\mathbf{v}_{n}^{(i)},\leavevmode\nobreak\ i=1,\dots,\ell$ , which may correspond to vertices previously computed by the image method (denoted earlier as $\widetilde{\mathbf{v}}_{n}^{(i)}$ ). The directions of propagation between vertices are represented by unit vectors $\widehat{\mathbf{k}}_{n}^{(i)},\leavevmode\nobreak\ i=0,\dots,\ell$ , defined as

\widehat{\mathbf{k}}_{n}^{(i)}=\frac{\mathbf{v}_{n}^{(i+1)}-\mathbf{v}_{n}^{(i% )}}{\left\lVert\mathbf{v}_{n}^{(i+1)}-\mathbf{v}_{n}^{(i)}\right\rVert_{2}},% \qquad i=0,\dots,\ell

(26)

where $\mathbf{v}_{n}^{(0)}=\mathbf{s}$ is the source and $\mathbf{v}_{n}^{(\ell+1)}=\mathbf{t}$ is the target. Note that these propagation directions may differ from those used during the initial SBR loop described in Section 3.1 due to the refinement by the image method.

Ensuring Energy Conservation

As explained in Section 3.1.2, at each interaction point, an interaction type is randomly chosen based on the probability distribution ${\cal Q}$ (17). Initially, the energy emitted by the source is distributed across $N_{S}$ samples according to the transmit antenna pattern (see Section A.3). Since only one interaction type is selected at each interaction point, ray tubes that do not match the chosen interaction type are discarded, as illustrated in Figure 20. The remaining ray tubes must therefore compensate for the discarded ones to maintain energy conservation. Specular reflection and refraction on planar surfaces only change the direction of the ray tube’s propagation. Given that there is a finite number of specular chains with a specific depth, the channel gain (5) observed at the target due to specular chains is accurate if all specular chains are identified by the path solver, which is assumed to be the case. This assumption will be reasonable if the number of samples $N_{S}$ is sufficiently large.

In the context of diffuse reflections, an incident ray tube scatters an amount of energy that is determined by its footprint on the diffusely reflective surface (see Section A.8). Due to the random sampling of interaction types, only a portion of the ray tubes that should reach a diffusely reflective surface are actually traced (see Figure 20), which requires compensating for the untraced ones. For a given path $p^{(\ell)}$ that includes a diffuse reflection and has a specular suffix index $\ell_{d}$ (10), we denote by $\Pr\left(\boldsymbol{\chi}_{n}^{(\ell_{d})}\right)=\prod_{\ell^{\prime}=1}^{% \ell_{d}}q\left(\chi_{n}^{(\ell^{\prime})}\right)$ the probability of sampling the prefix leading to the last diffuse reflection point, where $q\left(\cdot\right)$ is defined in (17). Consequently, on average, only a fraction $\Pr\left(\boldsymbol{\chi}_{n}^{(\ell_{d})}\right)$ of the ray tubes that should reach the diffusely reflective surface are actually traced. To account for this, the electric field is scaled by $\nicefrac{{1}}{{\sqrt{\Pr\left(\boldsymbol{\chi}_{n}^{(\ell_{d})}\right)}}}$ at the point of diffuse reflection.

Algorithm 2 Electric field computation

\mathbf{E}_{n}^{(0)}\leftarrow\textsc{tx\_pattern}(\widehat{\mathbf{k}}^{(0)}_% {n})

\triangleright

Electric field vector

\tau_{n}^{(0)}\leftarrow 0

\triangleright

Cumulative delay

r_{n}^{(0)}\leftarrow 0

\triangleright

Cumulative ray tube length

\gamma_{n}^{(0)}\leftarrow 1

\triangleright

Cumulative path probability

\ell^{\prime}\leftarrow 1

6: while not <Stop condition> do

\mathbf{T}_{n}^{(\ell^{\prime})}\leftarrow\textsc{Evaluate\_Material}(o_{n}^{(% \ell^{\prime})},\chi_{n}^{(\ell^{\prime})},\mathbf{E}_{n}^{(\ell^{\prime}-1)},% \widehat{\mathbf{k}}_{n}^{(\ell^{\prime}-1)},\widehat{\mathbf{k}}_{n}^{(\ell^{% \prime})})

\mathbf{E}_{n}^{(\ell^{\prime})}\leftarrow\mathbf{T}_{n}^{(\ell^{\prime})}% \mathbf{E}_{n}^{(\ell^{\prime}-1)}

\gamma_{n}^{(\ell^{\prime})}\leftarrow\gamma_{n}^{(\ell^{\prime}-1)}\cdot q_{n% }^{(\ell^{\prime})}\left(\chi_{n}^{(\ell^{\prime})}\right)

10: if

\chi_{n}^{(\ell^{\prime})}={\cal S}

then

11:

\mathbf{E}_{n}^{(\ell^{\prime})}\leftarrow\frac{\mathbf{E}_{n}^{(\ell^{\prime}% )}}{\sqrt{\gamma_{n}^{(\ell^{\prime})}}}

12:

\gamma_{n}^{(\ell^{\prime})}\leftarrow 1

13:

r_{n}^{(\ell^{\prime}-1)}\leftarrow 0

14: end if

15:

r_{n}^{(\ell^{\prime})}\leftarrow r_{n}^{(\ell^{\prime}-1)}+\left\lVert\mathbf% {v}_{n}^{(\ell^{\prime})}-\mathbf{v}_{n}^{(\ell^{\prime}-1)}\right\rVert_{2}

16:

\tau_{n}^{(\ell^{\prime})}\leftarrow\tau_{n}^{(\ell^{\prime}-1)}+\frac{\left% \lVert\mathbf{v}_{n}^{(\ell^{\prime})}-\mathbf{v}_{n}^{(\ell^{\prime}-1)}% \right\rVert_{2}}{c}

17:

\ell^{\prime}\leftarrow\ell^{\prime}+1

18: end while

19:

\tau_{n}\leftarrow\tau_{n}^{(\ell)}+\frac{\left\lVert\mathbf{t}-\mathbf{v}_{n}% ^{(\ell)}\right\rVert_{2}}{c}

20: if

\chi_{n}^{(\ell)}={\cal S}

then

21:

r_{n}^{(\ell-1)}\leftarrow 0

22: end if

23:

r_{n}^{(\ell)}\leftarrow r_{n}^{(\ell-1)}+\left\lVert\mathbf{t}-\mathbf{v}_{n}% ^{(\ell)}\right\rVert_{2}

24:

\mathbf{E}_{n}^{(r)}\leftarrow\textsc{rx\_pattern}(-\widehat{\mathbf{k}}_{n}^{% (\ell)})

25:

a_{n}\leftarrow\frac{\lambda}{4\pi r_{n}^{(\ell)}}\left(\left(\mathbf{E}_{n}^{% (\ell)}\right)^{\mathsf{T}}\mathbf{E}_{n}^{(r)}\right)

26: return

a_{n},\leavevmode\nobreak\ \tau_{n}

Path Coefficient and Delay Computation

Algorithm 2 details the computation of the electric field propagation along a valid path, leveraging the EM background presented in Appendix A. The algorithm begins by initializing the electric field vector $\mathbf{E}_{n}^{(0)}$ using the transmit antenna’s pattern evaluated in the initial propagation direction $\widehat{\mathbf{k}}_{n}^{(0)}$ . Note that the path solver supports dual polarization by considering two independent electric field vectors for each traced path. The variable $\tau_{n}^{(0)}$ tracks the cumulative path delay and is initialized to zero, and the variable $r_{n}^{(0)}$ tracks the ray tube length and is also initialized to zero. The variable $\gamma_{n}^{(0)}$ is used to track the cumulative path probability and is initialized to one (see previous paragraph). For each interaction point along the path, the electric field is updated by a transformation matrix $\mathbf{T}_{n}^{(\ell^{\prime})}$ , which is computed based on the intersected object $o_{n}^{(\ell^{\prime})}$ , interaction type $\chi_{n}^{(\ell^{\prime})}$ , incident field $\mathbf{E}_{n}^{(\ell^{\prime}-1)}$ , and both incident $\widehat{\mathbf{k}}_{n}^{(\ell^{\prime}-1)}$ and scattered $\widehat{\mathbf{k}}_{n}^{(\ell^{\prime})}$ directions. The cumulative path probability $\gamma_{n}^{(\ell^{\prime})}$ is updated by multiplying it with the interaction probability $q_{n}^{(\ell^{\prime})}\left(\chi_{n}^{(\ell^{\prime})}\right)$ . When a diffuse reflection occurs, the electric field is normalized by $\sqrt{\gamma_{n}^{(\ell^{\prime})}}$ to maintain energy conservation and $\gamma_{n}^{(\ell^{\prime})}$ is reset to one and the ray tube length $r_{n}^{(\ell^{\prime})}$ to zero. The ray tube length $r_{n}^{(\ell^{\prime})}$ is updated by adding the distance between the current and previous vertex, and the delay $\tau_{n}^{(\ell^{\prime})}$ is updated by adding the time it takes to travel that distance, where $c$ is the speed of light. After replaying all interactions, the final path segment to the target is incorporated into the delay and ray tube length. Finally, the receive antenna pattern is evaluated for the incoming direction, and the channel coefficient is computed following Section A.6.

Remark

The transformation matrices $\mathbf{T}_{n}^{(\ell^{\prime})}$ and the electric field vectors $\mathbf{E}_{n}^{(\ell^{\prime})}$ are complex-valued, but are implemented using their real-valued representations defined as

\begin{bmatrix}\Re\left(\mathbf{T}_{n}^{(\ell^{\prime})}\right)&-\Im\left(% \mathbf{T}_{n}^{(\ell^{\prime})}\right)\\ \Im\left(\mathbf{T}_{n}^{(\ell^{\prime})}\right)&\Re\left(\mathbf{T}_{n}^{(% \ell^{\prime})}\right)\end{bmatrix}

(27)

for matrices and

\begin{bmatrix}\Re\left(\mathbf{E}_{n}^{(\ell^{\prime})}\right)\\ \Im\left(\mathbf{E}_{n}^{(\ell^{\prime})}\right)\end{bmatrix}

(28)

for vectors, where $\Re{\left(\cdot\right)}$ and $\Im{\left(\cdot\right)}$ denote the real and imaginary parts, respectively.

3.3.1 Handling Multiple Antennas at the Transmitter and Receiver

As explained in the Introduction (Section 1), Sionna RT supports antenna arrays of any size at both the transmitter and receiver, utilizing either synthetic or non-synthetic arrays. In the case of non-synthetic arrays, each transmit antenna is treated as a source and each receive antenna as a target, meaning paths are traced from every transmit antenna, and the solver calculates paths for each transmit-receive antenna pair. For synthetic arrays, paths originate from a single “virtual” source at the center of the transmitter array, and a single “virtual” target is considered at the center of every receiver array. Channel coefficients are then computed by applying phase shifts synthetically. This method allows for efficient scaling with the number of antennas, but it is an approximation that is valid only when the transmitter and receiver array sizes are small relative to the distance from the radio devices to each other and to the scatterers. Formally, for synthetic arrays, the channel matrix for the $n$ -th path is expressed as:

\mathbf{H}_{\text{array},n}=a_{n}\mathbf{u}_{A}^{\text{rx}}\left(\mathbf{u}_{A% }^{\text{tx}}\right)^{\mathsf{T}}

(29)

where $a_{n}$ is the channel coefficient for the $n$ -th path, calculated using Algorithm 2. The array response vectors for the transmit ( $\mathbf{u}_{A}^{\text{tx}}$ ) and receive ( $\mathbf{u}_{A}^{\text{rx}}$ ) antenna arrays are defined as:

	$\displaystyle\mathbf{u}_{A}^{\text{tx}}$	$\displaystyle=\left[e^{j\frac{2\pi}{\lambda}\left(\widehat{\mathbf{k}}_{n}^{(0% )}\right)^{\mathsf{T}}\mathbf{d}_{1}^{\text{tx}}},\cdots,e^{j\frac{2\pi}{% \lambda}\left(\widehat{\mathbf{k}}_{n}^{(0)}\right)^{\mathsf{T}}\mathbf{d}_{N_% {A}^{\text{tx}}}^{\text{tx}}}\right]^{\mathsf{T}},$		(30)
	$\displaystyle\mathbf{u}_{A}^{\text{rx}}$	$\displaystyle=\left[e^{j\frac{2\pi}{\lambda}\left(-\widehat{\mathbf{k}}_{n}^{(% \ell)}\right)^{\mathsf{T}}\mathbf{d}_{1}^{\text{rx}}},\cdots,e^{j\frac{2\pi}{% \lambda}\left(-\widehat{\mathbf{k}}_{n}^{(\ell)}\right)^{\mathsf{T}}\mathbf{d}% _{N_{A}^{\text{rx}}}^{\text{rx}}}\right]^{\mathsf{T}}.$		(30)

Here, $\widehat{\mathbf{k}}_{n}^{(0)}$ ( $\widehat{\mathbf{k}}_{n}^{(\ell)}$ ) represents the departure (arrival) direction of the $n$ -th path, $N_{A}^{\text{tx}}$ ( $N_{A}^{\text{rx}}$ ) denotes the number of antennas in the transmit (receive) array, and $\mathbf{d}_{i}^{\text{tx}}$ ( $\mathbf{d}_{i}^{\text{rx}}$ ) indicates the relative position of the $i$ -th antenna element to the virtual source (target) in the transmit (receive) array.

3.3.2 Doppler Shifts

The impact of moving scene objects on the paths can be simulated in two ways. The first approach involves moving objects in small incremental steps along their trajectories and recomputing the propagation paths at each step. While this method provides the highest accuracy, it is computationally expensive. However, when trajectory lengths are small (not exceeding a few wavelengths), a more efficient approach may be used. In this alternative method, we assign a velocity vector $\boldsymbol{\varv}\left(o\right)\in\mathbb{R}^{3}$ to each object $o$ and compute the time evolution by accumulating Doppler shifts along the propagation paths.

Figure 21 illustrates this approach. When replaying the valid paths using Algorithm 2, the Doppler shift for each path $n$ is computed by accumulating the contributions from individual objects:

\nu_{n}^{(\ell^{\prime})}=\begin{cases}0&\text{if }\ell^{\prime}=0\\ \nu_{n}^{(\ell^{\prime}-1)}+\frac{\boldsymbol{\varv}\left(o_{n}^{(\ell^{\prime% })}\right)^{\mathsf{T}}\left(\widehat{\mathbf{k}}_{n}^{(\ell^{\prime})}-% \widehat{\mathbf{k}}_{n}^{(\ell^{\prime}-1)}\right)}{\lambda}&\text{if }1\leq% \ell^{\prime}<\ell\\ \end{cases}

(31)

where $\lambda$ is the wavelength. This results in the total Doppler shift $\nu_{n}^{(\ell)}$ caused by object movements along the path. The final Doppler shift observed at the target combines the accumulated shift from moving objects with the shifts caused by source and target motion:

\nu_{n}=\nu_{n}^{(\ell)}+\frac{\boldsymbol{\varv}\left(\mathbf{s}\right)^{% \mathsf{T}}\widehat{\mathbf{k}}_{n}^{(0)}}{\lambda}-\frac{\boldsymbol{\varv}% \left(\mathbf{t}\right)^{\mathsf{T}}\widehat{\mathbf{k}}_{n}^{(\ell)}}{\lambda}.

(32)

where $\boldsymbol{\varv}\left(\mathbf{s}\right)\in\mathbb{R}^{3}$ and $\boldsymbol{\varv}\left(\mathbf{t}\right)\in\mathbb{R}^{3}$ denote the velocity vectors of the source and target, respectively.

NVIDIA Sionna RT