# Calibrating Wireless Ray Tracing for Digital Twinning using Local Phase
Error Estimates

Clement Ruah, Osvaldo Simeone, Jakob Hoydis, and Bashir Al-Hashimi,  Clement
Ruah and Osvaldo Simeone are with the King’s Communications, Learning &
Information Processing (KCLIP) lab within the Centre for Intelligent
Information Processing Systems (CIIPS), Department of Engineering, King’s
College London, WC2R 2LS London, U.K. (e-mail<EMAIL_ADDRESS>osvaldo.simeone@kcl.ac.uk). Jakob Hoydis is with NVIDIA Corporation. Bashir M.
Al-Hashimi is with the Department of Engineering, King’s College London, WC2R
2LS London, U.K. The work of C. Ruah was supported by the Faculty of Natural,
Mathematical, and Engineering Sciences, King’s College London. The work of O.
Simeone and of J. Hoydis was partially supported by the European Union’s
Horizon Europe project CENTRIC (101096379). O. Simeone was also supported by
the Open Fellowships of the EPSRC (EP/W024101/1) by the EPSRC project
(EP/X011852/1), and by Project REASON, a UK Government funded project under
the Future Open Networks Research Challenge (FONRC) sponsored by the
Department of Science Innovation and Technology (DSIT).

###### Abstract

Embodying the principle of simulation intelligence, digital twin (DT) systems
construct and maintain a high-fidelity virtual model of a physical system.
This paper focuses on ray tracing (RT), which is widely seen as an enabling
technology for DTs of the radio access network (RAN) segment of next-
generation disaggregated wireless systems. RT makes it possible to simulate
channel conditions, enabling data augmentation and prediction-based
transmission. However, the effectiveness of RT hinges on the adaptation of the
electromagnetic properties assumed by the RT to actual channel conditions, a
process known as calibration. The main challenge of RT calibration is the fact
that small discrepancies in the geometric model fed to the RT software hinder
the accuracy of the predicted phases of the simulated propagation paths.
Existing solutions to this problem either rely on the channel power profile,
hence disregarding phase information, or they operate on the channel responses
by assuming the simulated phases to be sufficiently accurate for calibration.
This paper proposes a novel channel response-based scheme that, unlike the
state of the art, estimates and compensates for the phase errors in the RT-
generated channel responses. The proposed approach builds on the variational
expectation maximization algorithm with a flexible choice of the prior phase-
error distribution that bridges between a deterministic model with no phase
errors and a stochastic model with uniform phase errors. The algorithm is
computationally efficient, and is demonstrated, by leveraging the open-source
differentiable RT software available within the Sionna library, to outperform
existing methods in terms of the accuracy of RT predictions.

##  I Introduction

###  I-A Simulation Intelligence for Wireless Networks

A general trend in the sciences and engineering is the use of the
unprecedented processing power afforded by modern computing clusters to
implement _physics-driven simulators_ of the real world, with the aim of
facilitating scientific discovery, providing design insights, and reducing the
need for data collection. An important embodiment of this trend in
engineering, referred to as _simulation intelligence_ [1], is the adoption of
_digital twin_ (DT) systems. A DT software constructs and maintains a high-
fidelity virtual model of a physical system through two-way communications
between the real and the virtual worlds [2, 3, 4, 5, 6]. DTs provide a
_virtual sandbox_ in which operational policies can be designed and tested
before being deployed to the physical twin (PT).

DTs offer a natural framework for the deployment of software-based
disaggregated wireless networks that follow the open radio access network
(RAN) paradigm [7, 8, 9, 10, 11]. Among the key advantages of DTs in this
context is their capacity to synthesize data that can be leveraged to train,
as well as to feed at run time, artificial intelligence (AI) models. For
example, in the O-RAN architecture, real-time RAN intelligent controllers run
AI models that need to be tailored to the operating conditions of the
underlying physical network [7]. Data are generated at the DT by simulating
relevant conditions in the RAN. For instance, a DT may be used to estimate key
performance indicators (KPIs) on the RAN, such as throughput, on the basis of
information about the topology of the network and traffic conditions [12]. The
benefits of DT-generated data hinge on the adherence of the virtual model to
the physical system, and hence a key problem is the _calibration_ of the DT
using real-world data from the physical twin [13].

###  I-B Ray Tracing at the Digital Twin

_Ray tracing_ (RT) is widely seen as an enabling technology for DTs of the RAN
[14, 15, 12, 16, 17, 11, 18]. As illustrated in Fig. 1, given a three-
dimensional model specifying the shapes, positions, and electromagnetic
properties of the objects in the propagation environment, the RT simulates
multiple propagation paths between any given transmitter (Tx) and receiver
(Rx) through specular and diffuse reflection, refraction, and diffraction [19,
Ch. 4], [20, 21, 22]. By leveraging the predicted path components, a DT can
simulate the propagation conditions for a given network deployment. The
synthesized channels can then be used to train AI models that carry out tasks
such as beamforming [15, 16, 11, 23], user positioning [18], and channel
charting [24, 25] by simulating relevant traffic conditions and user
distributions. Given specific user locations, they may be also used to
directly optimize transmission strategies with minimal requirements in terms
of pilot overhead [16, 26].

Figure 1:  Taking as input the geometric properties of the scene, the
electromagnetic material parameters $\theta$, and the coordinates $c$ of
transmitter (Tx) and receiver (Rx), the ray tracer (RT) produces the features
$R(c|\theta)$ of a number $P$ of propagation paths. Based on this information,
the DT can obtain a model $H(c|\theta)$ of the channel conditions between
transmitter and receiver. To keep a faithful representation of its physical
twin (PT), during a calibration phase, the DT compares its model predictions
$H(c|\theta)$ to measured channel realizations
$\mathcal{D}=\\{H_{n},c_{n}\\}_{n=1}^{N}$ in order to optimize the material
parameters $\theta$. An accurate estimate of the ground-truth channel
conditions can then be used to monitor the PT, simulate alternate scenarios
(counterfactual analysis), or to optimize PT operations such as beamforming,
with reduced requirements on pilot transmissions and channel measurements
[16]. (The map and the 3D model in the the figure are built using geospatial
data from the OpenStreetMap database [27].)

The fidelity of the virtual simulation depends on the precision of the
geometric and electromagnetic properties, including permittivity, conductivity
and permeability [19, 28], fed to the DT as inputs. While specialized reports
[28] can provide approximate values for the electromagnetic properties of
generic materials, a more precise and bespoke selection is generally necessary
to ensure sufficiently reliable simulations. This optimization requires the
estimation of the material parameters of each object in the scene based on
measurements obtained from the ground-truth deployment scenario, in a process
known as _calibration_ [29, 30, 31, 32, 33, 34, 13].

###  I-C Calibration of Ray Tracing

The geometric properties of the deployment scenario of interest can be
generally determined up to a certain degree of precision based on in-situ
measurements, e.g., using LiDAR [35, 5, 36] or neural surface reconstruction
[37], or based-on geo-spatial data, e.g., via the Open Street Maps database
[27]. However, in practical scenarios, it is not uncommon to observe small
discrepancies, of the scale of a fraction of the carrier wavelength, between
the ground-truth geometric properties and their virtual counterpart. These
discrepancies hinder the accuracy of the predicted phases of the simulated
propagation paths, which are usually deemed to be unreliable [19, Ch. 7.5].

We distinguish two classes of existing calibration procedures, which make
distinct assumptions on the _phase error_ caused by imprecise geometric
information.

  * •

_Power profile-based_ schemes aim to find the material parameters that yield
the best match between the power-delay – or more generally power-angle-delay –
profiles simulated by the RT and the corresponding power profiles measured
from the real environment [29, 30, 31, 32, 33, 34]. Accordingly, power
profile-based schemes disregard phase information, retaining solely the
amplitudes of the measured and simulated channels.

  * •

_Channel response-based_ schemes find the material parameters that yield the
best match between simulated and recorded channel responses [13]. The state of
the art work [13] compares RT-generated and measured channel responses by
relying on the phase provided by the RT to be sufficiently accurate to enable
calibration. We refer to this class of channel response-based methods as
_phase error-oblivious_. Such an approach adopts an opposite modelling
assumption as compared to power profile-based methods, disregarding the
presence of phase errors in the RT simulation.

This paper proposes a novel channel response-based scheme that, unlike the
state of the art, compensates for the phase errors in the RT-generated channel
responses prior to carrying out a comparison with measured channel responses
during RT calibration. The main idea is to use measured channel responses to
estimate the phase error locally for each path and for each measurement. To
illustrate the potential benefits of the approach, we provide the following
toy example, which is depicted in Fig. 2.

(a) Ground-truth scenario

(b) Simulated scenario

(c) Ground-truth power-delay profiles

(d) Simulated power-delay profiles

(e) Calibration under high-bandwidth conditions

(f) Calibration under low-bandwidth conditions

Figure 2:  Toy example that illustrates the limitations of existing power
profile-based RT calibration methods, as well as of channel response-based
schemes that disregard path phase errors. In this example, the two propagation
paths between the transmitter (Tx) and the receiver (Rx) reflect on surfaces
with ground-truth reflectance $r^{*}$ in (a), and calibrated reflectance
$\hat{r}$ in the simulated scenario (b). Though the paths interfere
constructively at the receiver under ground-truth conditions, they are
predicted to interfere destructively in the simulated scenario due to an
inaccuracy $\Delta d$ of the order of the carrier’s wavelength $\lambda$ in
the geometric model. Due to this model error, phase error-oblivious schemes
provide inaccurate calibration results in both high and low-bandwidth regimes
(e)-(f). Furthermore, the model error causes the RT-generated power-delay
profiles in the low-bandwidth regime (d) to diverge from the ground-truth
power-delay profile (c), causing inaccurate calibration also for power
profile-based schemes.

In this example, as shown in Fig. 2(a) and Fig. 2(b), the goal of calibration
is to produce an estimate $\hat{r}$ of the non-negative ground-truth
reflectance value $r^{*}$, which represents the fraction of an incoming ray’s
energy that is not absorbed by the given surfaces upon reflection. As also
illustrated in the figures, the ground-truth propagation scenario and the
simulated scenario differ by a small amount of the order of the carrier
frequency’s wavelength $\lambda$. Owing to this small discrepancy, as seen in
Fig. 2(c) and Fig. 2(d), the two propagation paths, having different delays
$\tau_{1}$ and $\tau_{2}$, reach the receiver with the same phase in the
ground-truth scenario, while they have opposite phases in the simulated paths.

In an ideal situation in which the bandwidth, and hence the temporal
resolution, is large enough to separate the contributions of each individual
path, the power-delay profiles of the received signal (Fig. 2(c)) and of the
simulated paths (Fig. 2(d)) match perfectly by choosing the correct
reflectance value $\hat{r}=r^{*}$. Accordingly, as seen in Fig. 2(e), the
likelihood function maximized by a power profile-based scheme peaks at the
correct value $\hat{r}=r^{*}$.

In a potentially more relevant low-bandwidth regime, the temporal resolution
is too low to separate the two propagation paths, and the measured power
profile in Fig. 2(c) cannot be reproduced by the power profile generated by
the RT depicted in Fig. 2(d). In fact, as mentioned, the power profile
produced by the RT is characterized by destructive interference between the
two paths irrespective of the reflectance value $\hat{r}$, causing a mismatch
with the constructive interference exhibited by the measured power profile.
This results in an incorrectly estimated reflectance value $\hat{r}$ (Fig.
2(f)).

In both high and low-bandwidth regimes, a phase error-oblivious channel
response-based scheme would be unable to match measured and received responses
due to its reliance on the opposite phases for the two paths predicted by the
RT. As a result, this approach would return an incorrect value for the
estimated reflectance, as represented by the incorrect peak value for the
likelihood in Fig. 2(e) and Fig. 2(f). In contrast, as we will demonstrate in
this work, a channel response-based strategy that attempts to estimate and
mitigate simulated phase errors can provide well-calibrated electromagnetic
properties irrespective of the resolution of the physical system.

###  I-D Main Contributions

In this work, we present a novel channel response-based calibration scheme for
RT that estimates and mitigates per-path phase errors. Our specific
contributions are as follows.

  * •

We introduce a novel _phase error-aware_ calibration strategy for RT. The
proposed approach builds on the _variational expectation maximization_ (VEM)
algorithm [38] with a tailored choice of the variational distribution on the
phase errors. The algorithm is computationally efficient encompassing an
expectation step in closed form and a maximization step via gradient descent.

  * •

Leveraging the open-source differentiable RT made available through the Sionna
library [39, 13], we present a variety of experiments, evaluating calibration
accuracy in terms of the estimated electromagnetic parameters, as well as of
predicted power maps. The experiments confirm the significant advantages of
the proposed phase error-aware calibration strategy as compared to existing
power profile-based schemes and phase error-oblivious approaches.

The rest of the paper is organized as follows. Sec. II introduces the
considered communication system and the phase-error RT channel model. In Sec.
III, we describe the implemented calibration baselines for the phase error-
oblivious (Sec. III-A) and the power profile-based (Sec. III-B) schemes. Sec.
IV covers our proposed phase error-aware calibration scheme, which is compared
to the previously introduced baselines in Sec. V, where we study the accuracy
of the calibrated material parameters and predicted signal powers for a
synthetic scenario. Finally, Sec. VI concludes the paper.

##  II Setting and Problem Definition

In this work, we study the problem of calibrating a DT-based RT using channel
frequency response observations. This section first describes the studied
communication system and the RT output before introducing a stochastic phase-
error channel model that accounts for the errors of the RT-predicted phases.
We end this section by specifying the available data and the calibration
objective.

###  II-A Setting

Consider a transmitter (Tx) located at coordinates
$c^{\mathrm{Tx}}\in\mathbbm{R}^{3}$ and a receiver (Rx) located at coordinates
$c^{\mathrm{Rx}}\in\mathbbm{R}^{3}$ in a three-dimensional Cartesian space.
The transmitter and receiver are equipped with antenna arrays of
$N^{\mathrm{Tx}}$ and $N^{\mathrm{Rx}}$ antennas respectively. Using the paths
computed by the RT, the DT models the channel over a bandwidth
$B=[f^{\mathrm{min}},f^{\mathrm{max}}]$ with end frequencies
$f^{\mathrm{max}}>f^{\mathrm{min}}>0$ and central carrier frequency
$f^{c}=(f^{\mathrm{min}}+f^{\mathrm{max}})/2$ [Hz]. Assuming multi-carrier
transmission, the set of $S$ subcarrier frequencies is
$\\{f_{s}\\}_{s=1}^{S}$, where $f_{s}=f^{\mathrm{min}}+(s-1)\Delta f$ is the
frequency of the $s\mbox{-th}$ subcarrier; $\Delta f$ [Hz] is the subcarrier
spacing; and $S=\left\lfloor(f^{\mathrm{max}}-f^{\mathrm{min}})/\Delta
f\right\rfloor$ is the total number of subcarriers.

As illustrated in Fig. 1, the RT takes as input the geometric properties,
specifying the shapes and positions of the objects in the scene, the
coordinates $c=(c^{\mathrm{Rx}},c^{\mathrm{Tx}})$ of the transmitter and
receiver devices, and the materials properties, denoted as vector $\theta$,
which includes the permittivity, conductivity, permeability, scattering
coefficient and cross-polarization discrimination of the materials involved
[28, 21, 22]. Additional RT inputs not depicted in Fig. 1 are the antenna
patterns of the devices [40, Ch. 2.2], as well as the choice of a
parameterized scattering model [21]. Based on this information, the RT outputs
the parameters of a number $P^{c}$ of realizable propagation paths through
specular and diffuse reflection, refraction, and diffraction [19, Ch. 4], [20,
21, 22]. Note that the RT is only able to consider a finite number of samples
among the potentially infinite set diffusely scattered paths, whereas the
finite set of specular components is treated exhaustively up to a pre-
determined number of reflections or refractions.

In this work, we will only focus on the simulation of specular propagation
paths, and the parameter vector $\theta$ will encompass the permittivity,
conductivity and permeability of the materials in the scene. We leave the
study of diffusely scattered components, along with the calibration of
scattering coefficients and cross-polarization discrimination, for future
work. Accordingly, the RT output is assumed to be deterministic and given by
the function

$R(c|\theta)=\left\\{\alpha^{c}_{p}(\theta),\tau_{p}^{c},\Phi^{c,\mathrm{AoD}}_{p},\Phi^{c,\mathrm{AoA}}_{p}\right\\}_{p=1}^{P^{c}},$
(1)

which maps the input coordinates $c$ and material parameters $\theta$ into the
parameters of $P^{c}$ paths. By (1), each propagation path
$p\in\\{1,...,P^{c}\\}$ is described by a complex amplitude
$\alpha^{c}_{p}(\theta)\in\mathbbm{C}$, a delay $\tau_{p}^{c}>0$ [s], a pair
of angles of departure $\Phi^{c,\mathrm{AoD}}_{p}\in[0,\pi)\times[0,2\pi)$,
and a pair of angles of arrival
$\Phi^{c,\mathrm{AoA}}_{p}\in[0,\pi)\times[0,2\pi)$. The angles of departure
$\Phi^{c,\mathrm{AoD}}_{p}=(\Phi^{c,\mathrm{AoD}}_{\mathrm{el},p},\Phi^{c,\mathrm{AoD}}_{\mathrm{az},p})\in[0,\pi)\times[0,2\pi)$
include the elevation $\Phi^{c,\mathrm{AoD}}_{\mathrm{el},p}$ and azimuth
$\Phi^{c,\mathrm{AoD}}_{\mathrm{az},p}$ angles from the transmitter
perspective, and the angles of arrival
$\Phi^{c,\mathrm{AoA}}_{p}=(\Phi^{c,\mathrm{AoA}}_{\mathrm{el},p},\Phi^{c,\mathrm{AoA}}_{\mathrm{az},p})\in[0,\pi)\times[0,2\pi)$
encompass the elevation-azimuth angles at the receiver, as illustrated in Fig.
3. In order to keep the notation concise, in what follows, we omit the
position subscript $c$ from the path parameters $R(c|\theta)$ in (1), and from
the number of paths $P^{c}$.

Figure 3: Illustration of the definitions of angle path parameters for planar
antenna arrays. The red crosses represent the elements of each array. The
angles and directions of the departing and arriving signals, as well as the
positions of the array elements, are represented in the respective local
coordinates.

For a given receiver and transmitter position
$c=(c^{\mathrm{Rx}},c^{\mathrm{Tx}})$, signal propagation between each receive
and transmit antenna pair is assumed to follow the same path parameters
$R(c|\theta)$, while differing only in terms of phase. Accordingly, for each
simulated path $p\in\\{1,...,P\\}$, phase offsets between neighbouring antenna
elements are modeled as a $N^{\mathrm{Rx}}\times 1$ complex-valued receive
steering vector $a^{\mathrm{Rx}}(\Phi^{\mathrm{AoA}}_{p})$ and a
$N^{\mathrm{Tx}}\times 1$ transmit steering vector
$a^{\mathrm{Tx}}(\Phi^{\mathrm{AoD}}_{p})$, where the functions
$a^{\mathrm{Rx}}(\cdot)$ and $a^{\mathrm{Tx}}(\cdot)$ depend on the shape and
orientation of the receive and transmit array, respectively, and on the
transmission distance [40, Ch. 6].

For instance, in the far-field regime, the steering vectors
$a^{\mathrm{Rx}}(\Phi^{\mathrm{AoA}}_{p})$ can be approximated using a planar
wavefront model as

$a^{\mathrm{Rx}}(\Phi^{\mathrm{AoA}}_{p})=\begin{pmatrix}e^{j\frac{2\pi}{v}f^{c}(d^{\mathrm{Rx}}_{1})^{\top}u(\Phi^{\mathrm{AoA}}_{p})},\ldots,e^{j\frac{2\pi}{v}f^{c}(d^{\mathrm{Rx}}_{N^{\mathrm{Rx}}})^{\top}u(\Phi^{\mathrm{AoA}}_{p})}\end{pmatrix}^{\top},$
(2)

where $v=3\times 10^{8}$ m/s denotes the speed of light, while
$d^{\mathrm{Rx}}_{n}\in\mathbbm{R}^{3\times 1}$ is the position of the
$n\mbox{-th}$ receive array element and
$u(\Phi^{\mathrm{AoA}}_{p})\in\mathbbm{R}^{3\times 1}$ is a unitary-norm
vector pointing in the direction specified by $\Phi^{\mathrm{AoA}}_{p}$, both
expressed in the local coordinates of the receiver (see Fig. 3). The transmit
steering vector $a^{\mathrm{Tx}}(\Phi^{\mathrm{AoD}}_{p})$ can also be
computed in an analogous manner as

$a^{\mathrm{Tx}}(\Phi^{\mathrm{AoD}}_{p})=\begin{pmatrix}e^{j\frac{2\pi}{v}f^{c}(d^{\mathrm{Tx}}_{1})^{\top}u(\Phi^{\mathrm{AoD}}_{p})},\ldots,e^{j\frac{2\pi}{v}f^{c}(d^{\mathrm{Tx}}_{N^{\mathrm{Tx}}})^{\top}u(\Phi^{\mathrm{AoD}}_{p})}\end{pmatrix}^{\top},$
(3)

where $d^{\mathrm{Tx}}_{n}\in\mathbbm{R}^{3\times 1}$ is the position of the
$n\mbox{-th}$ transmit array element, and
$u(\Phi^{\mathrm{AoD}}_{p})\in\mathbbm{R}^{3\times 1}$ denotes the direction
specified by $\Phi^{\mathrm{AoD}}_{p}$, in the local transmitter coordinate
system (see Fig. 3). Furthermore, for curved wavefronts in near-field regimes,
the steering vectors can be approximated using a parabolic wavefront model
[41]. Alternatively, it is also possible to determine the vectors
$a^{\mathrm{Rx}}(\Phi^{\mathrm{AoA}}_{p})$ and
$a^{\mathrm{Tx}}(\Phi^{\mathrm{AoD}}_{p})$ numerically by running the RT
simulation for each antenna pair, i.e., at each coordinate
$c_{n^{\mathrm{Rx}},n^{\mathrm{Tx}}}=(c^{\mathrm{Rx}}+d^{\mathrm{Rx}}_{n^{\mathrm{Rx}}},c^{\mathrm{Tx}}+d^{\mathrm{Tx}}_{n^{\mathrm{Tx}}})$
for $n^{\mathrm{Rx}}\in\\{1,...,N^{\mathrm{Rx}}\\}$ and
$n^{\mathrm{Tx}}\in\\{1,...,N^{\mathrm{Tx}}\\}$ [13].

###  II-B Deterministic Channel Model

Using the multipath parameters $R(c|\theta)$ produced by the RT in (1), the DT
can reconstruct the frequency response of each path component
$p\in\\{1,...,P\\}$ for each subcarrier frequency $f_{s}\in B$ with
$s\in\\{1,...,S\\}$. This response is described by an
$N^{\mathrm{Rx}}N^{\mathrm{Tx}}\times 1$ vector of complex numbers, one for
each of the transmit-receive antenna pairs. It is given by

$\alpha_{p}(\theta)a^{s}(\tau_{p},\Phi^{\mathrm{AoA}}_{p},\Phi^{\mathrm{AoD}}_{p}),$
(4)

with

$a^{s}(\tau_{p},\Phi^{\mathrm{AoA}}_{p},\Phi^{\mathrm{AoD}}_{p})=e^{-j2\pi
f_{s}\tau_{p}}a^{\mathrm{Rx}}(\Phi^{\mathrm{AoA}}_{p})\otimes{a^{\mathrm{Tx}}(\Phi^{\mathrm{AoD}}_{p})}^{*},$
(5)

where ${\\{\cdot\\}}^{*}$ denotes the element-wise conjugate operation, and
$\otimes$ represents the Kronecker product [42]. Accordingly, each path is
characterized by a magnitude $|\alpha_{p}(\theta)|$ and by a phase
$(\angle\alpha_{p}(\theta)-2\pi f_{s}\tau_{p})$ dependent on the propagation
delay $\tau_{p}$ and on the phase $\angle\alpha_{p}(\theta)$ of the complex
amplitude $\alpha_{p}(\theta)$, in addition to the antenna-dependent
contributions of the steering vectors. The Kronecker product in
$a^{s}(\tau_{p},\Phi^{\mathrm{AoA}}_{p},\Phi^{\mathrm{AoD}}_{p})$ ensures that
the entry of the vector in (4) corresponding to any pair of transmit-receive
antennas is the product of the entries of the steering vectors for the
respective antennas.

Combining the contributions of all paths, the DT estimates the multi-path
channel at subcarrier frequency $f_{s}$ as the sum

$\begin{split}\hat{H}^{s}(c|\theta)&=\sum_{p=1}^{P}\alpha_{p}(\theta)a^{s}(\tau_{p},\Phi^{\mathrm{AoA}}_{p},\Phi^{\mathrm{AoD}}_{p})\\\
&=\sum_{p=1}^{P}\alpha_{p}(\theta)e^{-j2\pi
f_{s}\tau_{p}}a^{\mathrm{Rx}}(\Phi^{\mathrm{AoA}}_{p})\otimes{a^{\mathrm{Tx}}(\Phi^{\mathrm{AoD}}_{p})}^{*},\end{split}$
(6)

where the “hat” notation in (6) is a reminder that this is an estimate
produced by the DT on the basis of the RT output. Note that the more
conventional $N^{\mathrm{Rx}}\times N^{\mathrm{Tx}}$ matrix notation

$\sum_{p=1}^{P}\alpha_{p}(\theta)e^{-j2\pi
f_{s}\tau_{p}}a^{\mathrm{Rx}}(\Phi^{\mathrm{AoA}}_{p}){a^{\mathrm{Tx}}(\Phi^{\mathrm{AoD}}_{p})}^{\dagger},$
(7)

where ${\\{\cdot\\}}^{\dagger}$ represents the conjugate-transpose operation,
can be readily obtained by arranging the successive
$N^{\mathrm{Tx}}$-dimensional segments of vector $\hat{H}^{s}(c|\theta)$ into
the rows of a channel matrix.

Expanding from the per-subcarrier formulation in (4), we collect the
contributions of the $p\mbox{-th}$ path across all $S$ subcarriers into a
single $SN^{\mathrm{Rx}}N^{\mathrm{Tx}}\times 1$ vector
$\alpha_{p}(\theta)a(\tau_{p},\Phi^{\mathrm{AoA}}_{p},\Phi^{\mathrm{AoD}}_{p})$,
where

$a(\tau_{p},\Phi^{\mathrm{AoA}}_{p},\Phi^{\mathrm{AoD}}_{p})=w(\tau_{p})\otimes
a^{\mathrm{Rx}}(\Phi^{\mathrm{AoA}}_{p})\otimes{a^{\mathrm{Tx}}(\Phi^{\mathrm{AoD}}_{p})}^{*},$
(8)

represents the phase contributions emanating from all possible combinations of
antenna pairs, represented by the $N^{\mathrm{Rx}}N^{\mathrm{Tx}}\times 1$
vector
$a^{\mathrm{Rx}}(\Phi^{\mathrm{AoA}}_{p})\otimes{a^{\mathrm{Tx}}(\Phi^{\mathrm{AoD}}_{p})}^{*}$,
and subcarriers, represented by the $S\times 1$ vector

$w(\tau_{p})=\begin{pmatrix}e^{-j2\pi f_{1}\tau_{p}},...,e^{-j2\pi
f_{S}\tau_{p}}\end{pmatrix}^{\top}.$ (9)

Accordingly, we regroup the frequency responses
$\\{\hat{H}^{s}(c|\theta)\\}_{s=1}^{S}$ of all $S$ subcarriers into the
$SN^{\mathrm{Rx}}N^{\mathrm{Tx}}\times 1$ vector

$\begin{split}\hat{H}(c|\theta)&=(\hat{H}^{1}(c|\theta)^{\top},...,\hat{H}^{S}(c|\theta)^{\top})^{\top}\\\
&=\sum_{p=1}^{P}\alpha_{p}(\theta)a(\tau_{p},\Phi^{\mathrm{AoA}}_{p},\Phi^{\mathrm{AoD}}_{p}),\end{split}$
(10)

which constitutes the deterministic frequency response model assumed by the
DT.

###  II-C Phase Error Channel Model

The estimate $\hat{H}(c|\theta)$ in (10) hinges on the assumption that the
geometric properties provided by the DT are almost identical to the physical
environment of interest. In practice, as mentioned in Sec. I, even small
discrepancies between the geometry assumed by the RT and the actual physical
system, of a fraction of the studied wavelength, can significantly affect the
accuracy of the predicted phase $(\angle\alpha_{p}(\theta)-2\pi
f_{s}\tau_{p})$ of each path $p\in\\{1,...,P\\}$ [19, Ch. 7.5], resulting in
random small-scale fading effects that are not captured by the deterministic
model in (10). Therefore, even for high-fidelity DT models, one cannot depend
on the RT for a reliable estimate of the phases associated with each
propagation path [19, Ch. 7.5], [29, 31, 13].

In order to account for the inherent phase error of the DT model in (10), we
introduce a phase error term $z_{p}$ for each path $p\in\\{1,...,P\\}$ that is
distributed according to a von Mises distribution $\mathcal{VM}(0,\kappa_{0})$
with mean $0$ and concentration parameter $\kappa_{0}\geq 0$. Von Mises random
variables take values in the set $[-\pi,\pi)$, providing a natural model for
random phases [43, 44]. The probability density of a von Mises random variable
$X\sim\mathcal{VM}(\mu,\kappa)$ is defined as [45]

$\mathcal{VM}\left(x|\mu,\kappa\right)=\frac{1}{2\pi
I_{0}\left(\kappa\right)}e^{\kappa\cos\left(x-\mu\right)},$ (11)

for $x\in[-\pi,\pi)$, with parameters $\mu\in[-\pi,\pi)$ and $\kappa\geq 0$,
and where $I_{0}(\cdot)$ denotes the modified Bessel function of the first
kind of order $0$. The concentration parameter $\kappa$ determines the spread
of the distribution as illustrated in Fig. 4.

Figure 4: Von Mises probability density functions $\mathcal{VM}(0,\kappa)$ for
different values of the concentration parameters $\kappa$.

Incorporating the phase errors $\\{z_{p}\\}_{p=1}^{P}$ into the deterministic
model (10) yields the stochastic model

$\begin{split}&\hat{H}(c|\theta,\kappa_{0})=\sum_{p=1}^{P}e^{jz_{p}}\alpha_{p}(\theta)a(\tau_{p},\Phi^{\mathrm{AoA}}_{p},\Phi^{\mathrm{AoD}}_{p}),\end{split}$
(12)

where the phase errors
$\\{z_{p}\\}_{p=1}^{P}\overset{\mathrm{i.i.d.}}{\sim}\mathcal{VM}(0,\kappa_{0})$
are assumed to be independent and identically distributed (i.i.d.). Note that
all the subcarriers and antenna pairs share the same _phase error vector_

$Z=(z_{1},...,z_{P})^{T}\overset{\mathrm{i.i.d.}}{\sim}\mathcal{VM}(0,\kappa_{0})$
(13)

of i.i.d. phase error terms and, by extension, the same $P\times 1$ _phasor
error vector_

$e^{jZ}=(e^{jz_{1}},\ldots,e^{jz_{P}})^{\top},$ (14)

which corresponds to the element-wise application of the exponential function
to the $P\times 1$ phase-error vector $Z$. Accordingly, the stochastic model
(12) can be written as

$\hat{H}(c|\theta,\kappa_{0})=\hat{G}(c,\theta)e^{jZ},$ (15)

with the $SN^{\mathrm{Rx}}N^{\mathrm{Tx}}\times P$ matrix

$\hat{G}(c,\theta)=\begin{pmatrix}a(\tau_{1},\Phi^{\mathrm{AoA}}_{1},\Phi^{\mathrm{AoD}}_{1}),...,a(\tau_{P},\Phi^{\mathrm{AoA}}_{P},\Phi^{\mathrm{AoD}}_{P})\end{pmatrix}\mathrm{Diag}\left(\alpha(\theta)\right),$
(16)

where $\alpha(\theta)=(\alpha_{1}(\theta),...,\alpha_{P}(\theta))^{\top}$ is a
$P\times 1$ vector representing the complex amplitude of each of the $P$
predicted paths, and $\mathrm{Diag}(\cdot)$ denotes the diagonal matrix of the
given argument.

In practice, it may be useful to assume that different paths have different
concentration parameters $\kappa_{0}$. In particular, for line of sight (LoS)
paths one can set the corresponding phase error $z_{p}$ to zero, which
corresponds to an infinite concentration value $\kappa_{0}$. For other paths,
it is common to assume a uniform phase error, which is obtained by setting
$\kappa_{0}=0$, i.e., $\mathcal{VM}(z_{p}|0,0)=1/2\pi$ for
$z_{p}\in[-\pi,\pi)$. In this paper, we consider for simplicity of
presentation the i.i.d. model (13), but it is straightforward to generalize
the proposed technique to situations in which the phase errors are independent
but not identically distributed.

By decreasing the value of the concentration parameter $\kappa_{0}$, the
stochastic model (15) smoothly transitions from the deterministic-phase model
(10) for $\kappa_{0}\to+\infty$, which neglects phase errors, to a worst-case
phase-error RT model for $\kappa_{0}=0$, for which the distribution of the
phase error is uniform. Current applications of RT-simulated channels mostly
assume the uniform phase error setting, i.e., $\kappa_{0}=0$, for non-LoS
paths. Nevertheless, we keep the presentation general in order to evaluate the
potential benefits of more accurate ray tracing models.

###  II-D Data-Driven Calibration of the Ray Tracer

In addition to the geometric properties, the precision of the RT simulation
also depends on the accuracy of the estimated electromagnetic properties of
the materials. These are typically more complex to acquire than their
geometric counterpart [46, 47, 36, 27], and values estimated in technical
reports for generic materials do not capture the variety of properties found
in practical scenarios [28]. As a result, previous works [29, 30, 31, 32, 33,
34, 13] have proposed to include the unknown electromagnetic parameters,
namely the relative permittivity and conductivity, of all the material types
in the scene in a vector $\theta$, which is optimized based on data obtained
from the ground-truth scene at multiple locations in a process known as
_calibration_.

To this end, we assume the availability of a dataset

$\mathcal{D}=\left\\{\left(c_{n},H_{n}\right)\right\\}_{n=1}^{N}$ (17)

of $N$ _channel frequency response_ (CFR) observations
$H_{n}\in\mathbbm{C}^{SN^{\mathrm{Rx}}N^{\mathrm{Tx}}\times 1}$ for
$n\in\\{1,...,N\\}$, with each observation $H_{n}$ taken at receiver and
transmitter coordinates $c_{n}=(c^{\mathrm{Rx}}_{n},c^{\mathrm{Tx}}_{n})$ for
all $S$ subcarriers. Specifically, each observation
$H_{n}=((H^{1}_{n})^{\top},...,(H^{S}_{n})^{\top})^{\top}$ groups the observed
CFRs $H^{s}_{n}\in\mathbbm{C}^{N^{\mathrm{Rx}}N^{\mathrm{Tx}}\times 1}$ across
all the subcarrier frequencies $f_{s}$ indexed by $s\in\\{1,...,S\\}$.

The goal of calibration is to enhance the accuracy of the RT output in (1) for
locations $c$ at which channel data are available, as well as possibly for new
locations $c$ not included in the dataset, thus extending the DT simulation
capacity to new receiver and/or transmitter positions. That is, the DT aims to
provide a reliable simulation of the propagation and channel conditions
between any existing or new pair of receiver and transmitter devices at
coordinates $c$ by leveraging the output $R(c|\theta)$ of the RT. Such a
simulator can then be used by the DT to monitor the physical deployment of
interest, or to design optimal control policies without the need to collect
further data by training AI models in virtual space [26, 16, 11, 18, 24, 25].

##  III Calibration Baselines

In this section, we review two existing classes of calibration methods. The
first approach, proposed in [13], optimizes the electromagnetic parameters
$\theta$ given the channel observations $\mathcal{D}$ in (17) by assuming the
phase error-oblivious model $\hat{H}(c|\theta)$ in (10). We refer to this
channel response-based method as _phase error-oblivious calibration_. The
second approach, representing the line of work [29, 30, 31, 32, 33, 34],
calibrates the parameters $\theta$ by comparing the observed _channel power
profiles_ estimated from the dataset $\mathcal{D}$ to the corresponding
channel power profiles obtained from the RT under the uniform phase error
model $\hat{H}(c|\theta,0)$ in (15), i.e., with concentration parameter
$\kappa_{0}=0$. We refer to this power profile-based process as _uniform phase
error calibration_.

###  III-A Phase Error-Oblivious Calibration

As proposed in [13], the material parameters $\theta$ can be optimized by
comparing each observed CFR $H_{n}$ in the dataset $\mathcal{D}$ to its
respective simulated CFR in $\hat{H}(c_{n}|\theta)$ under the deterministic
model (10) for all observations $n\in\\{1,...,N\\}$. In this subsection, we
elaborate on the resulting phase error-oblivious calibration strategy.

Taking the measurement noise into account, each $n\mbox{-th}$ channel
observation $H_{n}$ is modelled by using an additive noise term as

$H_{n}=\hat{H}(c_{n}|\theta)+W_{n},$ (18)

for all $n\in\\{1,...,N\\}$, where the entries of the stochastic vector
$W_{n}\in\mathbbm{C}^{SN^{\mathrm{Rx}}N^{\mathrm{Tx}}\times 1}$ are i.i.d. as
circular Gaussian $\mathcal{CN}(0,\sigma^{2})$ variables with known noise
power $\sigma^{2}$. Accordingly, the distribution $P(H_{n}|c_{n},\theta)$ of
each observation $H_{n}$ is modelled as a circular Gaussian distribution

$P(H_{n}|c_{n},\theta)=\mathcal{CN}(H_{n}|\hat{H}(c_{n}|\theta),\sigma^{2}I)$
(19)

with mean $\hat{H}(c_{n}|\theta)$, where $I$ denotes the
$SN^{\mathrm{Rx}}N^{\mathrm{Tx}}\times SN^{\mathrm{Rx}}N^{\mathrm{Tx}}$
identity matrix.

The calibration process studied in [13] addresses the maximum likelihood (ML)
problem of finding the parameter $\theta^{*}$ that minimizes the cross-entropy
loss

$\mathcal{L}^{\mathrm{PEOC}}(\theta|\mathcal{D})=-\sum_{n=1}^{N}\log\left(P(H_{n}|c_{n},\theta)\right)$
(20)

for the available data $\mathcal{D}$. Using (19), the minimization of the
cross-entropy loss in (20) can be expressed as the least-squares problem

$\theta^{*}=\operatorname*{arg\,min}_{\theta}\sum_{n=1}^{N}\left\lVert
H_{n}-\hat{H}(c_{n}|\theta)\right\rVert^{2},$ (21)

where $\lVert X\rVert^{2}={X}^{\dagger}X$ denotes the squared 2-norm for a
complex vector $X$. Problem (21) is tackled in [13] through gradient descent.
This solution requires the differentiability of the predicted channel response
$\hat{H}(c_{n}|\theta)$ with respect to the electromagnetic parameters
$\theta$, which is ensured by Sionna RT [39, 13] .

###  III-B Uniform Phase Error Calibration

In earlier works [29, 30, 31, 32, 33, 34], RT calibration was conducted by
comparing the power profiles of the measured and simulated channels over time,
i.e., the power delay profiles [29, 31], and/or over space, i.e., the power-
angle-delay profiles [30, 32, 33, 34], while assuming uniform phase errors for
the RT-predicted path components. In this subsection, we briefly elaborate on
the resulting uniform phase error calibration approach.

Given an observation $H_{n}$ of the channel frequency response from the
dataset $\mathcal{D}$, an estimation of its power-angle-delay profile can be
obtained by projecting the CFR $H_{n}$ onto a channel vector matched to a
given triple $(\tau,\Phi^{\mathrm{AoD}},\Phi^{\mathrm{AoA}})$ of delay
$\tau\geq 0$, angle of departure
$\Phi^{\mathrm{AoD}}\in[0,\pi)\times[0,2\pi)$, and angle of arrival
$\Phi^{\mathrm{AoA}}\in[0,\pi)\times[0,2\pi)$. Such vector can be defined by
considering a normalized single-path component
$a(\tau,\Phi^{\mathrm{AoD}},\Phi^{\mathrm{AoA}})$ from (8) with the given path
parameters $(\tau,\Phi^{\mathrm{AoD}},\Phi^{\mathrm{AoA}})$ in lieu of
$(\tau_{p},\Phi^{\mathrm{AoD}}_{p},\Phi^{\mathrm{AoA}}_{p})$. Accordingly, the
_power-angle-delay profile_ for measured CFR $H_{n}$ is given as

$\mathrm{P}_{n}(\tau,\Phi^{\mathrm{AoD}},\Phi^{\mathrm{AoA}})=\left|\frac{1}{\sqrt{L}}{a(\tau,\Phi^{\mathrm{AoD}},\Phi^{\mathrm{AoA}})}^{\dagger}H_{n}\right|^{2},$
(22)

with $L=SN^{\mathrm{Rx}}N^{\mathrm{Tx}}$.

To enable calibration of the RT, the measured power-angle-delay profile (22)
must be compared to the simulated power-angle-delay profile obtained from the
RT under uniform phase errors. The latter can be obtained from the stochastic
model in (15) with concentration $\kappa_{0}=0$, in a manner similar to (22),
as

$\hat{\mathrm{P}}(c,\tau,\Phi^{\mathrm{AoD}},\Phi^{\mathrm{AoA}}|\theta)=\mathbbm{E}_{Z\overset{\mathrm{i.i.d.}}{\sim}\mathcal{VM}(0,0)}\left[\left|\frac{1}{\sqrt{L}}{a(\tau,\Phi^{\mathrm{AoD}},\Phi^{\mathrm{AoA}})}^{\dagger}\hat{G}(c,\theta)e^{jZ}\right|^{2}\right],$
(23)

where the average is taken with respect to the uniformly distributed phases
$Z=(z_{1},...,z_{P})^{\top}$. The average in (23) can be evaluated explicitly
under the assumption of i.i.d. uniform phase errors as (see Appendix A)

$\hat{\mathrm{P}}(c,\tau,\Phi^{\mathrm{AoD}},\Phi^{\mathrm{AoA}}|\theta)=\left\lVert\frac{1}{\sqrt{L}}{a(\tau,\Phi^{\mathrm{AoD}},\Phi^{\mathrm{AoA}})}^{\dagger}\hat{G}(c,\theta)\right\rVert^{2}.$
(24)

Optimization of the material parameters $\theta$ is carried by finding the
parameter

$\theta^{*}=\operatorname*{arg\,min}_{\theta}\mathcal{L}^{\mathrm{UPEC}}(\theta|\mathcal{D})$
(25)

that minimizes the discrepancy between the measured and simulated power
profiles for a set of pre-defined $M$ angle-delay triples
$\\{(\tau_{m},\Phi^{\mathrm{AoD}}_{m},\Phi^{\mathrm{AoA}}_{m})\\}_{m=1}^{M}$
under the squared loss metric

$\mathcal{L}^{\mathrm{UPEC}}(\theta|\mathcal{D})=\sum_{n=1}^{N}\sum_{m=1}^{M}\left|\mathrm{P}_{n}(\tau_{m},\Phi^{\mathrm{AoD}}_{m},\Phi^{\mathrm{AoA}}_{m})-\hat{\mathrm{P}}(c_{n},\tau_{m},\Phi^{\mathrm{AoD}}_{m},\Phi^{\mathrm{AoA}}_{m}|\theta)\right|^{2}.$
(26)

Problem (25) can be tackled through gradient descent, where differentiability
of the loss $\mathcal{L}^{\mathrm{UPEC}}(\theta|\mathcal{D})$ with respect to
the material parameters $\theta$ is ensured by the differentiability of the RT
output [13]. In this regard, we note that existing art has mostly relied on
less efficient strategies based on simulated annealing [29, 30, 31] or
exhaustive parameter search [32, 34].

The definition of criterion (26) requires the choice of the $M$ reference
triples
$\\{(\tau_{m},\Phi^{\mathrm{AoD}}_{m},\Phi^{\mathrm{AoA}}_{m})\\}_{m=1}^{M}$.
These may be selected from the predicted path triples
$\\{(\tau_{p},\Phi^{\mathrm{AoD}}_{p},\Phi^{\mathrm{AoA}}_{p})\\}_{p=1}^{P}$,
or as uniformly spaced values in a suitable space. For instance, one may
consider evenly spaced delays $\tau_{m}$ between the minimal and maximal
predicted times of arrival, while selecting the set of angles of departure
$\Phi^{\mathrm{AoD}}_{m}$ and arrival $\Phi^{\mathrm{AoA}}_{m}$ as
(approximately) evenly spaced points on a sphere via a spherical Fibonacci
lattice [48]. In practice, one may even restrict the points of the spherical
lattice to the surface spanned by the estimated elevation and azimuth angles
of departure or arrival, or to the relevant angles for the given antenna
patterns.

These choices of angle-delay triples are in line with previous works, which
have proposed to evenly spread delay points across a given delay-window [31],
or to select angle-delay points close to the angles and times of arrival of
the predicted paths [29, 30, 32, 33, 34]. That said, we note that the
calibration approaches in [30, 32, 33, 34] first carry out a pre-processing
step of the measured CFRs, which assumes the time-angle resolution of the
system to be sufficient to distinguish the different path contributions in the
measured data. Accordingly, such schemes first pair the predicted rays to
estimated measured rays by matching the predicted time of arrival, angle of
departure and angle of arrival of each simulated path to its nearest peak in
the measured power-angle-delay profile. The corresponding loss function is
then computed by comparing the predicted and measured power of each pair.

##  IV Phase Error-Aware Calibration

In this section, we introduce a novel calibration scheme that accounts for the
inherent phase uncertainty associated with the RT output $R(c|\theta)$ in (1)
via the stochastic channel model in (15). The proposed _phase error-aware
calibration_ scheme aims at extracting coherent information about the multi-
path phases for each channel realization based on the available data (17) with
the goal of improving the quality of the calibrated RT parameters $\theta$.
Therefore, the approach contrasts with the phase error-oblivious strategy of
[13], reviewed in the previous section, which hinges on the deterministic
model in (10), as well as with existing uniform phase error calibration
schemes, also described in previous section, which account only for the
received power.

###  IV-A Calibration Problem

In a manner similar to (18), the channel observations are modelled by
including an additive-noise term to the phase error-aware model in (15) as

$H_{n}=\hat{G}(c_{n},\theta)e^{jZ_{n}}+W_{n},$ (27)

where $W_{n}\sim\mathcal{CN}(0,\sigma^{2}I)$ is a
$SN^{\mathrm{Rx}}N^{\mathrm{Tx}}\times 1$ random vector of i.i.d. circular
Gaussian elements with zero mean and known power $\sigma^{2}$, and
$e^{jZ_{n}}=(e^{jz_{n,1}},...,e^{jz_{n,P}})^{\top}$ is the $P\times 1$ random
phasor vector accounting for phase errors. The phase errors
$Z_{n}=(z_{n,1},...,z_{n,P})^{\top}\overset{\mathrm{i.i.d.}}{\sim}\mathcal{VM}(0,\kappa_{0})$
are described by a prior concentration parameter $\kappa_{0}$. We first assume
parameter $\kappa_{0}$ to be fixed, and then we extend the proposed algorithm
to enable the estimation of the prior parameter $\kappa_{0}$ from data.

By (27), the probability distribution of each observation conditioned on the
phase-error vector $Z_{n}$ is given by

$P(H_{n}|Z_{n},c_{n},\theta)=\mathcal{CN}\left(H_{n}\Big{|}\hat{G}(c_{n},\theta)e^{jZ_{n}},\sigma^{2}I\right).$
(28)

The probability distribution $P(H_{n}|c_{n},\theta)$ of each measured CFR
$H_{n}$ is then obtained by marginalizing the joint distribution
$P(H_{n},Z_{n}|c_{n},\theta,\kappa_{0})=P(H_{n}|Z_{n},c_{n},\theta)P(Z_{n}|\kappa_{0})$
over the distribution $P(Z_{n}|\kappa_{0})$ of the phase-error vector $Z_{n}$,
i.e.,

$P(H_{n}|c_{n},\theta,\kappa_{0})=\int
P(H_{n}|Z_{n},c_{n},\theta)P(Z_{n}|\kappa_{0})dZ_{n},$ (29)

where integration is taken over the domain $[-\pi,\pi)^{\otimes P}$ of the
phase errors. Given the i.i.d. assumption, the phase-error vector is
distributed as
$P(Z_{n}|\kappa_{0})=\prod_{p=1}^{P}\mathcal{VM}(z_{n,p}|0,\kappa_{0})$.

The cross-entropy loss with respect to the RT parameters $\theta$ is then
given as

$\mathcal{L}^{\mathrm{PEAC}}(\theta|\mathcal{D})=-\sum_{n=1}^{N}\log\left(P(H_{n}|c_{n},\theta)\right).$
(30)

The resulting ML calibration problem is therefore defined as the minimization

$\theta^{*}=\operatorname*{arg\,min}_{\theta}\left\\{-\sum_{n=1}^{N}\log\left(\int
P(H_{n}|Z_{n},c_{n},\theta)P(Z_{n}|\kappa_{0})dZ_{n}\right)\right\\}.$ (31)

Problem (31) reduces to (21) when choosing an infinitely large prior
concentration parameter $\kappa_{0}$ for the phase noise.

###  IV-B Variational Expectation Maximization Algorithm

The direct solution of the problem (31) requires the expensive computation of
the negative log-likelihood function (30) by marginalizing over the latent
variables $Z_{n}$ (29) for all observations $n\in\\{1,...,N\\}$. To address
this computational challenge, we propose a computationally efficient
_variational expectation maximization_ (VEM) algorithm [38, Ch. 10] that
approximates the solution of the ML problem (31) by restricting the set of
possible solutions to a parametric family of distributions.

The proposed VEM approach iterates between the problem of optimizing the
parameter $\theta$ and that of estimating the phase errors $z_{n,p}$ for all
paths $p\in\\{1,...,P\\}$ and all observations $n\in\\{1,...,N\\}$. The main
underlying idea is that, while a priori all phase errors are assumed to have
the same distribution $\mathcal{VM}(0,\kappa_{0})$, each channel observation
$H_{n}$ is characterized by specific per-path phase errors, which may be
estimated and used to improve the calibration process.

The estimate of each phase error $z_{n,p}$ is represented by a von Mises
distribution $\mathcal{VM}\left(z_{n,p}|\mu_{n,p},\kappa_{n,p}\right)$ with
mean parameters $\mu_{n,p}$ and concentration parameter $\kappa_{n,p}$. The
mean represents the nominal estimate of the phase error, and the concentration
parameter quantifies the uncertainty of this estimate. As detailed below, both
parameters are estimated from the measured CFR $H_{n}$, as well as from the
current iterate $\theta$. Accordingly, the estimates of the phase errors for
all paths $p\in\\{1,...,P\\}$ associated with data point $n$ are described by
the _variational distribution_

$Q\left(Z_{n}|\mu_{n},\kappa_{n}\right)=\prod_{p=1}^{P}\mathcal{VM}\left(z_{n,p}|\mu_{n,p},\kappa_{n,p}\right),$
(32)

which assumes the estimated phase errors to be independent (but not
identically distributed). The adjective “variational” indicates that
distribution $Q\left(Z_{n}|\mu_{n},\kappa_{n}\right)$ depends on free
parameters that can be optimized, namely the mean parameter vector
$\mu_{n}=(\mu_{n,1},...,\mu_{n,P})^{\top}\in[-\pi,\pi)^{\otimes P}$ and the
concentration parameter vector
$\kappa_{n}=(\kappa_{n,1},...,\kappa_{n,P})^{\top}$ with $\kappa_{n,p}\geq 0$
for $p\in\\{1,...,P\\}$. As discussed, distribution (32) captures the coherent
information about the $p\mbox{-th}$ path that can be extracted from the data.
Considering all data points $n\in\\{1,...,N\\}$, this information is encoded
in the variational parameters $\\{\mu_{n},\kappa_{n}\\}_{n=1}^{N}$.

VEM approximates the solution of problem (31) by iteratively minimizing an
upper bound on the negative log-likelihood (31) over the material parameters
$\theta$ and over parameters $\\{\mu_{n},\kappa_{n}\\}_{n=1}^{N}$ of the
variational distribution (32). This upper bound, known as the _variational
free energy_ is given by (see, e.g., [49, 38])

$\mathcal{L}^{\mathrm{PEAC}}(\theta|\mathcal{D})\leq\sum_{n=1}^{N}\mathrm{F}_{n}(\mu_{n},\kappa_{n},\theta,\kappa_{0}),$
(33)

where

$\mathrm{F}_{n}\left(\mu_{n},\kappa_{n},\theta,\kappa_{0}\right)=\mathbbm{E}_{Z_{n}\sim
Q(Z_{n}|\mu_{n},\kappa_{n})}\left[\log\left(\frac{Q(Z_{n}|\mu_{n},\kappa_{n})}{P(H_{n},Z_{n}|c_{n},\theta,\kappa_{0})}\right)\right],$
(34)

and $P(H_{n},Z_{n}|c_{n},\theta,\kappa_{0})$ is the joint distribution defined
in (29). As detailed in Appendix B, using (28) and (32), the free energy (34)
can be written more explicitly as

$\begin{split}\mathrm{F}_{n}\left(\mu_{n},\kappa_{n},\theta,\kappa_{0}\right)=&\sum_{p=1}^{P}\left\\{\log\left(\frac{I_{0}(\kappa_{0})}{I_{0}(\kappa_{n,p})}\right)+b(\kappa_{n,p})\left(\kappa_{n,p}-\kappa_{0}\cos(\mu_{n,p})\right)\right\\}+L\log\left(\pi\sigma^{2}\right)\\\
&+\frac{1}{\sigma^{2}}\left[\left\lVert\hat{G}(c_{n},\theta)\mathrm{Diag}\left(B(\kappa_{n})\right)e^{j\mu_{n}}-H_{n}\right\rVert^{2}+L\left(\left|\alpha(\theta)\right|^{2}\right)^{\top}\left(\mathbbm{1}-B^{2}(\kappa_{n})\right)\right],\end{split}$
(35)

where $|\cdot|^{2}$ denotes the element-wise squared amplitude of a vector;
$\mathbbm{1}=(1,...,1)^{\top}$ denotes the all-ones $P\times 1$ vector;
$I_{\nu}(\cdot)$ is the modified Bessel function of the first kind of order
$\nu\geq 0$; $B(\kappa_{n})=(b(\kappa_{n,1}),...,b(\kappa_{n,P}))^{\top}$ is
the $P\times 1$ vector of Bessel ratios

$b(\kappa)=\frac{I_{1}(\kappa)}{I_{0}(\kappa)},$ (36)

for $\kappa\geq 0$; and
$B^{2}(\kappa_{n})=(b(\kappa_{n,1})^{2},...,b(\kappa_{n,P})^{2})^{\top}$ is
the $P\times 1$ vector of squared Bessel ratios. The Bessel ratio $b(\cdot)$
in (36) is a strictly increasing and continuous function in $[0,+\infty)$,
with $b(0)=0$ and $\lim_{\kappa\to+\infty}b(\kappa)=1$.

As in the classical expectation maximization (EM) algorithm, each iteration
$i\geq 0$ of VEM comprises of two steps, an _expectation step_ (E-step) and a
_maximization step_ (M-step), producing a sequence of iterates $\theta^{(i)}$
and $\\{\mu_{n}^{(i)},\kappa_{n}^{(i)}\\}_{n=1}^{N}$, with initializations
$\theta^{(0)}$ and $\\{\mu_{n}^{(0)},\kappa_{n}^{(0)}\\}_{n=1}^{N}$. In this
regard, the electromagnetic parameters $\theta^{(0)}$ may be initialized using
prior knowledge from previous observations or using known properties of the
materials in the scene [28]. The E-step estimates the statistics of the phase
errors $\\{\mu_{n}^{(i)},\kappa_{n}^{(i)}\\}_{n=1}^{N}$ across all data points
$(c_{n},H_{n})\in\mathcal{D}$ by using the current iterate $\theta^{(i)}$;
while the M-step updates the material parameters $\theta^{(i+1)}$ for fixed
updated phase error parameters
$\\{\mu_{n}^{(i+1)},\kappa_{n}^{(i+1)}\\}_{n=1}^{N}$. The overall proposed
phase error-aware calibration algorithm is summarized in Algorithm 1, and
details are provided next.

Algorithm 1 Phase Error-Aware Calibration

1:Training data $\mathcal{D}=\\{(c_{n},H_{n})\\}_{n=1}^{N}$; initial iterates
$(\theta^{(0)},\kappa_{0}^{(0)},\\{\mu_{n}^{(0)},\kappa_{n}^{(0)}\\}_{n=1}^{N})$
for the material parameters and for the per-data point estimated statistics of
the phase errors; number of iterations $I\geq 0$

2:Calibrated parameters $\theta^{(I)}$

3:for $i=0$ to $I-1$ do

4: Expectation Step :

5: for $n=1$ to $N$ do

6:
$\mu^{(i+1)}_{n}\leftarrow\angle\left\\{\left({\hat{G}(c_{n},\theta^{(i)})}^{\dagger}\hat{G}(c_{n},\theta^{(i)})\right)^{-1}\left(\frac{\sigma^{2}\kappa_{0}^{(i)}}{2}\mathbbm{1}+{\hat{G}(c_{n},\theta^{(i)})}^{\dagger}H_{n}\right)\right\\}$

7: for $p=1$ to $P$ do

8: if $\frac{L|\alpha_{p}(\theta^{(i)})|^{2}}{\sigma^{2}}>1$ then

9: $\kappa^{(i+1)}_{n,p}\leftarrow
2\sqrt{\frac{L|\alpha_{p}(\theta^{(i)})|^{2}}{\sigma^{2}}-1}\sqrt{\frac{L|\alpha_{p}(\theta^{(i)})|^{2}}{\sigma^{2}}}$

10: else

11: $\kappa^{(i+1)}_{n,p}\leftarrow 0$

12: end if

13: end for

14: end for

15: Maximization Step :

16:
$\theta^{(i+1)}\leftarrow\mathrm{GradientDescent}_{\theta}\left\\{\sum_{n=1}^{N}\mathrm{F}_{n}\left(\mu^{(i+1)}_{n},\kappa^{(i+1)}_{n},\theta,\kappa_{0}^{(i)}\right)\right\\}$

17: if
$(NP)^{-1}\sum_{n=1}^{N}\sum_{p=1}^{P}b(\kappa^{(i+1)}_{n,p})\cos(\mu^{(i+1)}_{n,p})\geq
0$ then

18: $\kappa_{0}^{(i+1)}\leftarrow
b^{-1}\left(\frac{1}{NP}\sum_{n=1}^{N}\sum_{p=1}^{P}b(\kappa^{(i+1)}_{n,p})\cos(\mu^{(i+1)}_{n,p})\right)$

19: else

20: $\kappa_{0}^{(i+1)}\leftarrow 0$

21: end if

22:end for

###  IV-C Expectation Step

As summarized by lines 5-13 of Algorithm 1, during the E-step of the proposed
phase error-aware calibration scheme, the material parameters $\theta^{(i)}$
are fixed, and minimization of the variational free energy in (35) is carried
out with respect to the local parameters $\mu_{n}$ and $\kappa_{n}$. This is
done by tackling the minimization problem

$\left(\mu_{n}^{(i+1)},\kappa_{n}^{(i+1)}\right)=\operatorname*{arg\,min}_{\mu_{n},\kappa_{n}}\left\\{\mathrm{F}_{n}\left(\mu_{n},\kappa_{n},\theta^{(i)},\kappa_{0}\right)\right\\}.$
(37)

This problem can be addressed, separately for each data point
$n\in\\{1,...,N\\}$, using any optimization algorithms such as gradient
descent. In the following, we present an approximate analytical solution that
will be tested experimentally in the next section.

As detailed in Appendix C, we propose to address a relaxed version of problem
(37) in which optimization is carried out over the complex phasor
$e^{j\mu_{n}}$ while removing the constraint that the magnitude of all entries
be equal to 1. For this relaxed problem, we obtain a stationary solution, and
then we project the solution into the space of unitary magnitude vectors to
recover an approximate estimate of the variational mean vector $\mu_{n}$. Note
that we do not make any claim of optimality for this solution.

To elaborate on the obtained approximate solution, let us assume that the
$L\times P$ matrix $\hat{G}(c_{n},\theta^{(i)})$ in (16) is of rank $P$. This
assumption requires that the number of channel coefficients measured for each
data points, $L$, is no smaller than the number, $P$, of paths, i.e., $L\geq
P$. In this common situation, full rank implies that the space-frequency
signatures of each of the $P$ paths, i.e., the columns of matrix
$\hat{G}(c_{n},\theta^{(i)})$, are linearly independent and hence
distinguishable using linear projections. Accordingly, this condition implies
that the system has sufficient resolution in the combined space-frequency
domain to resolve the different paths from the predicted channel response in
the given deployment scenario.

Under this assumption, by addressing the mentioned relaxed problem and
projecting back the solution into the space of phasors, we obtain the next
iterate

$\mu_{n}^{(i+1)}=\angle\left\\{\left({\hat{G}(c_{n},\theta^{(i)})}^{\dagger}\hat{G}(c_{n},\theta^{(i)})\right)^{-1}\left(\frac{\sigma^{2}\kappa_{0}}{2}\mathbbm{1}+{\hat{G}(c_{n},\theta^{(i)})}^{\dagger}H_{n}\right)\right\\},$
(38)

where $\angle(\cdot)$ represents the element-wise phase of a complex vector.
Furthermore, the concentration parameter vector $\kappa_{n,p}$ equals zero,
i.e., $\kappa^{(i+1)}_{n,p}=0$, if the inequality
$L|\alpha_{p}|^{2}/\sigma^{2}\leq 1$ holds, while it satisfies the fixed-point
equation

$\kappa_{n,p}^{(i+1)}=\frac{2L}{\sigma^{2}}|\alpha_{p}|^{2}B(\kappa_{n,p}^{(i+1)}).$
(39)

otherwise. While equation (39) can be addressed by fixed-point iterations, an
approximation can be obtained by bounding the Bessel ratio [50], yielding the
approximate closed-form solution

$\kappa_{n,p}^{(i+1)}=2\sqrt{\frac{L|\alpha_{p}|^{2}}{\sigma^{2}}-1}\sqrt{\frac{L|\alpha_{p}|^{2}}{\sigma^{2}}},$
(40)

for each path $p\in\\{1,...,P\\}$ and observation $n\in\\{1,...,N\\}$. Note
that, given the decomposition of the matrix $\hat{G}(c_{n},\theta^{(i)})$ in
(16), only the inverse of $\mathrm{Diag}(\alpha(\theta^{(i)}))$ needs to be
computed at each iteration. Accordingly, the E-step iteration can be computed
efficiently in $O(N\max(LP,P^{2}))$ operations, where the $O(\cdot)$ notation
refers to an asymptotic upper bound of the given expression [51, Ch. 3.1].

The condition

$\frac{L|\alpha_{p}|^{2}}{\sigma^{2}}>1,$ (41)

under which the fixed-point condition (39) applies, indicates that the signal-
to-noise ratio (SNR) associated with the observation of path $p$ is larger
than one. Accordingly, if the per-path SNR, i.e., the left-hand side of (41),
is smaller than one, the estimated concentration parameter
$\kappa_{n,p}^{(i+1)}$ corresponds to a uniform phase error, while, if it is
larger than one, the concentration parameter $\kappa_{n,p}^{(i+1)}$ is larger
than zero. For any strictly positive concentration $\kappa_{n,p}^{(i+1)}>0$,
the corresponding variational distribution
$\mathcal{VM}(\mu_{n,p}^{(i+1)},\kappa_{n,p}^{(i+1)})$ exhibits a preferred
direction along the mean phase $\mu_{n,p}^{(i+1)}$, and, therefore,
contributes to lowering the uncertainty of the phase error of path $p$ in
observation $n$.

###  IV-D Maximization Step

Let us consider first the case in which the concentration parameter
$\kappa_{0}$ is fixed. As summarized in the line 16 of Algorithm 1, during the
M-step of the proposed scheme, the previously selected phase error parameters
$\\{\mu_{n}^{(i+1)},\kappa_{n}^{(i+1)}\\}_{n=1}^{N}$, are given, and the
optimal material parameter iterate $\theta^{(i+1)}$ is obtained by tackling
the problem

$\theta^{(i+1)}=\operatorname*{arg\,min}_{\theta}\left\\{\sum_{n=1}^{N}\mathrm{F}_{n}\left(\mu_{n}^{(i+1)},\kappa_{n}^{(i+1)},\theta,\kappa_{0}\right)\right\\}.$
(42)

We suggest to use gradient descent for this purpose in a manner similar to the
phase error-oblivious scheme in [13]. In this case, the computational
complexity of the M-step is on par with the one of the considered baselines in
Sec. III.

When treating the prior concentration $\kappa_{0}$ as a parameter to be
jointly optimized with $\theta$ (lines 17-21 of Algorithm 1), the M-step
tackles the problem

$\left(\theta^{(i+1)},\kappa_{0}^{(i+1)}\right)=\operatorname*{arg\,min}_{\theta,\kappa_{0}}\left\\{\sum_{n=1}^{N}\mathrm{F}_{n}\left(\mu_{n}^{(i+1)},\kappa_{n}^{(i+1)},\theta,\kappa_{0}\right)\right\\}$
(43)

in lieu of (42). To address the minimization (43), parameter $\theta^{(i+1)}$
may still be updated through gradient descent; while, as derived in Appendix
D, an optimal value for the prior concentration iterate $\kappa_{0}^{(i+1)}$
is obtained as

$\kappa_{0}^{(i+1)}=b^{-1}\left(\frac{1}{NP}\sum_{n=1}^{N}\sum_{p=1}^{P}b(\kappa^{(i+1)}_{n,p})\cos(\mu^{(i+1)}_{n,p})\right),$
(44)

provided that the inequality
$(NP)^{-1}\sum_{n=1}^{N}\sum_{p=1}^{P}b(\kappa^{(i+1)}_{n,p})\cos(\mu^{(i+1)}_{n,p})\geq
0$ holds, and as $\kappa_{0}^{(i+1)}=0$ otherwise, where $b^{-1}(\cdot)$
denotes the inverse function of the Bessel ratio $b(\cdot)$ [52].

##  V Experiments

In this section, we present numerical experiments with the aim of offering
insights into the performance comparison between the existing calibration
schemes summarized in Sec. III-A and Sec. III-B and the proposed method
introduced in the previous section over a set of synthetic examples, allowing
us to control the evaluation setting. We adopt as performance metrics the
accuracy of calibration in terms of the errors on the estimated
electromagnetic material parameters $\theta$ and on the estimated received
signal power. Furthermore, we evaluate the RT capacity to predict received
powers for new coordinates of the receiver that were not covered by the
available data. The source code and 3D files of all the presented experiments
are available online [53].

###  V-A Setting

Using the Sionna RT module [39, 13], we simulate radio propagation for the
toy-example presented in Fig. 2, as well as for a urban scenario representing
the north side of King’s College London Strand campus, as depicted in Fig. 6.
Communication between the receive and transmit devices takes place at a
central carrier frequency $f^{c}=6$ GHz using a bandwidth $B=S\Delta f$ [Hz]
composed of $S$ evenly spaced subcarriers, with fixed spacing $\Delta f=30$
kHz. All antenna elements are simulated using an isotropic radiation pattern.
Antenna arrays are modelled using evenly spaced antenna elements, with a half-
wavelength spacing, and their steering vectors $a^{\mathrm{Rx}}(\cdot)$ and
$a^{\mathrm{Tx}}(\cdot)$ are computed using the planar wavefront assumption in
the far-field regime as presented in (2) and (3), respectively.

For each scenario, all surfaces are assumed to be made of the same non-
magnetic material, such that the properties of all the materials in the scene
can be entirely characterized by a single set of parameters
$\theta=\\{\epsilon,\gamma\\}$, where $\epsilon\geq 1$ denotes their relative
(real) _permittivity_ and $\gamma>0$ denotes their _conductivity_ [S/m].
Accordingly, channel observations $(c_{n},H_{n})\in\mathcal{D}$ are
synthetically generated by running the RT at the receive and transmit
locations $c_{n}=(c^{\mathrm{Rx}}_{n},c^{\mathrm{Tx}}_{n})$ under ground-truth
material parameters
$\theta^{\mathrm{true}}=\\{\epsilon^{\mathrm{true}},\gamma^{\mathrm{true}}\\}$,
which, following the ITU recommendation for concrete-like materials [20], are
set to $\epsilon^{\mathrm{true}}=5.31$ and $\gamma^{\mathrm{true}}=0.139$ S/m.
Channel responses $H_{n}$ are simulated using either the deterministic (Sec.
II-B) or the phase error (Sec. II-C) channel model, plus an additive Gaussian
noise term with known power $\sigma^{2}>0$. Accordingly, we define the SNR of
the observed data at coordinates $c$ as

$\mathrm{SNR}(c)=\frac{\mathrm{P}_{\mathrm{sig}}\left(c|\theta^{\mathrm{true}}\right)}{\sigma^{2}},$
(45)

where

$\mathrm{P}_{\mathrm{sig}}\left(c|\theta^{\mathrm{true}}\right)=\left\lVert\alpha(\theta^{\mathrm{true}})\right\rVert^{2}$
(46)

represents the total received power across all ground-truth propagation paths.

During calibration, the ground-truth material parameters
$\theta^{\mathrm{true}}$ are assumed to be unknown, and we leverage the given
dataset $\mathcal{D}$ to calibrate the material parameters estimate
$\theta=\\{\epsilon,\gamma\\}$. The geometric model available at the DT
presents small-scale differences with respect to the ground-truth scene used
to generate the data. These differences are generated either by modifying the
geometry of the ground-truth scene, as illustrated in Fig. 2(b), or by
injecting phase errors in the simulated paths responses (see Sec. V-C2).

Once the parameters $\theta=\\{\epsilon,\gamma\\}$ are calibrated using one of
the three schemes described in Sec. III and Sec. IV, we assess the relative
performance of calibration in terms of: (_i_) the normalized permittivity
error $|\epsilon-\epsilon^{\mathrm{true}}|/\epsilon^{\mathrm{true}}$; (_ii_)
the normalized conductivity error
$|\gamma-\gamma^{\mathrm{true}}|/\gamma^{\mathrm{true}}$; and (_iii_) the
normalized received power error

$\frac{\left|\hat{\mathrm{P}}_{\mathrm{sig}}\left(c|\theta\right)-\mathrm{P}_{\mathrm{sig}}\left(c|\theta^{\mathrm{true}}\right)\right|}{\mathrm{P}_{\mathrm{sig}}\left(c|\theta^{\mathrm{true}}\right)},$
(47)

where
$\hat{\mathrm{P}}_{\mathrm{sig}}(c|\theta)=\lVert\alpha(\theta)\rVert^{2}$ is
the total power received at position $c$ across all propagation paths
estimated using the DT model with calibrated parameters $\theta$. For all
calibration schemes, we initialize the parameters $\theta$ to be calibrated at
$\epsilon^{(0)}=3$ and $\gamma^{(0)}=0.1$ S/m. Throughout the experiments, the
phase error-aware scheme (Sec. IV) optimizes its prior concentration parameter
$\kappa_{0}$ as per the alternative M-step objective (43) with initialization
at $\kappa_{0}^{(0)}=0$. The uniform phase error scheme (Sec. III-B) uses a
set of $M=P$ projections in (26) defined by the delays and angles
$\\{(\tau_{p},\Phi^{\mathrm{AoD}}_{p},\Phi^{\mathrm{AoA}}_{p})\\}_{p=1}^{P}$
of the $P$ paths computed by the RT using the DT model.

###  V-B Toy Example

Following the discussion in Sec. I, we first compare the presented calibration
methods for the toy example presented in Fig. 2, with the receiver and
transmitter modeled as single isotropic antennas located on the plane at
positions $c^{\mathrm{Rx}}=(240,0)^{\top}$ and
$c^{\mathrm{Tx}}=(-240,0)^{\top}$, respectively, and the two horizontal walls
intersecting the vertical axis at positions $(0,100)^{\top}$ and
$(0,-180)^{\top}$, where the coordinates represent multiples of the carrier
wavelength $\lambda^{c}\approx 5$ cm. To model a two-dimensional scenario, we
discard the paths that reflect on the floor, and focus exclusively on the two
non-line of sight specular paths displayed in Fig. 2, which are caused by
reflections on the two horizontal walls. Accordingly, the two propagation
paths interfere constructively at the receiver coordinates, with a time of
arrival difference of $\Delta\tau\approx 13.3$ ns.

We consider a set of $N=50$ synthetically generated channel observations
$\\{H_{n}\\}_{n=1}^{N}\overset{\mathrm{i.i.d.}}{\sim}\mathcal{CN}(\hat{H}(c|\theta^{\mathrm{true}}),\sigma^{2})$
with $\sigma^{2}>0$, where $\hat{H}(c|\theta^{\mathrm{true}})$ is obtained
from the deterministic model in (10) by running the RT at the receive and
transmit locations $c=(c^{\mathrm{Rx}},c^{\mathrm{Tx}})$. During calibration,
we adopt the geometric model available at the DT which, as illustrated in Fig.
2(b), presents a small-scale discrepancy with respect to the ground-truth
geometry in the form of a shift of the position of the lower wall by an amount
of $\Delta d=0.4\lambda^{c}$ [m]. Accordingly, the lower wall intercepts the
vertical axis at coordinates $(0,-180.4)$ in the DT model, and the two
simulated paths interfere destructively at the receiver.

Figure 5:  Relative power estimation errors at the position $c$ used to
collect data for phase error-oblivious calibration (Sec. III-A), uniform phase
error calibration (Sec. III-B), and for the proposed phase error-aware
calibration (Sec. IV) as a function of the available bandwidth $B$, with a
signal-to-noise (SNR) ratio of $20$ dB. Lines represent the median error
across ten independent channel observation and calibration runs. Shaded areas
represent the first and third quartiles.

We plot the power estimation errors of all schemes as a function of the
bandwidth $B$ in Fig. 5 for an SNR equal to $20$ dB. To scale the bandwidth
$B$, the subcarrier spacing $\Delta f=30$ kHz is kept constant, while the
number of subcarriers $S$ is increased. Confirming the insights in Fig. 2, the
performance of the power profile-based scheme starts improving once the system
bandwidth approaches $B=1/\Delta\tau=75$ MHz, i.e., once its temporal
resolution is high enough to separate the contribution of each path based on
their arrival times at the receiver. Accordingly, the relative power
estimation error of the uniform phase error scheme decreases from $0$ dB at
bandwidths lower than $B<20$ MHz, to $-11$ dB at bandwidth $B=500$ MHz. As the
number of subcarriers increases for larger bandwidths, the phase error-
oblivious scheme is affected by a larger number of phase errors in the path
responses simulated using the DT model, yielding an error increase from $-6$
dB at lower bandwidths $B<5$ MHz, to $-1$ dB at $B=500$ MHz. On the contrary,
owing to the larger number of subcarriers at higher bandwidths, the phase
error-aware scheme has access to more information to estimate the phase
errors, outperforming the phase error-oblivious and uniform phase error
schemes for bandwidths of $B=2$ MHz and above. The relative power error of the
phase error-aware scheme stabilizes around $-27$ dB for bandwidths larger than
$B>50$ MHz, which correspond to time resolutions at which the system can
resolve the two propagation paths.

###  V-C Urban Setting

Throughout the rest of this section, we study a more realistic scenario that
simulates radio propagation on the north side of King’s College London Strand
campus, whose 3D representation, depicted in Fig. 6, is obtained using geo-
spatial data from the OpenStreetMap database [27]. We now consider $8\times 8$
two-dimensional antenna arrays at both receive and transmit locations, and we
fix the number of subcarriers to $S=64$. Non-line of sight propagation paths
are simulated through specular reflection only, with up to ten reflections per
path. We leave the evaluations under diffracted, refracted and diffusely
reflected paths for future work.

Figure 6:  3D model of King’s College London Strand campus built using data
from the OpenStreetMap database [27]. As an illustration of the operation of
the ray tracer (RT), we show RT-simulated wireless propagation paths (green
lines) for a transmitter (red dot) and receiver (blue dot) pair.

To enable calibration, we use synthetically generated channel observations for
a single pair of receiver and transmitter devices at coordinates
$c_{\mathrm{cal}}=(c^{\mathrm{Rx}}_{\mathrm{cal}},c^{\mathrm{Tx}}_{\mathrm{cal}})$
depicted in Fig. 6, with the receiver and transmitter at heights $1$ and $5$
meters, respectively. Specifically, we generate a synthetic dataset
$\mathcal{D}=\\{(c_{\mathrm{cal}},H_{n})\\}_{n=1}^{N}$ of $N=50$ channel
observations by running the RT in the same 3D model available at the DT with
ground-truth material parameters $\theta^{\mathrm{true}}$. We account for the
discrepancy between the ground-truth (simulated) deployment and the virtual
geometry assumed by the DT in one of the following ways.

####  V-C1 Receiver position mismatch

A first approach is to directly model small geometric discrepancies by
assuming that the position of the receiver is subject to uncertainty. Each
channel observation $H_{n}$ for $n\in\\{1,...,N\\}$ is obtained by placing the
receiver at the position

$\tilde{c}^{\mathrm{Rx}}_{n}=c^{\mathrm{Rx}}_{\mathrm{cal}}+\eta^{\mathrm{true}}\lambda^{c}u_{n},$
(48)

where $\lambda^{c}\approx 5$ cm is the carrier wavelength,
$u_{n}\in\mathbbm{R}^{3\times 1}$ is a uniformly sampled direction in the unit
sphere, i.e., $\lVert u_{n}\rVert=1$, and $\eta^{\mathrm{true}}\in[0,0.5]$
dictates the level of geometric discrepancy, with $\eta^{\mathrm{true}}=0$
indicating a perfect geometric model at the DT. Accordingly, the position
$c^{\mathrm{Rx}}_{\mathrm{cal}}$ assumed by the DT differs from the true
position $\tilde{c}^{\mathrm{Rx}}_{n}$ of observation $H_{n}$ due to the
random displacement $\eta^{\mathrm{true}}\lambda^{c}u_{n}$.

####  V-C2 Independent phase errors

The stochastic phase errors caused by the geometric displacement (48) are
generally correlated across different paths, since they stem from the same
uncertainty. In order to capture a worst case-scenario in terms of phase
errors caused by a geometric mismatch between the DT and the PT, we also study
a setting with independent phase errors. To this end, observations
$\\{H_{n}\\}_{n=1}^{N}$ are generated as independent samples from the phase-
error model
$\hat{H}(c_{\mathrm{cal}}|\theta^{\mathrm{true}},\kappa_{0}^{\mathrm{true}})$
in (27). The ground-truth concentration parameter $\kappa_{0}^{\mathrm{true}}$
controls the level of geometric discrepancy, with
$\kappa_{0}^{\mathrm{true}}=0$ corresponding to the most challenging case of
uniform phase errors, and $\kappa_{0}^{\mathrm{true}}\to+\infty$ corresponding
to a perfect DT model.

For both data generation methods, the SNR is defined as in (45), with the
caveat that the power $\mathrm{P}_{\mathrm{sig}}(c|\theta^{\mathrm{true}})$ at
the numerator is averaged over the positions given by the random displacements
in (48) for the former method.

###  V-D Calibration Performance

We start by comparing the three considered calibration schemes for geometric
discrepancies simulated through receiver position uncertainty (Fig. 7(a)), and
through independent phase errors (Fig. 7(b)), for an SNR of $20$ dB. Using the
three metrics described in Sec. V-A, Fig. 7 displays the performance of each
scheme in terms of normalized errors of the predicted permittivity $\epsilon$,
conductivity $\gamma$, and received power
$\hat{\mathrm{P}}_{\mathrm{sig}}(c_{\mathrm{cal}}|\theta)$ at coordinates
$c_{\mathrm{cal}}$. The standard deviation in Fig. 7(b) is a function of the
ground-truth concentration parameter $\kappa_{0}^{\mathrm{true}}$, and is
chosen here as it is easier to interpret as compared to the concentration
parameter. Note that, in the uniform phase error case, the standard deviation
equals $\pi/\sqrt{3}\mbox{ rad}\approx 104\mbox{ deg}$.

(a)

(b)

Figure 7:  Calibration relative errors in terms of estimation of the
electromagnetic properties (top) and of power estimation at the receiver
collecting data (bottom) for phase error-oblivious calibration (Sec. III-A),
uniform phase error calibration (Sec. III-B), and for the proposed phase
error-aware calibration (Sec. IV) as a function of the magnitude of the random
receiver displacements (a), and as a function of the standard deviation of the
independently sampled phase errors (b). Lines represent the median error
across ten independent channel observation and calibration runs, with a
signal-to-noise (SNR) ratio of $20$ dB. Shaded areas represent the first and
third quartiles.

A first general observation is that the three considered calibration schemes
display similar relative performance gaps for both discrepancy models. This
motivates the adoption of the worst case of independently and uniformly
distributed ($\kappa_{0}^{\mathrm{true}}=0$) phase errors at each path for the
following experiments. Furthermore, from the power error plots in the bottom
row of Fig. 7, we observe that the proposed phase error-aware scheme scheme
outperforms the uniform phase error baseline across all considered random
receiver displacements amplitudes and phase error standard deviation levels,
while also performing better than the phase error-oblivious scheme for
displacements amplitudes above $\eta^{\mathrm{true}}=0.1$ and phase error
standard deviations of $20$ degrees and higher. In particular, in the most
common regime of uniform phase errors, i.e., at a standard deviation of
approximately $104$ degrees, the proposed scheme outperforms the two baselines
in terms of power error, with a performance gap of about $20$ dB when compared
to the uniform phase error scheme, and about $32$ dB when compared to the
phase error-oblivious scheme. We also observe from the two top rows of Fig 7
that an accurate estimate of the true material parameters $\theta$ is not a
pre-requisite to achieve reliable power estimation, as multiple combinations
of permittivity $\epsilon$ and conductivity $\gamma$ values can result in the
same predicted power.

Another general observation is that the proposed phase error-aware scheme is
able to benefit from the increased diversity in the data set produced by
larger random displacements, in Fig. 7(a), or larger phase error standard
deviation, in Fig. 7(b). In fact, the figure illustrates how phase error-aware
calibration yields improved estimates of material parameters as the
displacement amplitude and the phase error standard deviations increase,
which, in turn, also improves the power estimation. In this regard, the
diversity in the calibration data set resulting from randomness in the
geometry or path phases has the effect of mitigating the inherent limited
resolution arising from a limited bandwidth by allowing the observation of
different randomized path combinations in the $N$ channel observations. For
example, with a single observation, i.e., $N=1$ (not shown), the power error
remains around $-31\pm 1$ dB across all phase error standard deviations, while
with $N=50$, as illustrated in Fig. 7(b), the power error drops from $-27$ dB
when there are no phase errors to $-36$ dB for uniform phase errors.

Figure 8:  Calibration relative errors in terms of estimation of the
electromagnetic properties (top) and of power estimation at the receiver
collecting data (bottom) for phase error-oblivious calibration (Sec. III-A),
uniform phase error calibration (Sec. III-B), and for the proposed phase
error-aware calibration (Sec. IV) as a function of the SNR for uniformly
sampled independent phase errors ($\kappa_{0}^{\mathrm{true}}=0$). Lines
represent the median error across ten independent channel observation and
calibration runs. Shaded areas represent the first and third quartiles.

Fig. 8 displays the calibration performance of all schemes under uniform
ground-truth phase errors as a function of the SNR. The performance gains of
the proposed phase error-aware scheme range from around $13$ dB at an SNR of
$0$ dB to around $25$ dB at an SNR of $30$ dB in terms of power error.
Moreover, the phase error-oblivious calibration scheme is unable to correctly
calibrate its material parameters under uniform phase errors, yielding
inaccurate material and power estimates at the calibration location
$c_{\mathrm{cal}}$ for all studied SNR levels, while the best benchmark in
terms of power estimation is given by the uniform phase error scheme.

###  V-E Generalization Performance for Received Power Prediction

We now study the implications of calibration on the operations of the DT for
the urban scenario in Sec. V-C. Specifically, we consider the task of
predicting received power at positions of the receiver that are different from
that used to collect the calibration data set. This allows the DT to build
power maps, which can be used for functionalities such as network optimization
and resource allocation.

To this end, we compare the average signal power predicted by the RT using the
calibrated parameters $\theta$ to the average ground-truth signal power
obtained by using the ground-truth parameters $\theta^{\mathrm{true}}$ for all
receiver positions within coverage of the transmitter. Calibration data
$\mathcal{D}$ is generated using independently sampled phase errors at each
path, as presented in Sec. V-C2, with concentration parameter
$\kappa_{0}^{\mathrm{true}}\geq 0$. We emphasize that the transmitter’s
position $c^{\mathrm{Tx}}_{\mathrm{cal}}$ is fixed across calibration and
testing, while the position of the receiver $c^{\mathrm{Rx}}_{\mathrm{cal}}$
displayed in Fig. 6 describes the conditions used to collect the calibration
data.

No Phase Error

Phase Error-Aware Calibration

(a)

Phase Error-Oblivious Calibration

(b)

Uniform Phase Error Calibration

(c)

Uniform Phase Error

(d)

(e)

(f)

Figure 9:  Relative estimation errors (in dB) for the average predicted power
of the signal sent by a single transmitter (red dot) under phase error-aware
calibration ((a) and (d)), phase error-oblivious calibration ((b) and (e)),
and uniform phase error calibration ((c) and (f)) in the absence of phase
errors $\kappa_{0}^{\mathrm{true}}\to+\infty$ ((a), (b) and (c)), and for
uniform phase errors $\kappa_{0}^{\mathrm{true}}=0$ ((d), (e) and (f)). The
estimation errors at each position are averaged across ten independent channel
observation and calibration procedures, with a signal-to-noise ratio (SNR) of
$20$ dB during calibration.

Given a new receiver position $c^{\mathrm{Rx}}_{\mathrm{new}}$, the ground-
truth received power
$\mathrm{P}_{\mathrm{sig}}(c_{\mathrm{new}}|\theta^{\mathrm{true}})$ and the
RT-predicted signal power
$\hat{\mathrm{P}}_{\mathrm{sig}}(c_{\mathrm{new}}|\theta)$ in (47) are
computed at receiver and transmitter position
$c_{\mathrm{new}}=(c^{\mathrm{Rx}}_{\mathrm{new}},c^{\mathrm{Tx}}_{\mathrm{cal}})$.
Accordingly, Fig. 9 displays the normalized power prediction error (47)
between predicted and true received powers in dB as a function of the position
of the receiver. For each row in Fig. 9, the relative power error is averaged
across a set of ten independent channel observation and calibration procedures
with SNR equal to $20$ dB. The top row in Fig. 9 shows the power prediction
error in the absence of phase errors during calibration, i.e.,
$\kappa_{0}^{\mathrm{true}}\to+\infty$, while the bottom row displays the
power errors under uniformly distributed phase errors, i.e.,
$\kappa_{0}^{\mathrm{true}}=0$.

In line with the observations made on the calibration performance of the phase
error-oblivious calibration strategy in Sec. V-D, the relative prediction
error under this scheme is satisfactory, below $-20$ dB, for zero phase error
(top row), while, for uniform phase errors (bottom row), the estimation error
rises above $-5$ dB in zones with LoS coverage and around $0$ dB in most
locations relying exclusively on non-LoS components, where propagation
conditions are more sensitive to calibration inaccuracies. In contrast, the
uniform phase error calibration baseline in Fig. 9(f) is able to keep the
power estimation error below $-15$ dB in locations with LoS components, and
around $-3$ dB in locations covered exclusively by non-LoS components.
Finally, the estimation errors of the proposed phase error-aware scheme in
Fig. 9(a) and 9(d) are below $-20$ dB across all covered receiver’s positions,
even in the worst case of uniform phase error.

The relative power predition error is further studied in Fig. 10, for which
the relative error is averaged over all covered receiver’s positions and
plotted as a function of the phase error standard deviation in Fig. 10(a), and
as a function of the SNR in Fig. 10(b). The figure confirms that the proposed
scheme outperforms both baselines in terms of its capacity to predict the
received power for new positions across a different phase error and SNR
levels, as long as the phase error is not negligible.

(a)

(b)

Figure 10:  Average estimation error of the predicted signal powers across all
new positions $c_{\mathrm{new}}$ within coverage of the transmitter under
ground-truth conditions as a function of the standard deviation of the phase
error (a), and as a function of the signal-to-noise ratio (SNR) (b); for the
phase error aware, phase error oblivious, and uniform phase error calibration
schemes. Lines represent the median error across ten independent channel
observation and calibration runs, with a $20$ dB SNR during calibration in
(a), and a uniform phase error during calibration in (b). Shaded areas
represent the first and third quartiles.

##  VI Conclusion

In this work, we have proposed a new phase error-aware method to calibrate the
material properties of RT simulations based on multi-carrier channel response
observations. The proposed scheme estimates RT-predicted phase-errors from
data, accounting for small-scale discrepancies between the assumed three-
dimensional DT model and the ground-truth geometry. The phase-error estimates
are used to correct the predicted interference pattern among propagation paths
during calibration of the material parameters. This is demonstrated to yield
more precise conductivity and permittivity estimates, as well as more accurate
signal power predictions, especially in the low-resolution regime
characterized by low-bandwidth conditions. We conducted quantitative
experiments for a single-material scene, showcasing the advantages of this
approach against current state-of-the-art power profile-based and channel
response-based schemes for practical conditions characterized by unreliable RT
phase predictions.

Future work may investigate the simulation of diffusely scattered propagation
paths, along with the calibration of scattering coefficients and cross-
polarization discrimination; the validation of the proposed approach on real-
world measurements; the impact of antenna patterns on phase uncertainty; and
the extension of the current study to a more diverse set of materials, such as
position-dependent material properties through neural network parameterization
[54].

##  Appendix A Average Signal Power under Von Mises Distributed Phases

Under the phase error model (27) with electromagnetic parameters $\theta$ and
phase-error concentration parameter $\kappa_{0}$, the average signal power at
a location $c_{n}\in\mathbbm{R}^{3}\times\mathbbm{R}^{3}$ per subcarrier and
antenna pair is computed as

$\begin{split}\hat{\mathrm{P}}_{\mathrm{sig}}\left(c_{n}|\theta,\kappa_{0}\right)=\frac{1}{L}\mathbbm{E}_{Z_{n}\overset{\mathrm{i.i.d.}}{\sim}\mathcal{VM}(0,\kappa_{0})}\left[\left\lVert\hat{G}(c_{n},\theta)e^{jZ_{n}}\right\rVert^{2}\right],\end{split}$
(49)

where $L=SN^{\mathrm{Rx}}N^{\mathrm{Tx}}$ and
$Z_{n}=(z_{n,1},...,z_{n,P})^{\top}$ with
$\\{z_{n,p}\\}_{p=1}^{P}\overset{\mathrm{i.i.d.}}{\sim}\mathcal{VM}(0,\kappa_{0})$.
Using the invariance of the trace to cyclic permutations[55], the expectation
term in (49) can be evaluated as

$\displaystyle\mathbbm{E}\left[\left\lVert\hat{G}(c_{n},\theta)e^{jZ_{n}}\right\rVert^{2}\right]$
$\displaystyle=\mathbbm{E}\left[\operatorname{Tr}\left\\{{\left(e^{jZ_{n}}\right)}^{\dagger}{\hat{G}(c_{n},\theta)}^{\dagger}\hat{G}(c_{n},\theta)e^{jZ_{n}}\right\\}\right]$
(50a)
$\displaystyle=\left\lVert\hat{G}(c_{n},\theta)\mathbbm{E}\left[e^{jZ_{n}}\right]\right\rVert^{2}+\operatorname{Tr}\left\\{\hat{G}(c_{n},\theta)\mathrm{Cov}\left[e^{jZ_{n}}\right]{\hat{G}(c_{n},\theta)}^{\dagger}\right\\},$
(50b)

where we have indicated $\mathbbm{E}[\cdot]$ for
$\mathbbm{E}_{Z_{n}\overset{\mathrm{i.i.d.}}{\sim}\mathcal{VM}(0,\kappa_{0})}[\cdot]$,
$\mathrm{Cov}(e^{jZ_{n}})$ denotes the $P\times P$ covariance matrix of the
random phasor $e^{jZ_{n}}$ for
$Z_{n}\overset{\mathrm{i.i.d.}}{\sim}\mathcal{VM}(0,\kappa_{0})$, and
$\operatorname{Tr}\\{\cdot\\}$ denotes the trace of a matrix. The mean and
variance terms in (50b) are then evaluated as [45]

$\begin{split}&\mathbbm{E}\left[e^{jZ_{n}}\right]=b\left(\kappa_{0}\right)\mathbbm{1}\\\
&\mathrm{Cov}\left[e^{jZ_{n}}\right]=\left(1-b(\kappa_{0})^{2}\right)I,\end{split}$
(51)

where $\mathbbm{1}=(1,...,1)^{\top}$ denotes the $P\times 1$ all-ones vector,
$I$ is the $P\times P$ identity matrix, and $b(\kappa_{0})$ is the is the
Bessel ratio in (36) for $\kappa\geq 0$.

Plugging (51) into (50), the average signal power in (49) can be evaluated as

$\hat{\mathrm{P}}_{\mathrm{sig}}\left(c_{n}|\theta,\kappa_{0}\right)=b(\kappa_{0})^{2}\frac{\lVert\hat{G}(c_{n},\theta)\mathbbm{1}\rVert^{2}}{L}+\left(1-b(\kappa_{0})^{2}\right)\left\lVert\alpha(\theta)\right\rVert^{2},$
(52)

where $\alpha(\theta)$ is the $P\times 1$ vector of complex paths coefficients
in (16). The final result in (52) writes the overall received power as the sum
of the average power predicted by the RT in the absence of phase errors, i.e.,
with $\kappa_{0}=+\infty$ and hence $b(\kappa_{0})=1$, and of the average
power under uniform phase errors for each propagation path, i.e., with
$\kappa_{0}=0$ and thus $b(\kappa_{0})=0$.

The result presented in (52) can be extended to any projection
$\hat{G}^{{}^{\prime}}(c_{n},\theta)={a}^{\dagger}\hat{G}(c_{n},\theta)$ of
the $L\times P$ matrix $\hat{G}(c_{n},\theta)$ by an $L\times 1$ vectors $a$
as

$\begin{split}\hat{\mathrm{P}}_{\mathrm{sig}}\left(c_{n}|\theta,\kappa_{0}\right)&=\frac{1}{L}\mathbbm{E}_{Z_{n}\overset{\mathrm{i.i.d.}}{\sim}\mathcal{VM}(0,\kappa_{0})}\left[\left\lVert\hat{G}^{{}^{\prime}}(c_{n},\theta)e^{jZ_{n}}\right\rVert^{2}\right]\\\
&=b(\kappa_{0})^{2}\frac{|{a}^{\dagger}\hat{G}(c_{n},\theta)\mathbbm{1}|^{2}}{L}+\left(1-b(\kappa_{0})^{2}\right)\left\lVert{a}^{\dagger}\hat{G}(c_{n},\theta)\right\rVert^{2}.\end{split}$
(53)

##  Appendix B Closed-Form Expression of the Free Energy

Following the notation in Sec. IV, this appendix summarizes the main steps to
derive a closed-form expression for the variational free energy in (35). Given
that phase errors are assumed to be independent, both under the prior
$P(Z_{n}|\kappa_{0})=\prod_{p=1}^{P}\mathcal{VM}(z_{n,p}|0,\kappa_{0})$ and
under the variational posterior distribution
$Q(Z_{n}|\mu_{n},\kappa_{n})=\prod_{p=1}^{P}\mathcal{VM}(z_{n,p}|\mu_{n,p},\kappa_{n,p})$,
the free energy decomposes as

$\begin{split}\mathrm{F}_{n}\left(\mu_{n},\kappa_{n},\theta,\kappa_{0}\right)=&\sum_{p=1}^{P}\left\\{\mathcal{H}\left(\mu_{n,p},\kappa_{n,p}\big{|}\big{|}0,\kappa_{0}\right)-\mathcal{H}\left(\mu_{n,p},\kappa_{n,p}\right)\right\\}\\\
&+\mathbbm{E}_{Z_{n}\sim
Q(Z_{n}|\mu_{n},\kappa_{n})}\left[-\log\left(P(H_{n}|Z_{n},c_{n},\theta)\right)\right],\end{split}$
(54)

where

$\begin{split}\mathcal{H}\left(\mu,\kappa\right)&=\mathbbm{E}_{z\sim\mathcal{VM}(\mu,\kappa)}\left[-\log\left(\mathcal{VM}\left(z|\mu,\kappa\right)\right)\right]\\\
&=\log\left(2\pi I_{0}(\kappa)\right)-\kappa b(\kappa)\end{split}$ (55)

denotes the differential entropy of the von Mises distribution
$\mathcal{VM}(\mu,\kappa)$ [56];

$\begin{split}\mathcal{H}\left(\mu,\kappa\big{|}\big{|}\mu_{0},\kappa_{0}\right)&=\mathbbm{E}_{z\sim\mathcal{VM}(\mu,\kappa)}\left[-\log\left(\mathcal{VM}\left(z|\mu_{0},\kappa_{0}\right)\right)\right]\\\
&=\log\left(2\pi
I_{0}(\kappa_{0})\right)-\kappa_{0}b(\kappa)\cos(\mu-\mu_{0})\end{split}$ (56)

is the cross-entropy between von Mises distributions
$\mathcal{VM}(\mu,\kappa)$ and $\mathcal{VM}(\mu_{0},\kappa_{0})$ [56]. Using
the distribution $P(H_{n}|Z_{n},c_{n},\theta)$ in (28), the remaining
expectation term in (54) can be expressed as

$\displaystyle\mathbbm{E}_{Z_{n}\sim
Q(Z_{n}|\mu_{n},\kappa_{n})}\left[-\log\left(P(H_{n}|Z_{n},c_{n},\theta)\right)\right]$
$\displaystyle=\log\left(\pi^{L}\sigma^{2L}\right)+\frac{1}{\sigma^{2}}\mathbbm{E}\left[\left\lVert
H_{n}-\hat{G}(c_{n},\theta)e^{jZ_{n}}\right\rVert^{2}\right]$ (57a)
$\displaystyle=L\log\left(\pi\sigma^{2}\right)+\frac{1}{\sigma^{2}}\left[\left\lVert\hat{G}(c_{n},\theta)\mathbbm{E}\left[e^{jZ_{n}}\right]-H_{n}\right\rVert^{2}+\operatorname{Tr}\left\\{\hat{G}(c_{n},\theta)\mathrm{Cov}(e^{jZ_{n}}){\hat{G}(c_{n},\theta)}^{\dagger}\right\\}\right]$
(57b)

with $\mathbbm{E}[\cdot]$ as a shorthand for $E_{Z_{n}\sim
Q(Z_{n}|\mu_{n},\kappa_{n})}[\cdot]$, and where

$\begin{split}\mathbbm{E}\left[e^{jZ_{n}}\right]&=\begin{pmatrix}b\left(\kappa_{n,1}\right)e^{j\mu_{n,1}}\ldots
b\left(\kappa_{n,P}\right)e^{j\mu_{n,P}}\end{pmatrix}^{\top}\\\
&=\mathrm{Diag}\left(B(\kappa_{n})\right)e^{j\mu_{n}}\end{split}$ (58)

is the $P\times 1$ mean vector and

$\begin{split}\mathrm{Cov}(e^{jZ_{n}})&=\mathrm{Diag}\left(1-b(\kappa_{n,1})^{2},...,1-b(\kappa_{n,P})^{2}\right)\\\
&=\mathrm{Diag}\left(\mathbbm{1}-B^{2}(\kappa_{n})\right)\end{split}$ (59)

is the $P\times P$ the covariance matrix of the phasor $e^{jZ_{n}}$ for
$Z_{n}\sim Q(Z_{n}|\mu_{n},\kappa_{n})$ [45]. Note that (57) follows from
(50), with $Z_{n}\sim Q(Z_{n}|\mu_{n},\kappa_{n})$ in lieu of $Z_{n}\sim
P(Z_{n}|\kappa_{0})$. Moreover, since the covariance matrix
$\mathrm{Cov}(e^{jZ_{n}})$ is diagonal, the trace term in (57) can be
efficiently computed as

$\begin{split}\operatorname{Tr}\left\\{\hat{G}(c_{n},\theta)\mathrm{Cov}(e^{jZ_{n}}){\hat{G}(c_{n},\theta)}^{\dagger}\right\\}&=L\left(\left|\alpha(\theta)\right|^{2}\right)^{\top}\left(\mathbbm{1}-B^{2}(\kappa_{n})\right).\end{split}$
(60)

All in all, the free energy in (35) can be expressed as

$\begin{split}\mathrm{F}_{n}\left(\mu_{n},\kappa_{n},\theta,\kappa_{0}\right)=&\sum_{p=1}^{P}\left\\{\log\left(\frac{I_{0}(\kappa_{0})}{I_{0}(\kappa_{n,p})}\right)+b(\kappa_{n,p})\left(\kappa_{n,p}-\kappa_{0}\cos(\mu_{n,p})\right)\right\\}+L\log\left(\pi\sigma^{2}\right)\\\
&+\frac{1}{\sigma^{2}}\left[\left\lVert\hat{G}(c_{n},\theta)\mathrm{Diag}\left(B(\kappa_{n})\right)e^{j\mu_{n}}-H_{n}\right\rVert^{2}+L\left(\left|\alpha(\theta)\right|^{2}\right)^{\top}\left(\mathbbm{1}-B^{2}(\kappa_{n})\right)\right].\end{split}$
(61)

##  Appendix C Derivation of the VEM Updates

In this appendix, we derive an approximation of a stationary point of the free
energy to be used as closed-form updates for the E-step of the VEM algorithm
presented in Sec. IV. In what follows, we adopt a Jacobian convention for
partial derivatives, such that, for
$y=(y_{1},...,y_{N_{y}})^{\top}\in\mathbbm{C}^{N_{y}\times 1}$ and
$x=(x_{1},...,x_{N_{x}})^{\top}\in\mathbbm{R}^{N_{x}\times 1}$,
$\frac{\partial y_{i}}{\partial x_{j}}$ is the element at the $i\mbox{-th}$
row and $j\mbox{-th}$ column of the Jacobian matrix $\frac{\partial
y}{\partial x}\in\mathbbm{C}^{N_{y}\times N_{x}}$ for $i\in\\{1,...,N_{y}\\}$
and $j\in\\{1,...,N_{x}\\}$. We first evaluate the gradients with respect to
the variational parameters $\mu_{n}$ and $\kappa_{n}$ for $n\in\\{1,...,N\\}$,
and then elaborate on the approximation of a stationary point with respect to
these parameters.

###  C-A Gradients

####  C-A1 Gradient for the Variational Means

From the variational free energy expression in (35), with fixed parameters
$\theta$, we evaluate the gradient with respect to the variational mean
parameter vector $\mu_{n}$ as

$\begin{split}\frac{\partial\mathrm{F}_{n}\left(\mu_{n},\kappa_{n},\theta,\kappa_{0}\right)}{\partial\mu_{n}}&=-\frac{\partial\sum_{p=1}^{P}b(\kappa_{n,p})\kappa_{0}\cos(\mu_{n,p})}{\partial\mu_{n}}+\frac{1}{\sigma^{2}}\frac{\partial\lVert\hat{G}(c_{n},\theta)\mathrm{Diag}\left(B(\kappa_{n})\right)e^{j\mu_{n}}-H_{n}\rVert^{2}}{\partial\mu_{n}}\\\
&=\Im\left\\{C(\mu_{n},\kappa_{n})\right\\}\mathrm{Diag}\left(B(\kappa_{n})\right),\end{split}$
(62)

with $1\times P$ row-vector

$\begin{split}C(\mu_{n},\kappa_{n})=&\kappa_{0}\left(e^{j\mu_{n}}\right)^{\top}+\frac{2}{\sigma^{2}}{H_{n}}^{\dagger}\hat{G}(c_{n},\theta)\mathrm{Diag}\left(e^{j\mu_{n}}\right)\\\
&-\frac{2}{\sigma^{2}}B(\kappa_{n})^{\top}{\left(\hat{G}(c_{n},\theta)\mathrm{Diag}\left(e^{j\mu_{n}}\right)\right)}^{\dagger}\left(\hat{G}(c_{n},\theta)\mathrm{Diag}\left(e^{j\mu_{n}}\right)\right),\end{split}$
(63)

where $\Im\\{\cdot\\}$ denotes the element-wise imaginary part of a complex
vector or matrix.

####  C-A2 Gradient for the Variational Concentration Parameters

The free energy gradient with respect to the concentration parameter vector
$\kappa_{n}$ is decomposed as

$\begin{split}\frac{\partial\mathrm{F}_{n}\left(\mu_{n},\kappa_{n},\theta,\kappa_{0}\right)}{\partial\kappa_{n}}&=-\frac{\partial\sum_{p=1}^{P}\log{I_{0}(\kappa_{n,p})}}{\partial\kappa_{n}}+\frac{\partial\sum_{p=1}^{P}b(\kappa_{n,p})\left(\kappa_{n,p}-\kappa_{0}\cos(\mu_{n,p})\right)}{\partial\kappa_{n}}+\\\
&\quad+\frac{1}{\sigma^{2}}\left[\frac{\partial\lVert\hat{G}(c_{n},\theta)\mathrm{Diag}\left(e^{j\mu_{n}}\right)B(\kappa_{n})-H_{n}\rVert^{2}}{\partial\kappa_{n}}+L\left(\left|\alpha(\theta)\right|^{2}\right)^{\top}\frac{\partial\left(\mathbbm{1}-B^{2}(\kappa_{n})\right)}{\partial\kappa_{n}}\right]\\\
&=\left(\kappa_{n}^{\top}-\frac{2L}{\sigma^{2}}B(\kappa_{n})^{\top}\mathrm{Diag}\left(\left|\alpha(\theta)\right|^{2}\right)-\Re\left\\{C(\mu_{n},\kappa_{n})\right\\}\right)\frac{\partial
B(\kappa_{n})}{\partial\kappa_{n}},\end{split}$ (64)

with row-vector $C(\mu_{n},\kappa_{n})$ defined in (63).

###  C-B Approximating the Stationary Points via Phasor Relaxation

In what follows, we address the subset of stationary points
$(\mu_{n}^{*},\kappa_{n}^{*})$ for the variational parameters that satisfy the
equality $C(\mu_{n}^{*},\kappa_{n}^{*})=0$. Note that this may be a strict
subset of the set of all stationary points, and that it may potentially be
empty. To proceed, we relax the problem by optimizing over an arbitrary
$P\times 1$ complex vector $\rho_{n}$ in lieu of the phasor $e^{j\mu_{n}}$ in
(63). The stationarity condition of mean gradient (62) is verified for all
solutions satisfying the equality $C(\rho_{n}^{*},\kappa_{n}^{*})=0$. Note
that we have replaced $\mu_{n}$ by $\rho_{n}$ in order to highlight the
adopted change of variables for optimization. Considering also the
stationarity condition for the concentration gradient (64), we obtain the
system of equations

$\left\\{\begin{array}[]{ll}C(\rho_{n}^{*},\kappa_{n}^{*})=0\\\
\left(\kappa_{n}^{*}-\frac{2L}{\sigma^{2}}\mathrm{Diag}\left(\left|\alpha(\theta)\right|^{2}\right)B(\kappa_{n}^{*})\right)\frac{\partial
B(\kappa_{n})}{\partial\kappa_{n}}\biggr{\rvert}_{\kappa_{n}^{*}}=0.\end{array}\right.$
(65)

Noting that the elements in the main diagonal of the matrix $\partial
B(\kappa_{n})/\partial\kappa_{n}$ are all non-null, the second equation in
(65) simplifies into $P$ independent scalar equations
$f_{p}(\kappa_{n,p}^{*})=0$, with

$f_{p}(\kappa_{n,p})=\kappa_{n,p}-\frac{2L}{\sigma^{2}}\left|\alpha_{p}(\theta)\right|^{2}b(\kappa_{n,p}),$
(66)

for $p\in\\{1,...,P\\}$. While $\kappa_{n,p}^{*}=0$ is a trivial solution to
$f_{p}(\kappa_{n,p}^{*})=0$, a strictly positive solution $\kappa_{n,p}^{*}>0$
can be found if and only if we have

$\frac{L|\alpha_{p}|^{2}}{\sigma^{2}}>1.$ (67)

This statement can be proven by studying the sign of $f_{p}(\kappa_{n,p})$ for
$\kappa_{n,p}\in[0,+\infty)$. When inequality (67) is verified, a solution
$\kappa_{n,p}^{*}>0$ to $f_{p}(\kappa_{n,p}^{*})=0$ can be bounded using the
Bessel ratio bounds in [50] as

$2\sqrt{\frac{L|\alpha_{p}|^{2}}{\sigma^{2}}-1}\sqrt{\frac{L|\alpha_{p}|^{2}}{\sigma^{2}}}<\kappa_{n,p}^{*}<2\sqrt{\frac{L|\alpha_{p}|^{2}}{\sigma^{2}}-\frac{1}{2}}\sqrt{\frac{L|\alpha_{p}|^{2}}{\sigma^{2}}},$
(68)

for $p\in\\{1,...,P\\}$. In practice, we use the lower bound in (68) to
approximate $\kappa_{n,p}^{*}$ for $p\in\\{1,...,P\\}$. Note that a tighter
bound in (68) can be obtained using sharper bounds of the Bessel ratio in
[50].

Furthermore, assuming the $L\times P$ matrix $\hat{G}(c_{n},\theta)$ is of
rank $P$, the first equation in (65) yields

$\mathrm{Diag}\left(B(\kappa_{n}^{*})\right)\rho^{*}_{n}=\left({\hat{G}(c_{n},\theta)}^{\dagger}\hat{G}(c_{n},\theta)\right)^{-1}\left(\frac{\sigma^{2}\kappa_{0}}{2}\mathbbm{1}+{\hat{G}(c_{n},\theta)}^{\dagger}H_{n}\right),$
(69)

where existence of the inverse
$({\hat{G}(c_{n},\theta)}^{\dagger}\hat{G}(c_{n},\theta))^{-1}$ stems from the
full-column rank assumption on $\hat{G}(c_{n},\theta)$. By (69), if all
entries in $B(\kappa_{n}^{*})$ are strictly positive, the element-wise phase
of vector $\rho^{*}_{n}$ satisfies

$\angle\left\\{\rho^{*}_{n}\right\\}=\angle\left\\{\left({\hat{G}(c_{n},\theta)}^{\dagger}\hat{G}(c_{n},\theta)\right)^{-1}\left(\frac{\sigma^{2}\kappa_{0}}{2}\mathbbm{1}+{\hat{G}(c_{n},\theta)}^{\dagger}H_{n}\right)\right\\}.$
(70)

Finally, we project this solution back into the space of phasors, producing
the proposed approximate solution (38).

##  Appendix D Optimizing the Prior Concentration Parameter

From (43), given that the function
$\sum_{n=1}^{N}\mathrm{F}_{n}(\mu_{n},\kappa_{n},\theta,\kappa_{0})$ is convex
with respect to $\kappa_{0}$, an optimal solution $\kappa_{0}^{*}$ can be
computed by evaluating a stationary point $\kappa_{0}^{*}$ as

$\frac{\partial\sum_{n=1}^{N}\mathrm{F}_{n}\left(\mu_{n},\kappa_{n},\theta,\kappa_{0}\right)}{\partial\kappa_{0}}\biggr{\rvert}_{\kappa_{0}^{*}}=0,$
(71)

where

$\begin{split}\frac{\partial\sum_{n=1}^{N}\mathrm{F}_{n}\left(\mu_{n},\kappa_{n},\theta,\kappa_{0}\right)}{\partial\kappa_{0}}&=\sum_{n=1}^{N}\sum_{p=1}^{P}\left\\{b(\kappa_{0})-b(\kappa_{n,p})\cos(\mu_{n,p})\right\\}.\end{split}$
(72)

This yields

$\kappa_{0}^{*}=b^{-1}\left(\frac{1}{NP}\sum_{n=1}^{N}\sum_{p=1}^{P}b(\kappa_{n,p})\cos(\mu_{n,p})\right),$
(73)

where $b^{-1}(\cdot)$ denotes the inverse function of the Bessel ratio
$b(\cdot)$ [52]. Given that $b^{-1}(\cdot)$ is defined on $[0,1)$, the result
in (73) holds only if
$0\leq(NP)^{-1}\sum_{n=1}^{N}\sum_{p=1}^{P}b(\kappa_{n,p})\cos(\mu_{n,p})<1$.
While the upper bound is verified for any set of parameters
$\\{\mu_{n},\kappa_{n}\\}_{n=1}^{N}$, the lower bound can be exceeded if some
of the variational mean estimates verify $|\mu_{n,p}|>\pi/2$ for
$\mu_{n,p}\in[-\pi,\pi)$. For such cases, i.e., when
$\sum_{n=1}^{N}\sum_{p=1}^{P}b(\kappa_{n,p})\cos(\mu_{n,p})<0$, the free
energy derivative in (72) is always positive, and
$\sum_{n=1}^{N}\mathrm{F}_{n}(\mu_{n},\kappa_{n},\theta,\kappa_{0})$ becomes a
strictly increasing function with respect to $\kappa_{0}$. Accordingly, we
extend the result in (73) to $\kappa_{0}^{*}=0$ when $\\\
(NP)^{-1}\sum_{n=1}^{N}\sum_{p=1}^{P}b(\kappa_{n,p})\cos(\mu_{n,p})<0$.

## References

  * [1] A. Lavin, D. Krakauer, H. Zenil, J. Gottschlich, T. Mattson, J. Brehmer, A. Anandkumar, S. Choudry, K. Rocki, A. G. Baydin _et al._ , “Simulation intelligence: Towards a new generation of scientific methods,” _arXiv preprint arXiv:2112.03235_ , 2021.
  * [2] W. Kritzinger, M. Karner, G. Traar, J. Henjes, and W. Sihn, “Digital twin in manufacturing: A categorical literature review and classification,” _IFAC-PapersOnLine_ , vol. 51, no. 11, pp. 1016–1022, 2018.
  * [3] I. Errandonea, S. Beltrán, and S. Arrizabalaga, “Digital twin for maintenance: A literature review,” _Computers in Industry_ , vol. 123, p. 103316, 2020.
  * [4] M. Liu, S. Fang, H. Dong, and C. Xu, “Review of digital twin about concepts, technologies, and industrial applications,” _Journal of Manufacturing Systems_ , vol. 58, pp. 346–361, 2021.
  * [5] A. Thelen, X. Zhang, O. Fink, Y. Lu, S. Ghosh, B. D. Youn, M. D. Todd, S. Mahadevan, C. Hu, and Z. Hu, “A comprehensive review of digital twin—part 1: modeling and twinning enabling technologies,” _Structural and Multidisciplinary Optimization_ , vol. 65, no. 12, pp. 1–55, 2022.
  * [6] ——, “A comprehensive review of digital twin—part 2: roles of uncertainty quantification and optimization, a battery digital twin, and perspectives,” _Structural and multidisciplinary optimization_ , vol. 66, no. 1, p. 1, 2023.
  * [7] M. Polese, L. Bonati, S. D’oro, S. Basagni, and T. Melodia, “Understanding O-RAN: Architecture, interfaces, algorithms, security, and research challenges,” _IEEE Communications Surveys & Tutorials_, 2023.
  * [8] L. U. Khan, W. Saad, D. Niyato, Z. Han, and C. S. Hong, “Digital-twin-enabled 6G: Vision, architectural trends, and future directions,” _IEEE Communications Magazine_ , vol. 60, no. 1, pp. 74–80, 2022.
  * [9] J. Jagannath, K. Ramezanpour, and A. Jagannath, “Digital twin virtualization with machine learning for IoT and beyond 5G networks: Research directions for security and optimal control,” in _Proceedings of the 2022 ACM Workshop on Wireless Security and Machine Learning_ , 2022, pp. 81–86.
  * [10] L. Hui, M. Wang, L. Zhang, L. Lu, and Y. Cui, “Digital twin for networking: A data-driven performance modeling perspective,” _IEEE Network_ , 2022.
  * [11] B. Salehi, U. Demir, D. Roy, S. Pradhan, J. Dy, S. Ioannidis, and K. Chowdhury, “Multiverse at the edge: Interacting real world and digital twins for wireless beamforming,” _arXiv preprint arXiv:2305.10350_ , 2023.
  * [12] D. Villa, M. Tehrani-Moayyed, C. P. Robinson, L. Bonati, P. Johari, M. Polese, S. Basagni, and T. Melodia, “Colosseum as a digital twin: Bridging real-world experimentation and wireless network emulation,” _arXiv preprint arXiv:2303.17063_ , 2023.
  * [13] J. Hoydis, F. A. Aoudia, S. Cammerer, M. Nimier-David, N. Binder, G. Marcus, and A. Keller, “Sionna RT: Differentiable ray tracing for radio propagation modeling,” _arXiv preprint arXiv:2303.11103_ , 2023.
  * [14] M. Pengnoo, M. T. Barros, L. Wuttisittikulkij, B. Butler, A. Davy, and S. Balasubramaniam, “Digital twin for metasurface reflector management in 6g terahertz communications,” _IEEE access_ , vol. 8, pp. 114 580–114 596, 2020.
  * [15] N. Zeulin, A. Ponomarenko-Timofeev, O. Galinina, and S. Andreev, “ML-assisted beam selection via digital twins for time-sensitive industrial IoT,” _IEEE Internet of Things Magazine_ , vol. 5, no. 1, pp. 36–40, 2022.
  * [16] S. Jiang and A. Alkhateeb, “Digital twin based beam prediction: Can we train in the digital world and deploy in reality?” _arXiv preprint arXiv:2301.07682_ , 2023.
  * [17] A. Alkhateeb, S. Jiang, and G. Charan, “Real-time digital twins: Vision and research directions for 6G and beyond,” _IEEE Communications Magazine_ , 2023.
  * [18] V.-H. Nguyen, V. Corlay, N. Gresset, and C. Ciochina, “Probabilistic ray-tracing aided positioning at mmwave frequencies,” _arXiv preprint arXiv:2308.08441_ , 2023.
  * [19] A. F. Molisch, _Wireless communications_. John Wiley & Sons, 2012.
  * [20] ITU, “Recommendation ITU-R P. 526-9: Propagation by diffraction,” 2019.
  * [21] V. Degli-Esposti, F. Fuschini, E. M. Vitucci, and G. Falciasecca, “Measurement and modelling of scattering from buildings,” _IEEE Transactions on Antennas and Propagation_ , vol. 55, no. 1, pp. 143–153, 2007.
  * [22] V. Degli-Esposti, V.-M. Kolmonen, E. M. Vitucci, and P. Vainikainen, “Analysis and modeling on co-and cross-polarized urban radio propagation for dual-polarized MIMO wireless systems,” _IEEE transactions on antennas and propagation_ , vol. 59, no. 11, pp. 4247–4256, 2011.
  * [23] K. Vuckovic, M. B. Mashhadi, F. Hejazi, N. Rahnavard, and A. Alkhateeb, “PARAMOUNT: Towards generalizable deep learning for mmwave beam selection using sub-6GHz channel measurements,” _IEEE Transactions on Wireless Communications_ , 2023.
  * [24] H. Viswanathan and P. E. Mogensen, “Communications in the 6G era,” _IEEE Access_ , vol. 8, pp. 57 063–57 074, 2020.
  * [25] P. Ferrand, M. Guillaud, C. Studer, and O. Tirkkonen, “Wireless channel charting: Theory, practice, and applications,” _IEEE Communications Magazine_ , vol. 61, no. 6, pp. 124–130, 2023.
  * [26] C. Ruah, O. Simeone, and B. Al-Hashimi, “A bayesian framework for digital twin-based control, monitoring, and data collection in wireless systems,” _IEEE Journal on Selected Areas in Communications_ , 2023.
  * [27] OpenStreetMap contributors, “Planet dump retrieved from https://planet.osm.org ,” https://www.openstreetmap.org, 2022.
  * [28] ITU, “Recommendation ITU-R P.2040-2: Effects of building materials and structures on radiowave propagation above about 100 MHz,” 2021.
  * [29] J. Jemai, P. C. Eggers, G. F. Pedersen, and T. Kurner, “Calibration of a UWB sub-band channel model using simulated annealing,” _IEEE Transactions on Antennas and Propagation_ , vol. 57, no. 10, pp. 3439–3443, 2009.
  * [30] S. Priebe, M. Jacob, and T. Kürner, “Calibrated broadband ray tracing for the simulation of wave propagation in mm and sub-mm wave indoor communication channels,” in _European Wireless 2012; 18th European Wireless Conference 2012_. VDE, 2012, pp. 1–10.
  * [31] D. He, B. Ai, K. Guan, L. Wang, Z. Zhong, and T. Kürner, “The design and applications of high-performance ray-tracing simulation platform for 5G and beyond wireless communications: A tutorial,” _IEEE communications surveys & tutorials_, vol. 21, no. 1, pp. 10–27, 2018.
  * [32] R. Charbonnier, C. Lai, T. Tenoux, D. Caudill, G. Gougeon, J. Senic, C. Gentile, Y. Corre, J. Chuang, and N. Golmie, “Calibration of ray-tracing with diffuse scattering against 28-GHz directional urban channel measurements,” _IEEE Transactions on Vehicular Technology_ , vol. 69, no. 12, pp. 14 264–14 276, 2020.
  * [33] O. Kanhere and T. S. Rappaport, “Calibration of NYURay, a 3D mmwave and sub-THz ray tracer using indoor, outdoor, and factory channel measurements,” _arXiv preprint arXiv:2302.12380_ , 2023.
  * [34] G. S. Bhatia, Y. Corre, and M. Di Renzo, “Tuning of ray-based channel model for 5G indoor industrial scenarios,” in _2023 IEEE International Mediterranean Conference on Communications and Networking (MeditCom)_. IEEE, 2023, pp. 311–316.
  * [35] F. Xue, W. Lu, Z. Chen, and C. J. Webster, “From LiDAR point cloud towards digital twin city: Clustering city objects based on Gestalt principles,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 167, pp. 418–431, 2020.
  * [36] X. Lin, L. Kundu, C. Dick, E. Obiodu, T. Mostak, and M. Flaxman, “6G digital twin networks: From theory to practice,” _IEEE Communications Magazine_ , 2023.
  * [37] Z. Li, T. Müller, A. Evans, R. H. Taylor, M. Unberath, M.-Y. Liu, and C.-H. Lin, “Neuralangelo: High-fidelity neural surface reconstruction,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2023, pp. 8456–8465.
  * [38] O. Simeone, _Machine Learning for Engineers_. Cambridge University Press, 2022.
  * [39] J. Hoydis, S. Cammerer, F. A. Aoudia, A. Vem, N. Binder, G. Marcus, and A. Keller, “Sionna: An open-source library for next-generation physical layer research,” _arXiv preprint arXiv:2203.11854_ , 2022.
  * [40] C. A. Balanis, _Antenna theory: analysis and design_. John wiley & sons, 2016.
  * [41] H. Do, N. Lee, and A. Lozano, “Parabolic wavefront model for line-of-sight MIMO channels,” _IEEE Transactions on Wireless Communications_ , 2023.
  * [42] C. F. Van Loan, “The ubiquitous Kronecker product,” _Journal of computational and applied mathematics_ , vol. 123, no. 1-2, pp. 85–100, 2000.
  * [43] M.-A. Badiu and J. P. Coon, “Communication through a large reflecting surface with phase errors,” _IEEE Wireless Communications Letters_ , vol. 9, no. 2, pp. 184–188, 2019.
  * [44] T. Wang, M.-A. Badiu, G. Chen, and J. P. Coon, “Outage probability analysis of RIS-assisted wireless networks with Von Mises phase errors,” _IEEE Wireless Communications Letters_ , vol. 10, no. 12, pp. 2737–2741, 2021.
  * [45] K. V. Mardia and P. E. Jupp, _Directional statistics_. Wiley Online Library, 2000, vol. 2.
  * [46] S. Royo and M. Ballesta-Garcia, “An overview of LiDAR imaging systems for autonomous vehicles,” _Applied sciences_ , vol. 9, no. 19, p. 4093, 2019.
  * [47] M. D. Bedford, A. Hrovat, G. Kennedy, T. Javornik, and P. Foster, “Modeling microwave propagation in natural caves using LiDAR and ray tracing,” _IEEE Transactions on Antennas and Propagation_ , vol. 68, no. 5, pp. 3878–3888, 2019.
  * [48] B. Keinert, M. Innmann, M. Sänger, and M. Stamminger, “Spherical Fibonacci mapping,” _ACM Transactions on Graphics (TOG)_ , vol. 34, no. 6, pp. 1–7, 2015.
  * [49] C. M. Bishop and N. M. Nasrabadi, _Pattern recognition and machine learning_. Springer, 2006, vol. 4, no. 4.
  * [50] D. Ruiz-Antolín and J. Segura, “A new type of sharp bounds for ratios of modified Bessel functions,” _Journal of Mathematical Analysis and Applications_ , vol. 443, no. 2, pp. 1232–1246, 2016.
  * [51] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, _Introduction to algorithms_. MIT press, 2022.
  * [52] G. W. Hill, “Evaluation and inversion of the ratios of modified bessel functions, $I_{1}(x)/I_{0}(x)$ and $I_{1.5}(x)/I_{0.5}(x)$,” _ACM Transactions on Mathematical Software (TOMS)_ , vol. 7, no. 2, pp. 199–208, 1981.
  * [53] C. Ruah, “Code Repository for Ray Tracer Calibration using Local Phase Error Estimates,” 2023. [Online]. Available: https://github.com/kclip/phase-aware-rt-calibration
  * [54] J. Hoydis, F. A. Aoudia, S. Cammerer, F. Euchner, M. Nimier-David, S. ten Brink, and A. Keller, “Learning radio environments by differentiable ray tracing,” _arXiv preprint arXiv:2311.18558_ , 2023.
  * [55] S. Lang, _Linear algebra_ , 3rd ed. Springer, 1987.
  * [56] U. Lund and S. Rao Jammalamadaka, “An entropy-based test for goodness of fit of the von Mises distribution,” _Journal of statistical computation and simulation_ , vol. 67, no. 4, pp. 319–332, 2000.