text:
\section{Introduction}
MmWave communications in 5G and beyond are useful for \acf{slam} applications, due to geometrical propagation channels, large bandwidths, and large antenna arrays \cite{nurmi2017multi}. Signals sent from the \acf{bs} reach the \acf{ue} via the propagation channel, which is determined by the geometric relationships among the propagation environment, the \ac{ue}, and the \ac{bs}. Large bandwidths and antenna arrays result in high temporal and  spatial resolutions \cite{patwari2005,larsen2009}. Therefore, state-of-the-art channel estimators can provide accurate estimates for multipath components by using the received signals, in terms of groups of channel gain, \ac{TOA}, \ac{AOA}, \ac{AOD}, and Doppler, which contain  the necessary information for \ac{slam} \cite{wymeersch20175g}. Although Doppler contains geometric information,
it is usually ignored.

The related works can be divided into two areas: works that exploit Doppler for radio positioning or mapping and works in the area of radio SLAM. Doppler has been used in most radars for mapping and tracking \cite{patole2017automotive}, but limited works have been done in radio scenarios. Doppler is used to localize radio emitters in \cite{amar2008localization}, but the proposed method can only be used for narrow-band signals, while \cite{Han2016} shows that the Doppler shift can provide more direction information for localization, and \cite{kakkavas2019performance} shows that the mobility can significantly improve the \ac{nlos}-only scenario in MIMO mmWave system,  indicating the involvement of the Doppler shift brings gain in the localization accuracy. Doppler is used for tracking \acp{ue}' positions and velocities in a Wi-Fi-based system in  \cite{qian2017widar}. However, these methods do not solve the \ac{slam} problem.
Several approaches have been proposed to address the mmWave radio \ac{slam} problem, including geometry-based methods \cite{wen20195g, yassin2018mosaic}, message passing-based methods \cite{Erik_BPSLAM_TWC2019,Rico_BPSLAM_JSTSP2019}, and \ac{rfs}-based methods \cite{kim20205g,ge20205GSLAM,ge2022computationally}. \ac{rfs}-based methods
can  handle uncertainties, as well as inherently deal with challenges of the unknown number of landmarks, unknown \acp{DA}, misdetections, and clutter measurements in radio \ac{slam}. Within these \ac{rfs}-based methods, the probability hypothesis density (PHD) filter is used in \cite{kim20205g}, which does not have an explicit enumeration of \acp{DA}, and the \ac{PMBM} filter is used in \cite{ge20205GSLAM,ge2022computationally}, which enumerates all possible DAs explicitly, thus allowing for improved performance. However, the Doppler is not considered in these works. To our best knowledge, the inclusion of Doppler in radio SLAM has not yet been conducted in the existing literature.

In this paper, we harness the Doppler component in bistatic mmWave radio \ac{slam} and analyze how Doppler benefits the \ac{slam} filter. The main contributions of this paper are summarized as follows: (\emph{i}) we derive the \ac{PEB}, \ac{HEB}, \ac{CEB} of the \ac{ue} and the \acp{LEB} of landmarks by computing the \acf{PCRB} for the system with utilizing Doppler;  (\emph{ii}) we analyze the effect of the quality of the Doppler measurement on the bounds of the \ac{ue} and the average \ac{LEB} of landmarks; (\emph{iii}) we extend our previous work in \cite{ge2022computationally} by involving the Doppler as a dimension of the measurement in the  \ac{ek}-\ac{PMB} \ac{slam} filter, and validate the benefits of exploiting Doppler in the system through numerical experiments in the mmWave network context.

\subsubsection*{Notations}
Scalars (e.g., $x$) are denoted in italic, vectors (e.g., $\boldsymbol{x}$) in bold lower-case letters  with $\left\Vert\boldsymbol{x}\right\Vert$ representing its L2-norm, matrices (e.g., $\boldsymbol{X}$) in bold capital letters, sets  (e.g., $\mathcal{X}$) in calligraphic with $\left|\mathcal{X}\right|$ representing its cardinality. The transpose is denoted by $(\cdot)^{\mathsf{T}}$, the Hermitian transpose is denoted by $(\cdot)^{\mathsf{H}}$, the union of mutually disjoint sets is denoted by $\uplus$, the Kronecker product is denoted by $\otimes$, the expectation is denoted by $\mathbb{E}[\cdot]$, $\mathcal{N}(\boldsymbol{u},\boldsymbol{\Sigma})$ denotes a multivariate Gaussian distribution with mean $\boldsymbol{u}$ and covariance $\boldsymbol{\Sigma}$, and $d_{\boldsymbol{x}}=\text{dim}(\boldsymbol{x})$.

\vspace{-1mm}
\section{System model}
In this section, the \ac{ue} model, the environment model, the signal model, and the measurement model for a mmWave radio downlink scenario as shown in Fig.\,\ref{fig:scenario} are briefly introduced.

\begin{figure}
\centering
    \includegraphics[width=0.98\linewidth]{scenario_details_1.pdf}
\caption{A mmWave downlink scenario with the environment of a \ac{bs}, a \ac{ue}, a reflecting surface, and a small object. The reflecting surface can be modeled as a \acf{va}, which is the reflection of the \ac{bs} with respect to the surface. The small object is modeled as a \acf{sp}. The \ac{bs} sends signals to the \ac{ue} via \acf{los} path and/or \ac{nlos} paths, shown as the blue line and  purple lines, respectively. The channel parameters of each path depend on the underlying geometry.}
\label{fig:scenario}
\end{figure}

\subsection{State Models}
In this paper, a single-user scenario is considered, where the \ac{ue} does a constant turn-rate movement around a known \ac{bs} on the $x-y$ plane.
We denote the dynamic state of the \ac{ue} at time step $k$ as $\boldsymbol{s}_{k}=[x_{k},y_{k},\varpi_{k},b_{k}]^{\mathsf{T}}$, containing the \ac{ue} position on $x$- and $y$-axis, heading, and clock bias, respectively. The \ac{ue} dynamics can be expressed as \cite{roth2014}
\begin{equation}\label{eq:dynamic_model}
     \underset{\boldsymbol{s}_{k}}{\underbrace{\begin{bmatrix} x_{k} \\ y_{k} \\ \varpi_{k} \\ b_{k} \end{bmatrix}}} = \underset{\boldsymbol{f}\left(\boldsymbol{s}_{k-1} \right)}{\underbrace{\begin{bmatrix} x_{k-1} + \tfrac{2v}{\omega}  \sin \! \left(\tfrac{\omega T}{2}\right)  \cos \! \left(\varpi_{k-1} + \tfrac{\omega T}{2}\right) \\ y_{k-1} + \tfrac{2v}{\omega}\sin \! \left(\tfrac{\omega T}{2}\right) \sin \! \left(\varpi_{k-1} + \tfrac{\omega T}{2}\right) \\ \varpi_{k-1} + \omega T \\ b_{k-1} \end{bmatrix}}} + \boldsymbol{\eta}_{k},
\end{equation}
where $v$, $\omega$, and $T$ are known control inputs, denoting the speed,  the turn-rate,  and the sampling interval, respectively, and  $\boldsymbol{\eta}_k \sim \mathcal{N}\left(\boldsymbol{0},\boldsymbol{Q}\right)$.

The considered scenario contains a single known \ac{bs}, and multiple unknown reflecting surfaces and small objects, which are modeled as \acp{va} and \acp{sp}, respectively. The \ac{bs} sends downlink signals to the \ac{ue}, which can reach the \ac{ue} directly, termed as the \ac{los} path, and/or can be reflected by the reflecting surfaces or scattered by the small objects and reach the  \ac{ue}, termed as \ac{nlos} paths  (see Fig.\,\ref{fig:scenario}). In this paper, we assume that there is at most one path associated to each landmark every time step. The landmark state  can be modeled as $\boldsymbol{x} = [\boldsymbol{x}^{\textsf{T}}_{\text{LM}},m]^{\textsf{T}}$,  with $\boldsymbol{x}_{\text{LM}} \in \mathbb{R}^{3}$ denoting the landmark location, and $m\in\{\text{BS},\text{VA},\text{SP}\}$ denoting the landmark type. Therefore, the map of the environment can be modeled by a set of landmark $\mathcal{X}=\{\boldsymbol{x}^{1},\dots, \boldsymbol{x}^{I}\}$, with ${I}$ representing the total number of landmarks.

\subsection{Signal Model}
\label{Signal model}
The \ac{bs} sends downlink signals to the \ac{ue} at every time step with a period of $T$ seconds. These signals can reach the \ac{ue} via \ac{los} path and/or \ac{nlos} paths, and the received signal for OFDM symbol $n$, at subcarrier $\kappa$ and time step $k$ can be expressed as~\cite{heath2016overview}
\begin{align}
    \boldsymbol{y}_{\kappa,n,k} =\boldsymbol{W}_{n,k}^{\mathsf{H}}&\sum _{i=1}^{I_{k}}g_{k}^{i}\boldsymbol{a}_{\text{R}}(\boldsymbol{\theta}_{k}^{i})\boldsymbol{a}_{\text{T}}^{\mathsf{H}}(\boldsymbol{\phi}_{k}^{i}) e^{-\jmath 2\pi \kappa \Delta f \tau_{k}^{i}}\nonumber\\&\times e^{\jmath 2\pi n T_{\text{mea}}  d_{k}^{i}/\varsigma} \boldsymbol{p}_{n,k} + \boldsymbol{W}_{n,k}^{\mathsf{H}}\boldsymbol{n}_{\kappa,n,k}, \label{eq:FreqObservation}
\end{align}
where $\boldsymbol{y}_{\kappa,n,k}$ denotes the received signal, $\boldsymbol{p}_{n,k}$ denotes the pre-coded pilot signal,  $\boldsymbol{n}_{\kappa,n,k}$ denotes a white Gaussian noise, $\boldsymbol{W}_{n,k}$ denotes a combining matrix, $I_{k}$ denotes the number of all visible landmarks at time $k$, $\boldsymbol{a}_{\text{R}}(\cdot)$ and $\boldsymbol{a}_{\text{T}}(\cdot)$ denote the steering vectors of the receiver and transmitter antenna arrays, respectively, $\Delta f$ denotes the subcarrier spacing, $T_{\text{mea}}$ denotes transmission duration of each symbol, and $\varsigma$ denotes the wavelength.
Each path $i$ can be described by a complex gain $g_{k}^{i}$, a \ac{TOA} $\tau_{k}^{i}$, an \ac{AOA} pair $\boldsymbol{\theta}_{k}^{i}$ in azimuth and elevation,  an \ac{AOD} pair $\boldsymbol{\phi}_{k}^{i}$ in azimuth and elevation, and a Doppler $d_{k}^{i}$ which we express in [m/s]. These parameters can be  estimated by a parametric
channel estimation algorithm from $\boldsymbol{y}_{\kappa,n,k}$, for example, by \cite{richter2005estimation,alkhateeb2014channel,jiang2021high}, which is out of the scope of this paper, and the estimation results are utilized directly. Adding Doppler can make paths more resolvable, as there may be some paths not resolvable in delay or angle domains, but resolvable in the Doppler domain.

\subsection{Measurement Model}
At time step $k$, the channel estimator provides a set of measurements $\mathcal{Z}_{k}=\{\boldsymbol{z}_{k}^{1},\dots, \boldsymbol{z}_{k}^{\hat{{I}}_{k}} \}$, where usually $\hat{{I}}_{k} \neq {{I}}_{k}$, as some visible landmarks may be misdetected and there could be some clutter measurements. Please note that the source of each element in $\mathcal{Z}_{k}$ is unknown.
If the measurement noise is zero-mean Gaussian, the measurement originating from landmark $\boldsymbol{x}^{i}$ at time step $k$ can be described as follows
\begin{align}
    \boldsymbol{z}^{i}_{k}=\boldsymbol{h}(\boldsymbol{x}^{i},\boldsymbol{s}_{k})+\boldsymbol{\epsilon}_k^i,\label{pos_to_channelestimation}
\end{align}

where
$\boldsymbol{h}(\boldsymbol{x}^{i},\boldsymbol{s}_{k})=[\tau_{k}^{i},(\boldsymbol{\theta}_{k}^{i})^{\mathsf{T}},(\boldsymbol{\phi}_{k}^{i})^{\mathsf{T}}, d_{k}^{i}]^{\mathsf{T}}$ represents the nonlinear function that transforms the geometric information to the channel parameters, and $\boldsymbol{\epsilon}_k \sim \mathcal{N}\left(\boldsymbol{0},\boldsymbol{R}_{k}^{i}\right)$

These channel parameters depend on the geometric relationships among the \ac{bs}, the \ac{ue} and the landmarks. Specifically, \ac{TOA} $\tau_{k}^{i}$ can be defined as
\begin{align}
    \tau_{k}^{i} = \begin{cases} \Vert\boldsymbol{x}_{\text{BS}}-\boldsymbol{x}_{\text{UE},k}\Vert/c + b_{k} & m^{i}=\text{BS} \\  (\Vert\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{UE},k}\Vert + \Vert\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{BS}}\Vert)/c +b_{k} & m^{i}\neq\text{BS},\end{cases} \label{delay}
\end{align}
where $\boldsymbol{x}_{\text{UE},k}= [x_{k},y_{k},0]^{\textsf{T}}$ denotes the 3D position of the \ac{ue}, $\boldsymbol{x}_{\text{inc},k}^{i}$ is the incidence point of the $i$-th path on the corresponding landmark at time step $k$, which can be determined by $\boldsymbol{x}_{\text{LM}}^{i}$ and $\boldsymbol{x}_{\text{UE},k}$, and $c$ is the speed of light. As $\boldsymbol{\theta}_{k}^{i}$ is determined by the arrival direction of the signal, which can be calculated by
\begin{align}
    \boldsymbol{q}_{\text{AOA},k}^{i} = \begin{cases} (\boldsymbol{x}_{\text{BS}}-\boldsymbol{x}_{\text{UE},k})/\Vert\boldsymbol{x}_{\text{BS}}-\boldsymbol{x}_{\text{UE},k}\Vert & m^i=\text{BS} \\ (\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{UE},k})/\Vert\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{UE},k}\Vert & m^i\neq\text{BS},\end{cases}
\end{align}
we can define $\boldsymbol{\theta}_{k}^{i}$ as
\begin{align}
    &\theta_{\text{az},k}^{i} = \mathrm{arctan2}([\boldsymbol{q}_{\text{AOA},k}^{i}]_{2},[\boldsymbol{q}_{\text{AOA},k}^{i}]_{1})-\varpi_{k},\label{AOAAZ}\\
    &\theta_{\text{el},k}^{i} = \arcsin([\boldsymbol{q}_{\text{AOA},k}^{i}]_{3},\Vert\boldsymbol{q}_{\text{AOA},k}^{i}\Vert).\label{AOAEL}
\end{align}
Similarly, the departure direction can be calculated by  \begin{align}
    \boldsymbol{q}_{\text{AOD},k}^{i} = \begin{cases} (\boldsymbol{x}_{\text{UE},k}-\boldsymbol{x}_{\text{BS}})/\Vert(\boldsymbol{x}_{\text{UE},k}-\boldsymbol{x}_{\text{BS}})\Vert & m^i=\text{BS} \\ (\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{BS}})/\Vert\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{BS}}\Vert & m^i\neq\text{BS},\end{cases}
\end{align}
and $\boldsymbol{\phi}_{k}^{i}$ can be defined as
\begin{align}
    &\phi_{\text{az},k}^{i} = \mathrm{arctan2}([\boldsymbol{q}_{\text{AOD},k}^{i}]_{2},[\boldsymbol{q}_{\text{AOD},k}^{i}]_{1}),\label{AODAZ}\\
    &\phi_{\text{el},k}^{i} = \arcsin([\boldsymbol{q}_{\text{AOD},k}^{i}]_{3},\Vert\boldsymbol{q}_{\text{AOD},k}^{i}\Vert).\label{AODEL}
\end{align}
The Doppler $d_{k}^{i}$ can be calculated by projecting the \ac{ue} velocity on the direction of $\boldsymbol{q}_{\text{AOA},k}^{i}$, as displayed in Fig.\,\ref{fig:scenario}.  It can be computed by
\begin{align}
    &d_{k}^{i} = \boldsymbol{v}_{k}^{\textsf{T}}\boldsymbol{q}_{\text{AOA},k}^{i},\label{Doppler}
\end{align}
where $\boldsymbol{v}_{k}$ denotes the 3D velocity of the \ac{ue} at time $k$, which can be defined as
\begin{align}
    &\boldsymbol{v}_{k} = [v \cos \varpi_{k}, v\sin \varpi_{k}, 0]^{\mathsf{T}},\label{velocity}
\end{align}
since we assume the \ac{ue} only moves on the $x-y$ plane. Therefore, $d_{k}^{i}$ is positive, when the \ac{ue} approaches the $i$-th landmark, and it is negative, when the \ac{ue} moves away from the $i$-th landmark.

\vspace{-1mm}

\section{Proposed Doppler-assisted SLAM}
In this section,  the \ac{PCRB} is briefly introduced,  the \ac{PEB}, \ac{HEB}, \ac{CEB} and \ac{LEB} are derived, and the contribution of Doppler to the \ac{slam} system is analyzed.

\subsection{\ac{PCRB}} \label{sec:performance_bounds}
The \ac{PCRB} is the lower bound that is analogous to the \ac{CRB} but takes a Bayesian perspective and assumes that there is a prior on the parameters, indicating that the \ac{MSE} of an estimator should always be larger than the inverse of the \ac{PIM} \cite[Ch. 4.2]{van2004detection}. Suppose $\boldsymbol{\alpha}$ denotes a vector of the measurements, $\boldsymbol{\beta}$ denotes an $r$-dimensional estimated random parameter, and $\boldsymbol{g}(\boldsymbol{\alpha})$ denotes an estimate of $\boldsymbol{\beta}$, which is a function of $\boldsymbol{\alpha}$. The \ac{PCRB} on the estimation error has the form
\begin{equation}\label{eq:cramer_rao_bound}
\mathbb{E}[ \left(\boldsymbol{g}(\boldsymbol{\alpha}) - \boldsymbol{\beta}\right)\left(\boldsymbol{g}(\boldsymbol{\alpha}) - \boldsymbol{\beta}\right)^{\mathsf{T}} ] \geq \boldsymbol{J}^{-1},
\end{equation}
where $\boldsymbol{J}$ is the $r\times r$ \ac{PIM} with elements
\begin{equation}\label{eq:cramer_rao_bound_elements}
[\boldsymbol{J}]_{\mu,\nu}=\mathbb{E}\lbrack-\frac{\partial^{2}\log f(\boldsymbol{\alpha},\boldsymbol{\beta})}{\partial {\beta}_{\mu} \partial {\beta}_{\nu}}\rbrack \qquad \mu,\nu=1,\cdots,r,
\end{equation}
with $f(\boldsymbol{\alpha},\boldsymbol{\beta})$ denoting the joint density. The inequality in \eqref{eq:cramer_rao_bound} indicates that the \ac{MSE} of an estimator is larger than $\boldsymbol{J}^{-1}$ in the positive semidefinite sense. The \ac{PIM} is a counterpart to the \ac{FIM} for the \ac{PCRB}, and can be decomposed into two  parts  $ \boldsymbol{J} = \boldsymbol{J}_{\text{data}} + \boldsymbol{J}_{\text{prior}}$,
where $\boldsymbol{J}_{\text{data}}$  is the standard \ac{FIM}, which contains information obtained from the measurements, and $\boldsymbol{J}_{\text{prior}}$ is the priori information matrix, which contains prior information. The elements are given by (for $\mu,\nu=1,\cdots,r$) \vspace{-3mm}
\begin{align}
    \label{eq:cramer_rao_bound_elements_part1}
[\boldsymbol{J}_{\text{data}}]_{\mu,\nu}&=\mathbb{E}\lbrack-\frac{\partial^{2}\log f(\boldsymbol{\alpha}|\boldsymbol{\beta})}{\partial {\beta}_{\mu} \partial {\beta}_{\nu}}\rbrack, \\
\label{eq:cramer_rao_bound_elements_part2}
[\boldsymbol{J}_{\text{prior}}]_{\mu,\nu}&=\mathbb{E}\lbrack-\frac{\partial^{2}\log f(\boldsymbol{\beta})}{\partial {\beta}_{\mu} \partial {\beta}_{\nu}}\rbrack.
\end{align}

\subsection{Performance Bounds}

To compute the error bounds for the considered problem, we need to construct a  complete state of the system, denoted as $\breve{\boldsymbol{s}}_k = [(\boldsymbol{s}_k)^{\mathsf{T}}, \; (\boldsymbol{x}^{1}_{\text{LM},k})^{\mathsf{T}}, \; \ldots, \; (\boldsymbol{x}^{I}_{\text{LM},k})^{\mathsf{T}}]^{\mathsf{T}}$, the transition function of the complete state, denoted as $\breve{\boldsymbol{f}}(\breve{\boldsymbol{s}}_k) = [(\boldsymbol{f}(\boldsymbol{s}_{k})^{\mathsf{T}}, \; (\boldsymbol{x}^{1}_{\text{LM},k})^{\mathsf{T}}, \; \ldots, \; (\boldsymbol{x}^{I}_{\text{LM},k})^{\mathsf{T}}]^{\mathsf{T}}$, and the measurement function of the complete state given the ground-truth \ac{DA}, denoted as $\breve{\boldsymbol{h}}(\breve{\boldsymbol{s}}_k) = [(\bar{\boldsymbol{h}}(\boldsymbol{s}_{k},\boldsymbol{x}^{1}_{\text{LM},k})^{\mathsf{T}}, \; \ldots, \; (\bar{\boldsymbol{h}}(\boldsymbol{s}_{k},\boldsymbol{x}^{I}_{\text{LM},k})^{\mathsf{T}}]^{\mathsf{T}}$, with $\bar{\boldsymbol{h}}(\boldsymbol{s}_{k},\boldsymbol{x}^{i}_{\text{LM},k})=\boldsymbol{0}_{d_{\boldsymbol{z}} \times 1}$ when the landmark is not detected, and $\bar{\boldsymbol{h}}(\boldsymbol{s}_{k},\boldsymbol{x}^{i}_{\text{LM},k})$ is the same as $\boldsymbol{h}(\boldsymbol{s}_{k},\boldsymbol{x}^{i}_{\text{LM},k})$ when the landmark is detected. Then, if the ground-truth \ac{DA} is given, the \ac{PIM} can be recursively updated by \cite{tichavsky1998,EKPHD2021Ossi}
\begin{equation}\label{eq:fim_recursive}
\boldsymbol{J}_{k} =
   \underset{\boldsymbol{J}_{\text{data},k}}{\underbrace{\breve{\boldsymbol{H}}_{k}^{\T} \breve{\boldsymbol{R}}_k^{-1}\breve{\boldsymbol{H}}_{k}}} +
   \underset{\boldsymbol{J}_{\text{prior},k}}{\underbrace{\left(\breve{\boldsymbol{Q}}_k + \breve{\boldsymbol{F}}_{k} \boldsymbol{J}_{k-1}^{-1} \breve{\boldsymbol{F}}_{k}^{\T} \right)^{-1}}}.
\end{equation}
Here $\breve{\boldsymbol{H}}_{k}$ and $\breve{\boldsymbol{F}}_{k}$ are the Jacobian matrices of $\breve{\boldsymbol{h}}(\breve{\boldsymbol{s}}_k)$ and $\breve{\boldsymbol{f}}(\breve{\boldsymbol{s}}_k)$ with respect to $\breve{\boldsymbol{s}}_k$, evaluated at the true state, and $\breve{\boldsymbol{Q}}_k = \text{blkdiag}(\boldsymbol{Q},\boldsymbol{0}_{3I \times 3I})$ and $\breve{\boldsymbol{R}}_k = \text{blkdiag}(\boldsymbol{R}_k^1,\cdots,\boldsymbol{R}_k^I)$ denote the process and measurement noise covariances, respectively. Then, the \ac{PEB} can be computed as $\text{PEB}_k = \sqrt{[ \boldsymbol{J}_{k} ]_{1,1}^{-1} + [ \boldsymbol{J}_{k} ]_{2,2}^{-1}}$, the \ac{HEB} can be computed as $\text{HEB}_k = \sqrt{[ \boldsymbol{J}_{k} ]_{3,3}^{-1}}$, the \ac{CEB} can be computed as $\text{CEB}_k = \sqrt{[ \boldsymbol{J}_{k} ]_{4,4}^{-1}}$, and the \ac{LEB} of $i$-th landmark can be computed as $\text{LEB}_k^i = \sqrt{\sum_{\mu = 3i+2}^{3i + 4} [ \boldsymbol{J}_{k} ]_{\mu,\mu}^{-1}}$.

\subsection{Contribution of Doppler} \label{sec:contribution_Doppler}
From \eqref{eq:fim_recursive}, we can observe that the Doppler of the current measurements only contribute to $\boldsymbol{J}_{\text{data},k}$, as Doppler-related components are only involved in $\breve{\boldsymbol{H}}_{k}$ and $\breve{\boldsymbol{R}}_k$.  As $\breve{\boldsymbol{H}}_{k}$ is the Jacobian matrices of $\breve{\boldsymbol{h}}(\breve{\boldsymbol{s}}_k)$ with respect to $\breve{\boldsymbol{s}}_k$, it can be correspondingly decomposed into blocks
\begin{align}
    \breve{\boldsymbol{H}}_{k}  = \left[ \begin{array}{c : c c c c} \boldsymbol{A}_{k}^{1} & \boldsymbol{B}_{k}^{1}&\boldsymbol{0}_{6 \times 3}  &\cdots & \boldsymbol{0}_{6 \times 3}\\ \boldsymbol{A}_{k}^{2} & \boldsymbol{0}_{6 \times 3} & \boldsymbol{B}_{k}^{2} & \ddots & \boldsymbol{0}_{6 \times 3} \\ \vdots & \vdots& \ddots & \ddots & \vdots  \\ \boldsymbol{A}_{k}^{I} & \boldsymbol{0}_{6 \times 3}& \cdots& \boldsymbol{0}_{6 \times 3}& \boldsymbol{B}_{k}^{I}  \end{array}\right].
\end{align}
Here, $\boldsymbol{A}_{k}^{i}$ and $\boldsymbol{B}_{k}^{i}$ are the Jacobian matrices of $\bar{\boldsymbol{h}}(\boldsymbol{s}_{k},\boldsymbol{x}^{i}_{\text{LM}})$ with respect to $\boldsymbol{s}_{k}$ and $\boldsymbol{x}^{i}_{\text{LM}}$, respectively, which are $\boldsymbol{0}_{6\times4}$ and $\boldsymbol{0}_{6\times3}$ when the $i$-th landmark is not detected, and are
\begin{align}
    \boldsymbol{A}_{k}^{i} =\left[ \begin{array}{c}
    \tilde{\boldsymbol{A}}_{k}^{i}\\ \hdashline  (\tilde{\boldsymbol{a}}_{k}^{i})^{\mathsf{T}}
    \end{array}\right], \quad \boldsymbol{B}_{k}^{i} =\left[ \begin{array}{c}
    \tilde{\boldsymbol{B}}^{i}_{k}\\ \hdashline \vspace{-4mm} \\ (\tilde{\boldsymbol{b}}_{k}^{i})^{\mathsf{T}}
    \end{array}\right],
\end{align}
when the landmark is detected, where $ \tilde{\boldsymbol{A}}_{k}^{i}$ is the Jacobian of $[\tau_{k}^{i},(\boldsymbol{\theta}_{k}^{i})^{\mathsf{T}},(\boldsymbol{\phi}_{k}^{i})^{\mathsf{T}}]^{\mathsf{T}}$ with respect to $\boldsymbol{s}_{k}$, $ \tilde{\boldsymbol{B}}^{i}_{k}$ is the Jacobian of $[\tau_{k}^{i},(\boldsymbol{\theta}_{k}^{i})^{\mathsf{T}},(\boldsymbol{\phi}_{k}^{i})^{\mathsf{T}}]^{\mathsf{T}}$ with respect to $\boldsymbol{x}^{i}_{\text{LM}}$, and
\begin{align}
    (\tilde{\boldsymbol{a}}_{k}^{i})^{\mathsf{T}} & =[
    \begin{array}{cccc}
         \frac{\partial d_{k}^{i}}{\partial x_{k}} & \frac{\partial d_{k}^{i}}{\partial y_{k}} &\frac{\partial d_{k}^{i}}{\partial \varpi_{k}} & \frac{\partial d_{k}^{i}}{\partial b_{k}}
    \end{array}
    ],\\
    (\tilde{\boldsymbol{b}}_{k}^{i})^{\mathsf{T}} &  =[
    \begin{array}{cccc}
        \frac{\partial d_{k}^{i}}{\partial x_{\text{LM}}^{i}} & \frac{\partial d_{k}^{i}}{\partial y_{\text{LM}}^{i}} &\frac{\partial d_{k}^{i}}{\partial z_{\text{LM}}^{i}}
     \end{array}
    ].
\end{align}
\begin{comment}
\begin{align}
    &\boldsymbol{A}_{k}^{i} =\left[ \begin{array}{c}
    \tilde{\boldsymbol{A}}_{k}^{i}\\ \hdashline \vspace{-4mm} \\ (\tilde{\boldsymbol{a}}_{k}^{i})^{\mathsf{T}}
    \end{array}\right]= \left[ \begin{array}{c c c c} \frac{\partial \tau_{k}^{i}}{\partial x_{k}} & \frac{\partial \tau_{k}^{i}}{\partial y_{k}} &\frac{\partial \tau_{k}^{i}}{\partial \varpi_{k}} & \frac{\partial \tau_{k}^{i}}{\partial b_{k}}\\
    \frac{\partial \boldsymbol{\theta}_{k}^{i}}{\partial x_{k}} & \frac{\partial \boldsymbol{\theta}_{k}^{i}}{\partial y_{k}} &\frac{\partial \boldsymbol{\theta}_{k}^{i}}{\partial \varpi_{k}} & \frac{\partial \boldsymbol{\theta}_{k}^{i}}{\partial b_{k}}\\
    \frac{\partial \boldsymbol{\phi}_{k}^{i}}{\partial x_{k}} & \frac{\partial \boldsymbol{\phi}_{k}^{i}}{\partial y_{k}} &\frac{\partial \boldsymbol{\phi}_{k}^{i}}{\partial \varpi_{k}} & \frac{\partial \boldsymbol{\phi}_{k}^{i}}{\partial b_{k}} \vspace{1mm}\\ \hdashline \vspace{1.5mm}
    \frac{\partial d_{k}^{i}}{\partial x_{k}} & \frac{\partial d_{k}^{i}}{\partial y_{k}} &\frac{\partial d_{k}^{i}}{\partial \varpi_{k}} & \frac{\partial d_{k}^{i}}{\partial b_{k}}
    \end{array}\right], \\ \vspace{1mm}
    &\boldsymbol{B}_{k}^{i} =\left[ \begin{array}{c}
    \tilde{\boldsymbol{B}}^{i}_{k}\\ \hdashline \vspace{-4mm} \\ (\tilde{\boldsymbol{b}}_{k}^{i})^{\mathsf{T}}
    \end{array}\right] = \left[ \begin{array}{c c c} \frac{\partial \tau_{k}^{i}}{\partial x_{\text{LM}}^{i}} & \frac{\partial \tau_{k}^{i}}{\partial y_{\text{LM}}^{i}} &\frac{\partial \tau_{k}^{i}}{\partial z_{\text{LM}}^{i}}\\
    \frac{\partial \boldsymbol{\theta}_{k}^{i}}{\partial x_{\text{LM}}^{i}} & \frac{\partial \boldsymbol{\theta}_{k}^{i}}{\partial y_{\text{LM}}^{i}} &\frac{\partial \boldsymbol{\theta}_{k}^{i}}{\partial z_{\text{LM}}^{i}} \\
    \frac{\partial \boldsymbol{\phi}_{k}^{i}}{\partial x_{\text{LM}}^{i}} & \frac{\partial \boldsymbol{\phi}_{k}^{i}}{\partial y_{\text{LM}}^{i}} &\frac{\partial \boldsymbol{\phi}_{k}^{i}}{\partial z_{\text{LM}}^{i}} \vspace{1mm}\\ \hdashline \vspace{1.5mm}
    \frac{\partial d_{k}^{i}}{\partial x_{\text{LM}}^{i}} & \frac{\partial d_{k}^{i}}{\partial y_{\text{LM}}^{i}} &\frac{\partial d_{k}^{i}}{\partial z_{\text{LM}}^{i}}
    \end{array}\right],
\end{align}
when the $i$-th landmark is detected, with $\tilde{\boldsymbol{A}}_{k}$ and $\tilde{\boldsymbol{B}}_{k}$ denoting the parts for delay and angles, and $\tilde{\boldsymbol{a}}^{i}_{k}$ and $\tilde{\boldsymbol{b}}^{i}_{k}$ denoting the parts for Doppler.
\end{comment}
Furthermore, we assume $\boldsymbol{R}_{k}^{i}=\text{diag}(\tilde{\boldsymbol{R}}_{k}^{i},(\sigma_{\text{d},k}^{i})^{2})$.

  According to \eqref{eq:fim_recursive} and the inverse of a block matrix, $\boldsymbol{J}_{\text{data},k}$ is only non-zero in the diagonal block related to the UE state
\begin{align}

&[\boldsymbol{J}_{\text{data},k}]_{1:4,1:4}= \sum_{i=1}^{I}(\boldsymbol{A}_{k}^{i})^{\mathsf{T}}(\boldsymbol{R}_{k}^{i})^{-1}\boldsymbol{A}_{k}^{i}\nonumber\\&=
\sum_{i=1}^{I}\underbrace{(\tilde{\boldsymbol{A}}_{k}^{i})^{\mathsf{T}}(\tilde{\boldsymbol{R}}_{k}^{i})^{-1}\tilde{\boldsymbol{A}}_{k}^{i}}_{\text{non-Doppler related}}+\sum_{i=1}^{I}\underbrace{\tilde{\boldsymbol{a}}_{k}^{i}(\sigma_{\text{d},k}^{i})^{-2}(\tilde{\boldsymbol{a}}_{k}^{i})^{\mathsf{T}}}_{\text{Doppler related}},\label{eq:cramer_rao_bound_elements_data_decomposeA}
\end{align}
and the diagonal blocks related to each landmark
\begin{align}
&[\boldsymbol{J}_{\text{data},k}]_{3i+2:3i+4,3i+2:3i+4}= (\boldsymbol{B}_{k}^{i})^{\mathsf{T}}(\boldsymbol{R}_{k}^{i})^{-1}\boldsymbol{B}_{k}^{i}\nonumber\\&= \underbrace{
(\tilde{\boldsymbol{B}}_{k}^{i})^{\mathsf{T}}(\tilde{\boldsymbol{R}}_{k}^{i})^{-1}\tilde{\boldsymbol{B}}_{k}^{i}}_{\text{non-Doppler related}}+\underbrace{\tilde{\boldsymbol{b}}_{k}^{i}(\sigma_{\text{d},k}^{i})^{-2}(\tilde{\boldsymbol{b}}_{k}^{i})^{\mathsf{T}}}_{\text{Doppler related}}.\label{eq:cramer_rao_bound_elements_data_decomposeB}
\end{align}
Clearly, Doppler provides non-negative information to both UE and landmark states. Due to the block-diagonal structure of
$\boldsymbol{J}_{\text{data},k}$, \eqref{eq:cramer_rao_bound_elements_data_decomposeA} and \eqref{eq:cramer_rao_bound_elements_data_decomposeB} can be interpreted as the equivalent FIMs of the UE and landmark states.

To gain further insights, we expand
$\tilde{\boldsymbol{a}}^i_{k}$ and $\tilde{\boldsymbol{b}}^i_{k}$ as
\begin{align}
    & \tilde{\boldsymbol{a}}_{k}^{i}=\big[-\frac{\left([\boldsymbol{v}_{k}]_{1:2}-d_{k}^{i}[\boldsymbol{q}_{\text{AOA},k}^{i}]_{1:2}\right)^{\mathsf{T}}}{\Vert\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{UE},k}\Vert}, \boldsymbol{v}_{\perp,k}^{\textsf{T}}\boldsymbol{q}_{\text{AOA},k}^{i},0 \big]^{\mathsf{T}}, \\
    & \tilde{\boldsymbol{b}}_{k}^{i} =\frac{\partial \boldsymbol{x}_{\text{inc},k}^{i}}{\partial \boldsymbol{x}_{\text{LM}}^{i}}\frac{\boldsymbol{v}_{k}-d_{k}^{i}\boldsymbol{q}_{\text{AOA},k}^{i}}{\Vert\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{UE},k}\Vert},
\end{align}
where $\boldsymbol{v}_{\perp,k}=[ -v\sin \varpi_{k}, v\cos \varpi_{k}, 0]^{\mathsf{T}}$ (i.e., $\boldsymbol{v}^\top_{\perp,k}\boldsymbol{v}_{k}=0$). Then,  the additive Doppler-related parts in \eqref{eq:cramer_rao_bound_elements_data_decomposeA} and \eqref{eq:cramer_rao_bound_elements_data_decomposeB} can be summarized as follows:
\begin{itemize}
\item \emph{\ac{ue} position:} as $\boldsymbol{v}_{k}-d_{k}^{i}\boldsymbol{q}_{\text{AOA},k}^{i}=\boldsymbol{v}_{k}-\boldsymbol{v}_{k}^{\textsf{T}}\boldsymbol{q}_{\text{AOA},k}^{i}\boldsymbol{q}_{\text{AOA},k}^{i}$, the  information brought by Doppler is along the direction of the velocity, when the velocity is orthogonal to a path's $\boldsymbol{q}_{\text{AOA},k}^{i}$ direction, i.e., when $\boldsymbol{v}_{k}^{\textsf{T}}\boldsymbol{q}_{\text{AOA},k}^{i}=0$. On the other hand, when the velocity is parallel to a path's $\boldsymbol{q}_{\text{AOA},k}^{i}$ direction, the Doppler of that path does not provide direct Fisher information.  A similar argument holds for landmark locations.

\item \emph{Heading $\varpi_{k}$:} the direct contribution is from the projection of $\boldsymbol{v}_{\perp,k}$ on the direction of $\boldsymbol{q}_{\text{AOA},k}^{i}$ of each path. Hence, if the \ac{ue} velocity is parallel to a path's $\boldsymbol{q}_{\text{AOA},k}^{i}$ direction,
the direct Fisher information  is 0.

\item \emph{Clock bias $b_{k}$:} there is no direct Fisher information contribution. However, the direct contributions to the other dimensions still benefit estimate of the clock bias.
\end{itemize}
Therefore, if the \ac{ue} moves alongside the same or the opposite direction of the  $\boldsymbol{q}_{\text{AOA},k}^{i}$, the corresponding group of measurement  does not  has any Doppler-related information. Moreover, a smaller $(\sigma_{\text{d},k}^{i})^{2}$ results in larger additional Doppler-related information.

\begin{comment}

As $\boldsymbol{J}_{\text{data},k}$ is block-diagonal, the equivalent \ac{FIM} of $\boldsymbol{s}_{k}$ and the $i$-th landmark is the same as \eqref{eq:cramer_rao_bound_elements_data_decompose} and \eqref{eq:cramer_rao_bound_elements_data_decompose2}
\begin{align}
&\boldsymbol{J}_{\text{UE},k}^{\text{E}}=[\boldsymbol{J}_{\text{data},k}]_{1:4,1:4},\\
    &\boldsymbol{J}_{\text{LM},k}^{i,\text{E}}=[\boldsymbol{J}_{\text{data},k}]_{3i+2:3i+4,3i+2:3i+4} \qquad \quad  i=1,\cdots,I.
\end{align}
Clearly, the equivalent \acp{FIM} are non-decreasing, compared to the case without Doppler, as positive semi-definitive parts were added in \eqref{eq:cramer_rao_bound_elements_data_decompose} and \eqref{eq:cramer_rao_bound_elements_data_decompose2}. This is because the complete state remains the same, and an additional dimension of measurement is added. Moreover, a smaller $(\sigma_{\text{d},k}^{i})^{2}$ results in larger $\boldsymbol{J}_{\text{UE},k}^{\text{E}}$ and $\boldsymbol{J}_{\text{LM},k}^{i,\text{E}}$. As the \ac{PIM} is recursively updated by \eqref{eq:fim_recursive}, if $\boldsymbol{J}_{\text{data},k}$ non-decreases, $\boldsymbol{J}_{k}$ also non-decreases, so that the bounds non-increases.

\end{comment}

\section{EK-PMB(M) SLAM filter}
The \ac{slam} framework in this paper follows the EK-PMB \ac{slam} filter proposed in \cite{ge2022computationally}. In this section, the basics of the PMB(M) density and the EK-PMB(M) SLAM filter are summarized in the following.

\subsection{Basics of PMB(M) Density}
The map of the environment $\mathcal{X}$ is formulated as a \ac{PMBM} \ac{rfs}, which can be viewed as the union of two disjoint \acp{rfs}: $\mathcal{X}_{\mathrm{U}}$ for the set of undetected landmarks, which are the landmarks that have never been detected, and $\mathcal{X}_{\mathrm{D}}$ for the set of detected landmarks, which are the the landmarks that have been detected at least once before. Its density follows \cite{garcia2018poisson}
\begin{equation}
    f(\mathcal{X})=\sum_{\mathcal{X}_{\mathrm{U}}\biguplus\mathcal{X}_{\mathrm{D}}=\mathcal{X}}f_{\mathrm{P}}(\mathcal{X}_{\mathrm{U}})f_{\mathrm{MBM}}(\mathcal{X}_{\mathrm{D}}),\label{PMBM}
\end{equation}
where $\mathcal{X}_{\mathrm{U}}$ and  $\mathcal{X}_{\mathrm{D}}$ are modeled as a \ac{PPP} and a \ac{MBM}, respectively, with densities following
\begin{align}
    f_{\mathrm{P}}(\mathcal{X}_{\mathrm{U}})&=e^{-\int\lambda(\boldsymbol{x})\mathrm{d}\boldsymbol{x}}\prod_{\boldsymbol{x} \in \mathcal{X}_{\mathrm{U}} }\lambda(\boldsymbol{x}),\label{PPP} \\
     f_{\mathrm{MBM}}(\mathcal{X}_{\mathrm{D}})&\propto \sum_{j \in \mathbb{I}}w^{j}\sum_{\mathcal{X}^{1}\biguplus \dots \biguplus \mathcal{X}^{|\mathcal{X}_{\mathrm{D}}|}=\mathcal{X}_{\mathrm{D}}}\prod_{i=1}^{|\mathcal{X}_{\mathrm{D}}|}f^{j,i}_{\mathrm{B}}(\mathcal{X}^{i}).\label{MBM}
\end{align}
Here, $\lambda(\cdot)$ denotes the intensity function of the \ac{PPP} density, and  $\mathbb{I}$ is the index set of all global hypotheses with weights satisfying $\sum_{j\in\mathbb{I}}w^{j}=1,w^{j}\ge 0$, and the global hypotheses in \ac{slam} correspond to different \acp{DA} \cite{williams2015marginal}. Each individual component in \eqref{MBM} is termed as a Bernoulli process, and $f_{\mathrm{B}}^{j,i}(\cdot)$ denotes the Bernoulli density of the $i$-th landmark under the $j$-th global hypothesis, following
\begin{equation}
f^{j,i}_{\mathrm{B}}(\mathcal{X}^{i})=
\begin{cases}
1-r^{j,i} \quad& \mathcal{X}^{i}=\emptyset \\ r^{j,i}f^{j,i}(\boldsymbol{x}) \quad & \mathcal{X}^{i}=\{\boldsymbol{x}\} \\ 0 \quad & \mathrm{otherwise},
\end{cases}
\end{equation}
where $f^{j,i}(\cdot)$ denotes the state density, and $r^{j,i} \in [0,1]$ denotes the probability that the corresponding landmark exists.
We note that if there is only one mixture component in the \ac{MBM}, then \eqref{MBM} reduces to a \ac{MB}, and \eqref{PMBM} reduces to a \ac{PMB}.

\subsection{ EK-PMB(M) SLAM Filter Recursion} \label{EK-PMB}
The EK-PMB(M) SLAM filter follows the Bayesian filtering recursion with \acp{rfs}, and it computes the joint posterior $f(\boldsymbol{s}_{k+1},\mathcal{X}|\mathcal{Z}_{1:k+1})$  by \cite{ge2022computationally}
\begin{align}
     f(\boldsymbol{s}_{k+1},\mathcal{X}|\mathcal{Z}_{1:k+1}&) \propto   \ell(\mathcal{Z}_{k+1}|\boldsymbol{s}_{k+1},\mathcal{X}) f(\mathcal{X}|\mathcal{Z}_{1:k}) \nonumber \\ &\times f(\boldsymbol{s}_{k+1}|\mathcal{Z}_{1:k}), \label{joint_posterior}
\end{align}
with $f(\boldsymbol{s}_{k+1}|\mathcal{Z}_{1:k})=\int f(\boldsymbol{s}_{k}|\mathcal{Z}_{1:k})f(\boldsymbol{s}_{k+1}|\boldsymbol{s}_{k}) \text{d} \boldsymbol{s}_{k}$ and $\ell(\mathcal{Z}_{k+1}|\boldsymbol{s}_{k+1},\mathcal{X})$ denoting the RFS likelihood function, given by  \cite[eqs.\,(5)--(6)]{garcia2018poisson}. Instead of tracking the joint density, the filter keeps track of the marginal  \ac{ue} $f(\boldsymbol{s}_{k}|\mathcal{Z}_{1:k})$ and map $f(\mathcal{X}|\mathcal{Z}_{1:k})$ posteriors by marginalizing out the map and the \ac{ue}  state from the joint posterior, respectively. The marginal posteriors are given by
\begin{align}
     f(\boldsymbol{s}_{k+1}|\mathcal{Z}_{1:k+1})  &= \int  f(\boldsymbol{s}_{k+1},\mathcal{X}|\mathcal{Z}_{1:k+1})\delta \mathcal{X},  \label{eq:UpMarObjMarg1}\\
    f(\mathcal{X}|\mathcal{Z}_{1:k+1}) &= \int f(\boldsymbol{s}_{k+1},\mathcal{X}|\mathcal{Z}_{1:k+1}) \mathrm{d}\boldsymbol{s}_{k+1},\label{eq:UpMarObjMarg2}
\end{align}
where $\int \psi(\mathcal{X})\delta \mathcal{X}$ refers to the set integral \cite[eq.~(4)]{williams2015marginal}.

The EK-PMB(M) \ac{slam} filter proposed in \cite{ge2022computationally}  determines $\gamma\ge 1$ most likely \acp{DA} for each prior global hypothesis with corresponding weights, and the joint posterior of the \ac{ue} state and the map conditioned on each \ac{DA} is computed by utilizing an \ac{ekf}-like update step. Then the joint posterior in \eqref{joint_posterior} can be acquired by the weighted summation of the densities for \acp{DA}, followed by \eqref{eq:UpMarObjMarg1} and \eqref{eq:UpMarObjMarg2} to compute the \ac{ue} state and the map. To avoid the exponential increase of \acp{DA}, we can approximate the PMBM density to a PMB density at the end of each update step, which reduces the complexity significantly.

It is important to mention that the weights of $\gamma\ge 1$ most likely \acp{DA} for each prior global hypothesis are computed by using measurements \cite[eq.\,(29)]{ge2022computationally}. Since more information is provided with the involvement of Doppler, finding the correct \acp{DA} becomes more likely,
as correct local hypothesis weights computed in \cite[eqs.\,(22)--(24)]{ge2022computationally} become more prominent.

\section{Results}

\subsection{Simulation Environment}
We consider a 5G downlink scenario with a single \ac{bs} located at $[0 \, \text{m},0 \, \text{m},40\, \text{m}]^{\mathsf{T}} $, 4 \acp{va} located at  $[\pm 200 \, \text{m},0 \, \text{m},40 \, \text{m}]^{\mathsf{T}}$, $[0 \, \text{m}, \pm 200 \, \text{m},40 \, \text{m}]^{\mathsf{T}}$, and 4 \acp{sp} located at  $[\pm 99 \, \text{m},0 \, \text{m},10 \, \text{m}]^{\mathsf{T}}$, $[0 \, \text{m}, \pm 99 \, \text{m},10 \, \text{m}]^{\mathsf{T}}$. The \ac{ue} does a counterclockwise  constant turn-rate movement around the \ac{bs} according to \eqref{eq:dynamic_model}, with $v = 22.22 \text{ m/s}$,  $\omega = \pi/10 \text{ rad/s}$, $T=0.5 \text{ s}$, and $\mathbf{Q}= \text{diag}(0.2^2 \text{ m}^2, \,0.2^2 \text{ m}^2, \,0.001^2 \text{ rad}^2, \,0.2^2 \text{ m}^2)$, and it takes $K = 40$ samples for the \ac{ue} to circle the road once. The measurement covariance matrix is set as  $\boldsymbol{R}_{k}^{i} = \text{blkdiag}(10^{-2} \text{ m}^2, 2.5\times 10^{-3} \cdot \mathbf{I}_4  \text{ rad}^2, \sigma_{\textsf{d}}^{2})$. We initialize the \ac{ue} at $[70.7285 \, \text{m},0 \, \text{m},\pi/2 \, \text{rad}, 300 \, \text{m}]^{\mathsf{T}}$ with covariance as $\text{diag}[ 0.3 \, \text{m}^{2},0.3 \, \text{m}^{2} \, 0.0052 \, \text{rad}^{2},0.3 \, \text{m}^{2}]$. The \ac{bs} is a priori known to the \ac{ue}. We implemented the EK-PMB SLAM filter with considering $\gamma=10$ best \acp{DA} every time step, and compared the results of cases using different levels of $\sigma_{\textsf{d}}$ and without considering Doppler. The mapping  performance is quantified by the \ac{gospa} distance \cite{rahmathullah2017generalized} for both \acp{va} and \acp{sp}, separately, the positioning performance is evaluated by the root mean squared error (RMSE), and we also compare bounds of different cases. More details on parametric settings of the filter can be found in \cite{ge2022computationally}.  The results are averaged over 100 Monte Carlo  simulations.

\subsection{Results and Discussion}
We firstly analyze how different levels of $\sigma_{\textsf{d}}$ affect the bounds. Fig.\,\ref{Fig.bounds} shows how the bounds
change with $\sigma_{\textsf{d}}$, compared with the bounds without considering Doppler. We observe that involving Doppler as a dimension of the measurement improves the positioning and mapping performances of the \ac{slam} system, as the error bounds are lower than in the case where the Doppler measurement is ignored. This effect is stronger when Doppler measurements are more accurate, i.e., when $\sigma_{\textsf{d}}$ is smaller.

\begin{figure}[t]
\center
\input{bounds}\vspace{-5mm}
\caption{The \ac{PEB},  the \ac{HEB}, the \ac{CEB} and the average \ac{LEB} change with $\sigma_{\textsf{d}}$. The benchmarks are the bounds of the case without considering Doppler.}\vspace{-4mm}
\label{Fig.bounds}
\end{figure}

\begin{figure}[t]
\center
\input{mapping}\vspace{-5mm}
\caption{Comparison of mapping performances for VAs between two cases: with and without considering Doppler.}\vspace{-5mm}
\label{Fig.mapping}
\end{figure}

To validate the benefits of Doppler in mmWave radio SLAM,  we then implement the EK-PMB SLAM filter on two cases: 1) involving Doppler into the measurement with $\sigma_{\textsf{d}}=0.1 \text{ m/s}$ and 2) without considering Doppler, and compare the mapping performance for both \acp{va} and \acp{sp}, and the positioning performance between two cases, as illustrated in Fig.\,\ref{Fig.mapping}, Fig.\,\ref{Fig.mapping_SP} and Fig.\,\ref{Fig.bar}, respectively. Fig.\,\ref{Fig.mapping} and Fig.\,\ref{Fig.mapping_SP} demonstrate that the SLAM filter can map the environment for both cases, and the mapping accuracy improves with more measurements being received, as overall, all GOSPA distances gradually decrease over time.  Clearly, when Doppler is involved, the SLAM filter has better mapping performance, as the solid lines are lower than the dashed lines in both Fig.\,\ref{Fig.mapping} and Fig.\,\ref{Fig.mapping_SP}. Fig.\,\ref{Fig.bar} indicates that considering Doppler results in better \ac{ue} state estimates, as lower RMSEs can be acquired. The reasons are that considering Doppler as a dimension of the measurement brings lower bounds. In addition, it also helps the SLAM filter to solve the \ac{DA} problem, which improves the average weight for the correct \ac{DA} from 0.6763 to 0.8622, with weights of all selected $\gamma=10$ \acp{DA} summed to 1 every time step.

\begin{figure}
\center
\input{mapping_SP}
\caption{Comparison of mapping performances for SPs between two cases: with and without considering Doppler.}
\label{Fig.mapping_SP}
\end{figure}

\begin{figure}
\center
\input{Doppler for 5G SLAM/positioning}
\caption{Comparison of \ac{ue} state estimation between two cases: with and without considering Doppler.}
\label{Fig.bar}
\end{figure}

\section{Conclusions}
In this paper, we exploited Doppler as a part of measurement in bistatic mmWave radio SLAM, formulated the EK-PMB SLAM filter on the new measurement model, and provided the \ac{PCRB} for the model. Our results theoretically indicate that the involvement of Doppler helps the SLAM filter to acquire better mapping and positioning performance than the case without considering Doppler. The better Doppler observations are, the lower \acp{PCRB} and the more accurate  map and \ac{ue} state estimates will be. The implementation of the EK-PMB SLAM filter validates the theoretical benefits of involving  Doppler. Future work will include extending the \ac{ue} model to a model with unknown speed, the inclusion of high-dimensional channel estimation, as well as extending the \ac{slam} problem to a \ac{slat} problem.

\scriptsize{
\section*{Acknowledgment}
This work was partially supported by the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by Knut and Alice Wallenberg Foundation, and the Vinnova 5GPOS project under grant 2019-03085, by the Swedish Research Council under grant 2018-03705, and by the Academy of Finland under the grants \#315858, \#328214, \#319994, \#323244, \#346622.
}
\balance

\section{Introduction}
MmWave communications in 5G and beyond are useful for \acf{slam} applications, due to geometrical propagation channels, large bandwidths, and large antenna arrays \cite{nurmi2017multi}. Signals sent from the \acf{bs} reach the \acf{ue} via the propagation channel, which is determined by the geometric relationships among the propagation environment, the \ac{ue}, and the \ac{bs}. Large bandwidths and antenna arrays result in high temporal and  spatial resolutions \cite{patwari2005,larsen2009}. Therefore, state-of-the-art channel estimators can provide accurate estimates for multipath components by using the received signals, in terms of groups of channel gain, \ac{TOA}, \ac{AOA}, \ac{AOD}, and Doppler, which contain  the necessary information for \ac{slam} \cite{wymeersch20175g}. Although Doppler contains geometric information,
it is usually ignored.

The related works can be divided into two areas: works that exploit Doppler for radio positioning or mapping and works in the area of radio SLAM. Doppler has been used in most radars for mapping and tracking \cite{patole2017automotive}, but limited works have been done in radio scenarios. Doppler is used to localize radio emitters in \cite{amar2008localization}, but the proposed method can only be used for narrow-band signals, while \cite{Han2016} shows that the Doppler shift can provide more direction information for localization, and \cite{kakkavas2019performance} shows that the mobility can significantly improve the \ac{nlos}-only scenario in MIMO mmWave system,  indicating the involvement of the Doppler shift brings gain in the localization accuracy. Doppler is used for tracking \acp{ue}' positions and velocities in a Wi-Fi-based system in  \cite{qian2017widar}. However, these methods do not solve the \ac{slam} problem.
Several approaches have been proposed to address the mmWave radio \ac{slam} problem, including geometry-based methods \cite{wen20195g, yassin2018mosaic}, message passing-based methods \cite{Erik_BPSLAM_TWC2019,Rico_BPSLAM_JSTSP2019}, and \ac{rfs}-based methods \cite{kim20205g,ge20205GSLAM,ge2022computationally}. \ac{rfs}-based methods
can  handle uncertainties, as well as inherently deal with challenges of the unknown number of landmarks, unknown \acp{DA}, misdetections, and clutter measurements in radio \ac{slam}. Within these \ac{rfs}-based methods, the probability hypothesis density (PHD) filter is used in \cite{kim20205g}, which does not have an explicit enumeration of \acp{DA}, and the \ac{PMBM} filter is used in \cite{ge20205GSLAM,ge2022computationally}, which enumerates all possible DAs explicitly, thus allowing for improved performance. However, the Doppler is not considered in these works. To our best knowledge, the inclusion of Doppler in radio SLAM has not yet been conducted in the existing literature.

In this paper, we harness the Doppler component in bistatic mmWave radio \ac{slam} and analyze how Doppler benefits the \ac{slam} filter. The main contributions of this paper are summarized as follows: (\emph{i}) we derive the \ac{PEB}, \ac{HEB}, \ac{CEB} of the \ac{ue} and the \acp{LEB} of landmarks by computing the \acf{PCRB} for the system with utilizing Doppler;  (\emph{ii}) we analyze the effect of the quality of the Doppler measurement on the bounds of the \ac{ue} and the average \ac{LEB} of landmarks; (\emph{iii}) we extend our previous work in \cite{ge2022computationally} by involving the Doppler as a dimension of the measurement in the  \ac{ek}-\ac{PMB} \ac{slam} filter, and validate the benefits of exploiting Doppler in the system through numerical experiments in the mmWave network context.

\subsubsection*{Notations}
Scalars (e.g., $x$) are denoted in italic, vectors (e.g., $\boldsymbol{x}$) in bold lower-case letters  with $\left\Vert\boldsymbol{x}\right\Vert$ representing its L2-norm, matrices (e.g., $\boldsymbol{X}$) in bold capital letters, sets  (e.g., $\mathcal{X}$) in calligraphic with $\left|\mathcal{X}\right|$ representing its cardinality. The transpose is denoted by $(\cdot)^{\mathsf{T}}$, the Hermitian transpose is denoted by $(\cdot)^{\mathsf{H}}$, the union of mutually disjoint sets is denoted by $\uplus$, the Kronecker product is denoted by $\otimes$, the expectation is denoted by $\mathbb{E}[\cdot]$, $\mathcal{N}(\boldsymbol{u},\boldsymbol{\Sigma})$ denotes a multivariate Gaussian distribution with mean $\boldsymbol{u}$ and covariance $\boldsymbol{\Sigma}$, and $d_{\boldsymbol{x}}=\text{dim}(\boldsymbol{x})$.

\vspace{-1mm}
\section{System model}
In this section, the \ac{ue} model, the environment model, the signal model, and the measurement model for a mmWave radio downlink scenario as shown in Fig.\,\ref{fig:scenario} are briefly introduced.

\begin{figure}
\centering
    \includegraphics[width=0.98\linewidth]{scenario_details_1.pdf}
\caption{A mmWave downlink scenario with the environment of a \ac{bs}, a \ac{ue}, a reflecting surface, and a small object. The reflecting surface can be modeled as a \acf{va}, which is the reflection of the \ac{bs} with respect to the surface. The small object is modeled as a \acf{sp}. The \ac{bs} sends signals to the \ac{ue} via \acf{los} path and/or \ac{nlos} paths, shown as the blue line and  purple lines, respectively. The channel parameters of each path depend on the underlying geometry.}
\label{fig:scenario}
\end{figure}

\subsection{State Models}
In this paper, a single-user scenario is considered, where the \ac{ue} does a constant turn-rate movement around a known \ac{bs} on the $x-y$ plane.
We denote the dynamic state of the \ac{ue} at time step $k$ as $\boldsymbol{s}_{k}=[x_{k},y_{k},\varpi_{k},b_{k}]^{\mathsf{T}}$, containing the \ac{ue} position on $x$- and $y$-axis, heading, and clock bias, respectively. The \ac{ue} dynamics can be expressed as \cite{roth2014}
\begin{equation}\label{eq:dynamic_model}
     \underset{\boldsymbol{s}_{k}}{\underbrace{\begin{bmatrix} x_{k} \\ y_{k} \\ \varpi_{k} \\ b_{k} \end{bmatrix}}} = \underset{\boldsymbol{f}\left(\boldsymbol{s}_{k-1} \right)}{\underbrace{\begin{bmatrix} x_{k-1} + \tfrac{2v}{\omega}  \sin \! \left(\tfrac{\omega T}{2}\right)  \cos \! \left(\varpi_{k-1} + \tfrac{\omega T}{2}\right) \\ y_{k-1} + \tfrac{2v}{\omega}\sin \! \left(\tfrac{\omega T}{2}\right) \sin \! \left(\varpi_{k-1} + \tfrac{\omega T}{2}\right) \\ \varpi_{k-1} + \omega T \\ b_{k-1} \end{bmatrix}}} + \boldsymbol{\eta}_{k},
\end{equation}
where $v$, $\omega$, and $T$ are known control inputs, denoting the speed,  the turn-rate,  and the sampling interval, respectively, and  $\boldsymbol{\eta}_k \sim \mathcal{N}\left(\boldsymbol{0},\boldsymbol{Q}\right)$.

The considered scenario contains a single known \ac{bs}, and multiple unknown reflecting surfaces and small objects, which are modeled as \acp{va} and \acp{sp}, respectively. The \ac{bs} sends downlink signals to the \ac{ue}, which can reach the \ac{ue} directly, termed as the \ac{los} path, and/or can be reflected by the reflecting surfaces or scattered by the small objects and reach the  \ac{ue}, termed as \ac{nlos} paths  (see Fig.\,\ref{fig:scenario}). In this paper, we assume that there is at most one path associated to each landmark every time step. The landmark state  can be modeled as $\boldsymbol{x} = [\boldsymbol{x}^{\textsf{T}}_{\text{LM}},m]^{\textsf{T}}$,  with $\boldsymbol{x}_{\text{LM}} \in \mathbb{R}^{3}$ denoting the landmark location, and $m\in\{\text{BS},\text{VA},\text{SP}\}$ denoting the landmark type. Therefore, the map of the environment can be modeled by a set of landmark $\mathcal{X}=\{\boldsymbol{x}^{1},\dots, \boldsymbol{x}^{I}\}$, with ${I}$ representing the total number of landmarks.

\subsection{Signal Model}
\label{Signal model}
The \ac{bs} sends downlink signals to the \ac{ue} at every time step with a period of $T$ seconds. These signals can reach the \ac{ue} via \ac{los} path and/or \ac{nlos} paths, and the received signal for OFDM symbol $n$, at subcarrier $\kappa$ and time step $k$ can be expressed as~\cite{heath2016overview}
\begin{align}
    \boldsymbol{y}_{\kappa,n,k} =\boldsymbol{W}_{n,k}^{\mathsf{H}}&\sum _{i=1}^{I_{k}}g_{k}^{i}\boldsymbol{a}_{\text{R}}(\boldsymbol{\theta}_{k}^{i})\boldsymbol{a}_{\text{T}}^{\mathsf{H}}(\boldsymbol{\phi}_{k}^{i}) e^{-\jmath 2\pi \kappa \Delta f \tau_{k}^{i}}\nonumber\\&\times e^{\jmath 2\pi n T_{\text{mea}}  d_{k}^{i}/\varsigma} \boldsymbol{p}_{n,k} + \boldsymbol{W}_{n,k}^{\mathsf{H}}\boldsymbol{n}_{\kappa,n,k}, \label{eq:FreqObservation}
\end{align}
where $\boldsymbol{y}_{\kappa,n,k}$ denotes the received signal, $\boldsymbol{p}_{n,k}$ denotes the pre-coded pilot signal,  $\boldsymbol{n}_{\kappa,n,k}$ denotes a white Gaussian noise, $\boldsymbol{W}_{n,k}$ denotes a combining matrix, $I_{k}$ denotes the number of all visible landmarks at time $k$, $\boldsymbol{a}_{\text{R}}(\cdot)$ and $\boldsymbol{a}_{\text{T}}(\cdot)$ denote the steering vectors of the receiver and transmitter antenna arrays, respectively, $\Delta f$ denotes the subcarrier spacing, $T_{\text{mea}}$ denotes transmission duration of each symbol, and $\varsigma$ denotes the wavelength.
Each path $i$ can be described by a complex gain $g_{k}^{i}$, a \ac{TOA} $\tau_{k}^{i}$, an \ac{AOA} pair $\boldsymbol{\theta}_{k}^{i}$ in azimuth and elevation,  an \ac{AOD} pair $\boldsymbol{\phi}_{k}^{i}$ in azimuth and elevation, and a Doppler $d_{k}^{i}$ which we express in [m/s]. These parameters can be  estimated by a parametric
channel estimation algorithm from $\boldsymbol{y}_{\kappa,n,k}$, for example, by \cite{richter2005estimation,alkhateeb2014channel,jiang2021high}, which is out of the scope of this paper, and the estimation results are utilized directly. Adding Doppler can make paths more resolvable, as there may be some paths not resolvable in delay or angle domains, but resolvable in the Doppler domain.

\subsection{Measurement Model}
At time step $k$, the channel estimator provides a set of measurements $\mathcal{Z}_{k}=\{\boldsymbol{z}_{k}^{1},\dots, \boldsymbol{z}_{k}^{\hat{{I}}_{k}} \}$, where usually $\hat{{I}}_{k} \neq {{I}}_{k}$, as some visible landmarks may be misdetected and there could be some clutter measurements. Please note that the source of each element in $\mathcal{Z}_{k}$ is unknown.
If the measurement noise is zero-mean Gaussian, the measurement originating from landmark $\boldsymbol{x}^{i}$ at time step $k$ can be described as follows
\begin{align}
    \boldsymbol{z}^{i}_{k}=\boldsymbol{h}(\boldsymbol{x}^{i},\boldsymbol{s}_{k})+\boldsymbol{\epsilon}_k^i,\label{pos_to_channelestimation}
\end{align}

where
$\boldsymbol{h}(\boldsymbol{x}^{i},\boldsymbol{s}_{k})=[\tau_{k}^{i},(\boldsymbol{\theta}_{k}^{i})^{\mathsf{T}},(\boldsymbol{\phi}_{k}^{i})^{\mathsf{T}}, d_{k}^{i}]^{\mathsf{T}}$ represents the nonlinear function that transforms the geometric information to the channel parameters, and $\boldsymbol{\epsilon}_k \sim \mathcal{N}\left(\boldsymbol{0},\boldsymbol{R}_{k}^{i}\right)$

These channel parameters depend on the geometric relationships among the \ac{bs}, the \ac{ue} and the landmarks. Specifically, \ac{TOA} $\tau_{k}^{i}$ can be defined as
\begin{align}
    \tau_{k}^{i} = \begin{cases} \Vert\boldsymbol{x}_{\text{BS}}-\boldsymbol{x}_{\text{UE},k}\Vert/c + b_{k} & m^{i}=\text{BS} \\  (\Vert\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{UE},k}\Vert + \Vert\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{BS}}\Vert)/c +b_{k} & m^{i}\neq\text{BS},\end{cases} \label{delay}
\end{align}
where $\boldsymbol{x}_{\text{UE},k}= [x_{k},y_{k},0]^{\textsf{T}}$ denotes the 3D position of the \ac{ue}, $\boldsymbol{x}_{\text{inc},k}^{i}$ is the incidence point of the $i$-th path on the corresponding landmark at time step $k$, which can be determined by $\boldsymbol{x}_{\text{LM}}^{i}$ and $\boldsymbol{x}_{\text{UE},k}$, and $c$ is the speed of light. As $\boldsymbol{\theta}_{k}^{i}$ is determined by the arrival direction of the signal, which can be calculated by
\begin{align}
    \boldsymbol{q}_{\text{AOA},k}^{i} = \begin{cases} (\boldsymbol{x}_{\text{BS}}-\boldsymbol{x}_{\text{UE},k})/\Vert\boldsymbol{x}_{\text{BS}}-\boldsymbol{x}_{\text{UE},k}\Vert & m^i=\text{BS} \\ (\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{UE},k})/\Vert\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{UE},k}\Vert & m^i\neq\text{BS},\end{cases}
\end{align}
we can define $\boldsymbol{\theta}_{k}^{i}$ as
\begin{align}
    &\theta_{\text{az},k}^{i} = \mathrm{arctan2}([\boldsymbol{q}_{\text{AOA},k}^{i}]_{2},[\boldsymbol{q}_{\text{AOA},k}^{i}]_{1})-\varpi_{k},\label{AOAAZ}\\
    &\theta_{\text{el},k}^{i} = \arcsin([\boldsymbol{q}_{\text{AOA},k}^{i}]_{3},\Vert\boldsymbol{q}_{\text{AOA},k}^{i}\Vert).\label{AOAEL}
\end{align}
Similarly, the departure direction can be calculated by  \begin{align}
    \boldsymbol{q}_{\text{AOD},k}^{i} = \begin{cases} (\boldsymbol{x}_{\text{UE},k}-\boldsymbol{x}_{\text{BS}})/\Vert(\boldsymbol{x}_{\text{UE},k}-\boldsymbol{x}_{\text{BS}})\Vert & m^i=\text{BS} \\ (\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{BS}})/\Vert\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{BS}}\Vert & m^i\neq\text{BS},\end{cases}
\end{align}
and $\boldsymbol{\phi}_{k}^{i}$ can be defined as
\begin{align}
    &\phi_{\text{az},k}^{i} = \mathrm{arctan2}([\boldsymbol{q}_{\text{AOD},k}^{i}]_{2},[\boldsymbol{q}_{\text{AOD},k}^{i}]_{1}),\label{AODAZ}\\
    &\phi_{\text{el},k}^{i} = \arcsin([\boldsymbol{q}_{\text{AOD},k}^{i}]_{3},\Vert\boldsymbol{q}_{\text{AOD},k}^{i}\Vert).\label{AODEL}
\end{align}
The Doppler $d_{k}^{i}$ can be calculated by projecting the \ac{ue} velocity on the direction of $\boldsymbol{q}_{\text{AOA},k}^{i}$, as displayed in Fig.\,\ref{fig:scenario}.  It can be computed by
\begin{align}
    &d_{k}^{i} = \boldsymbol{v}_{k}^{\textsf{T}}\boldsymbol{q}_{\text{AOA},k}^{i},\label{Doppler}
\end{align}
where $\boldsymbol{v}_{k}$ denotes the 3D velocity of the \ac{ue} at time $k$, which can be defined as
\begin{align}
    &\boldsymbol{v}_{k} = [v \cos \varpi_{k}, v\sin \varpi_{k}, 0]^{\mathsf{T}},\label{velocity}
\end{align}
since we assume the \ac{ue} only moves on the $x-y$ plane. Therefore, $d_{k}^{i}$ is positive, when the \ac{ue} approaches the $i$-th landmark, and it is negative, when the \ac{ue} moves away from the $i$-th landmark.

\vspace{-1mm}

\section{Proposed Doppler-assisted SLAM}
In this section,  the \ac{PCRB} is briefly introduced,  the \ac{PEB}, \ac{HEB}, \ac{CEB} and \ac{LEB} are derived, and the contribution of Doppler to the \ac{slam} system is analyzed.

\subsection{\ac{PCRB}} \label{sec:performance_bounds}
The \ac{PCRB} is the lower bound that is analogous to the \ac{CRB} but takes a Bayesian perspective and assumes that there is a prior on the parameters, indicating that the \ac{MSE} of an estimator should always be larger than the inverse of the \ac{PIM} \cite[Ch. 4.2]{van2004detection}. Suppose $\boldsymbol{\alpha}$ denotes a vector of the measurements, $\boldsymbol{\beta}$ denotes an $r$-dimensional estimated random parameter, and $\boldsymbol{g}(\boldsymbol{\alpha})$ denotes an estimate of $\boldsymbol{\beta}$, which is a function of $\boldsymbol{\alpha}$. The \ac{PCRB} on the estimation error has the form
\begin{equation}\label{eq:cramer_rao_bound}
\mathbb{E}[ \left(\boldsymbol{g}(\boldsymbol{\alpha}) - \boldsymbol{\beta}\right)\left(\boldsymbol{g}(\boldsymbol{\alpha}) - \boldsymbol{\beta}\right)^{\mathsf{T}} ] \geq \boldsymbol{J}^{-1},
\end{equation}
where $\boldsymbol{J}$ is the $r\times r$ \ac{PIM} with elements
\begin{equation}\label{eq:cramer_rao_bound_elements}
[\boldsymbol{J}]_{\mu,\nu}=\mathbb{E}\lbrack-\frac{\partial^{2}\log f(\boldsymbol{\alpha},\boldsymbol{\beta})}{\partial {\beta}_{\mu} \partial {\beta}_{\nu}}\rbrack \qquad \mu,\nu=1,\cdots,r,
\end{equation}
with $f(\boldsymbol{\alpha},\boldsymbol{\beta})$ denoting the joint density. The inequality in \eqref{eq:cramer_rao_bound} indicates that the \ac{MSE} of an estimator is larger than $\boldsymbol{J}^{-1}$ in the positive semidefinite sense. The \ac{PIM} is a counterpart to the \ac{FIM} for the \ac{PCRB}, and can be decomposed into two  parts  $ \boldsymbol{J} = \boldsymbol{J}_{\text{data}} + \boldsymbol{J}_{\text{prior}}$,
where $\boldsymbol{J}_{\text{data}}$  is the standard \ac{FIM}, which contains information obtained from the measurements, and $\boldsymbol{J}_{\text{prior}}$ is the priori information matrix, which contains prior information. The elements are given by (for $\mu,\nu=1,\cdots,r$) \vspace{-3mm}
\begin{align}
    \label{eq:cramer_rao_bound_elements_part1}
[\boldsymbol{J}_{\text{data}}]_{\mu,\nu}&=\mathbb{E}\lbrack-\frac{\partial^{2}\log f(\boldsymbol{\alpha}|\boldsymbol{\beta})}{\partial {\beta}_{\mu} \partial {\beta}_{\nu}}\rbrack, \\
\label{eq:cramer_rao_bound_elements_part2}
[\boldsymbol{J}_{\text{prior}}]_{\mu,\nu}&=\mathbb{E}\lbrack-\frac{\partial^{2}\log f(\boldsymbol{\beta})}{\partial {\beta}_{\mu} \partial {\beta}_{\nu}}\rbrack.
\end{align}

\subsection{Performance Bounds}

To compute the error bounds for the considered problem, we need to construct a  complete state of the system, denoted as $\breve{\boldsymbol{s}}_k = [(\boldsymbol{s}_k)^{\mathsf{T}}, \; (\boldsymbol{x}^{1}_{\text{LM},k})^{\mathsf{T}}, \; \ldots, \; (\boldsymbol{x}^{I}_{\text{LM},k})^{\mathsf{T}}]^{\mathsf{T}}$, the transition function of the complete state, denoted as $\breve{\boldsymbol{f}}(\breve{\boldsymbol{s}}_k) = [(\boldsymbol{f}(\boldsymbol{s}_{k})^{\mathsf{T}}, \; (\boldsymbol{x}^{1}_{\text{LM},k})^{\mathsf{T}}, \; \ldots, \; (\boldsymbol{x}^{I}_{\text{LM},k})^{\mathsf{T}}]^{\mathsf{T}}$, and the measurement function of the complete state given the ground-truth \ac{DA}, denoted as $\breve{\boldsymbol{h}}(\breve{\boldsymbol{s}}_k) = [(\bar{\boldsymbol{h}}(\boldsymbol{s}_{k},\boldsymbol{x}^{1}_{\text{LM},k})^{\mathsf{T}}, \; \ldots, \; (\bar{\boldsymbol{h}}(\boldsymbol{s}_{k},\boldsymbol{x}^{I}_{\text{LM},k})^{\mathsf{T}}]^{\mathsf{T}}$, with $\bar{\boldsymbol{h}}(\boldsymbol{s}_{k},\boldsymbol{x}^{i}_{\text{LM},k})=\boldsymbol{0}_{d_{\boldsymbol{z}} \times 1}$ when the landmark is not detected, and $\bar{\boldsymbol{h}}(\boldsymbol{s}_{k},\boldsymbol{x}^{i}_{\text{LM},k})$ is the same as $\boldsymbol{h}(\boldsymbol{s}_{k},\boldsymbol{x}^{i}_{\text{LM},k})$ when the landmark is detected. Then, if the ground-truth \ac{DA} is given, the \ac{PIM} can be recursively updated by \cite{tichavsky1998,EKPHD2021Ossi}
\begin{equation}\label{eq:fim_recursive}
\boldsymbol{J}_{k} =
   \underset{\boldsymbol{J}_{\text{data},k}}{\underbrace{\breve{\boldsymbol{H}}_{k}^{\T} \breve{\boldsymbol{R}}_k^{-1}\breve{\boldsymbol{H}}_{k}}} +
   \underset{\boldsymbol{J}_{\text{prior},k}}{\underbrace{\left(\breve{\boldsymbol{Q}}_k + \breve{\boldsymbol{F}}_{k} \boldsymbol{J}_{k-1}^{-1} \breve{\boldsymbol{F}}_{k}^{\T} \right)^{-1}}}.
\end{equation}
Here $\breve{\boldsymbol{H}}_{k}$ and $\breve{\boldsymbol{F}}_{k}$ are the Jacobian matrices of $\breve{\boldsymbol{h}}(\breve{\boldsymbol{s}}_k)$ and $\breve{\boldsymbol{f}}(\breve{\boldsymbol{s}}_k)$ with respect to $\breve{\boldsymbol{s}}_k$, evaluated at the true state, and $\breve{\boldsymbol{Q}}_k = \text{blkdiag}(\boldsymbol{Q},\boldsymbol{0}_{3I \times 3I})$ and $\breve{\boldsymbol{R}}_k = \text{blkdiag}(\boldsymbol{R}_k^1,\cdots,\boldsymbol{R}_k^I)$ denote the process and measurement noise covariances, respectively. Then, the \ac{PEB} can be computed as $\text{PEB}_k = \sqrt{[ \boldsymbol{J}_{k} ]_{1,1}^{-1} + [ \boldsymbol{J}_{k} ]_{2,2}^{-1}}$, the \ac{HEB} can be computed as $\text{HEB}_k = \sqrt{[ \boldsymbol{J}_{k} ]_{3,3}^{-1}}$, the \ac{CEB} can be computed as $\text{CEB}_k = \sqrt{[ \boldsymbol{J}_{k} ]_{4,4}^{-1}}$, and the \ac{LEB} of $i$-th landmark can be computed as $\text{LEB}_k^i = \sqrt{\sum_{\mu = 3i+2}^{3i + 4} [ \boldsymbol{J}_{k} ]_{\mu,\mu}^{-1}}$.

\subsection{Contribution of Doppler} \label{sec:contribution_Doppler}
From \eqref{eq:fim_recursive}, we can observe that the Doppler of the current measurements only contribute to $\boldsymbol{J}_{\text{data},k}$, as Doppler-related components are only involved in $\breve{\boldsymbol{H}}_{k}$ and $\breve{\boldsymbol{R}}_k$.  As $\breve{\boldsymbol{H}}_{k}$ is the Jacobian matrices of $\breve{\boldsymbol{h}}(\breve{\boldsymbol{s}}_k)$ with respect to $\breve{\boldsymbol{s}}_k$, it can be correspondingly decomposed into blocks
\begin{align}
    \breve{\boldsymbol{H}}_{k}  = \left[ \begin{array}{c : c c c c} \boldsymbol{A}_{k}^{1} & \boldsymbol{B}_{k}^{1}&\boldsymbol{0}_{6 \times 3}  &\cdots & \boldsymbol{0}_{6 \times 3}\\ \boldsymbol{A}_{k}^{2} & \boldsymbol{0}_{6 \times 3} & \boldsymbol{B}_{k}^{2} & \ddots & \boldsymbol{0}_{6 \times 3} \\ \vdots & \vdots& \ddots & \ddots & \vdots  \\ \boldsymbol{A}_{k}^{I} & \boldsymbol{0}_{6 \times 3}& \cdots& \boldsymbol{0}_{6 \times 3}& \boldsymbol{B}_{k}^{I}  \end{array}\right].
\end{align}
Here, $\boldsymbol{A}_{k}^{i}$ and $\boldsymbol{B}_{k}^{i}$ are the Jacobian matrices of $\bar{\boldsymbol{h}}(\boldsymbol{s}_{k},\boldsymbol{x}^{i}_{\text{LM}})$ with respect to $\boldsymbol{s}_{k}$ and $\boldsymbol{x}^{i}_{\text{LM}}$, respectively, which are $\boldsymbol{0}_{6\times4}$ and $\boldsymbol{0}_{6\times3}$ when the $i$-th landmark is not detected, and are
\begin{align}
    \boldsymbol{A}_{k}^{i} =\left[ \begin{array}{c}
    \tilde{\boldsymbol{A}}_{k}^{i}\\ \hdashline  (\tilde{\boldsymbol{a}}_{k}^{i})^{\mathsf{T}}
    \end{array}\right], \quad \boldsymbol{B}_{k}^{i} =\left[ \begin{array}{c}
    \tilde{\boldsymbol{B}}^{i}_{k}\\ \hdashline \vspace{-4mm} \\ (\tilde{\boldsymbol{b}}_{k}^{i})^{\mathsf{T}}
    \end{array}\right],
\end{align}
when the landmark is detected, where $ \tilde{\boldsymbol{A}}_{k}^{i}$ is the Jacobian of $[\tau_{k}^{i},(\boldsymbol{\theta}_{k}^{i})^{\mathsf{T}},(\boldsymbol{\phi}_{k}^{i})^{\mathsf{T}}]^{\mathsf{T}}$ with respect to $\boldsymbol{s}_{k}$, $ \tilde{\boldsymbol{B}}^{i}_{k}$ is the Jacobian of $[\tau_{k}^{i},(\boldsymbol{\theta}_{k}^{i})^{\mathsf{T}},(\boldsymbol{\phi}_{k}^{i})^{\mathsf{T}}]^{\mathsf{T}}$ with respect to $\boldsymbol{x}^{i}_{\text{LM}}$, and
\begin{align}
    (\tilde{\boldsymbol{a}}_{k}^{i})^{\mathsf{T}} & =[
    \begin{array}{cccc}
         \frac{\partial d_{k}^{i}}{\partial x_{k}} & \frac{\partial d_{k}^{i}}{\partial y_{k}} &\frac{\partial d_{k}^{i}}{\partial \varpi_{k}} & \frac{\partial d_{k}^{i}}{\partial b_{k}}
    \end{array}
    ],\\
    (\tilde{\boldsymbol{b}}_{k}^{i})^{\mathsf{T}} &  =[
    \begin{array}{cccc}
        \frac{\partial d_{k}^{i}}{\partial x_{\text{LM}}^{i}} & \frac{\partial d_{k}^{i}}{\partial y_{\text{LM}}^{i}} &\frac{\partial d_{k}^{i}}{\partial z_{\text{LM}}^{i}}
     \end{array}
    ].
\end{align}
\begin{comment}
\begin{align}
    &\boldsymbol{A}_{k}^{i} =\left[ \begin{array}{c}
    \tilde{\boldsymbol{A}}_{k}^{i}\\ \hdashline \vspace{-4mm} \\ (\tilde{\boldsymbol{a}}_{k}^{i})^{\mathsf{T}}
    \end{array}\right]= \left[ \begin{array}{c c c c} \frac{\partial \tau_{k}^{i}}{\partial x_{k}} & \frac{\partial \tau_{k}^{i}}{\partial y_{k}} &\frac{\partial \tau_{k}^{i}}{\partial \varpi_{k}} & \frac{\partial \tau_{k}^{i}}{\partial b_{k}}\\
    \frac{\partial \boldsymbol{\theta}_{k}^{i}}{\partial x_{k}} & \frac{\partial \boldsymbol{\theta}_{k}^{i}}{\partial y_{k}} &\frac{\partial \boldsymbol{\theta}_{k}^{i}}{\partial \varpi_{k}} & \frac{\partial \boldsymbol{\theta}_{k}^{i}}{\partial b_{k}}\\
    \frac{\partial \boldsymbol{\phi}_{k}^{i}}{\partial x_{k}} & \frac{\partial \boldsymbol{\phi}_{k}^{i}}{\partial y_{k}} &\frac{\partial \boldsymbol{\phi}_{k}^{i}}{\partial \varpi_{k}} & \frac{\partial \boldsymbol{\phi}_{k}^{i}}{\partial b_{k}} \vspace{1mm}\\ \hdashline \vspace{1.5mm}
    \frac{\partial d_{k}^{i}}{\partial x_{k}} & \frac{\partial d_{k}^{i}}{\partial y_{k}} &\frac{\partial d_{k}^{i}}{\partial \varpi_{k}} & \frac{\partial d_{k}^{i}}{\partial b_{k}}
    \end{array}\right], \\ \vspace{1mm}
    &\boldsymbol{B}_{k}^{i} =\left[ \begin{array}{c}
    \tilde{\boldsymbol{B}}^{i}_{k}\\ \hdashline \vspace{-4mm} \\ (\tilde{\boldsymbol{b}}_{k}^{i})^{\mathsf{T}}
    \end{array}\right] = \left[ \begin{array}{c c c} \frac{\partial \tau_{k}^{i}}{\partial x_{\text{LM}}^{i}} & \frac{\partial \tau_{k}^{i}}{\partial y_{\text{LM}}^{i}} &\frac{\partial \tau_{k}^{i}}{\partial z_{\text{LM}}^{i}}\\
    \frac{\partial \boldsymbol{\theta}_{k}^{i}}{\partial x_{\text{LM}}^{i}} & \frac{\partial \boldsymbol{\theta}_{k}^{i}}{\partial y_{\text{LM}}^{i}} &\frac{\partial \boldsymbol{\theta}_{k}^{i}}{\partial z_{\text{LM}}^{i}} \\
    \frac{\partial \boldsymbol{\phi}_{k}^{i}}{\partial x_{\text{LM}}^{i}} & \frac{\partial \boldsymbol{\phi}_{k}^{i}}{\partial y_{\text{LM}}^{i}} &\frac{\partial \boldsymbol{\phi}_{k}^{i}}{\partial z_{\text{LM}}^{i}} \vspace{1mm}\\ \hdashline \vspace{1.5mm}
    \frac{\partial d_{k}^{i}}{\partial x_{\text{LM}}^{i}} & \frac{\partial d_{k}^{i}}{\partial y_{\text{LM}}^{i}} &\frac{\partial d_{k}^{i}}{\partial z_{\text{LM}}^{i}}
    \end{array}\right],
\end{align}
when the $i$-th landmark is detected, with $\tilde{\boldsymbol{A}}_{k}$ and $\tilde{\boldsymbol{B}}_{k}$ denoting the parts for delay and angles, and $\tilde{\boldsymbol{a}}^{i}_{k}$ and $\tilde{\boldsymbol{b}}^{i}_{k}$ denoting the parts for Doppler.
\end{comment}
Furthermore, we assume $\boldsymbol{R}_{k}^{i}=\text{diag}(\tilde{\boldsymbol{R}}_{k}^{i},(\sigma_{\text{d},k}^{i})^{2})$.

  According to \eqref{eq:fim_recursive} and the inverse of a block matrix, $\boldsymbol{J}_{\text{data},k}$ is only non-zero in the diagonal block related to the UE state
\begin{align}

&[\boldsymbol{J}_{\text{data},k}]_{1:4,1:4}= \sum_{i=1}^{I}(\boldsymbol{A}_{k}^{i})^{\mathsf{T}}(\boldsymbol{R}_{k}^{i})^{-1}\boldsymbol{A}_{k}^{i}\nonumber\\&=
\sum_{i=1}^{I}\underbrace{(\tilde{\boldsymbol{A}}_{k}^{i})^{\mathsf{T}}(\tilde{\boldsymbol{R}}_{k}^{i})^{-1}\tilde{\boldsymbol{A}}_{k}^{i}}_{\text{non-Doppler related}}+\sum_{i=1}^{I}\underbrace{\tilde{\boldsymbol{a}}_{k}^{i}(\sigma_{\text{d},k}^{i})^{-2}(\tilde{\boldsymbol{a}}_{k}^{i})^{\mathsf{T}}}_{\text{Doppler related}},\label{eq:cramer_rao_bound_elements_data_decomposeA}
\end{align}
and the diagonal blocks related to each landmark
\begin{align}
&[\boldsymbol{J}_{\text{data},k}]_{3i+2:3i+4,3i+2:3i+4}= (\boldsymbol{B}_{k}^{i})^{\mathsf{T}}(\boldsymbol{R}_{k}^{i})^{-1}\boldsymbol{B}_{k}^{i}\nonumber\\&= \underbrace{
(\tilde{\boldsymbol{B}}_{k}^{i})^{\mathsf{T}}(\tilde{\boldsymbol{R}}_{k}^{i})^{-1}\tilde{\boldsymbol{B}}_{k}^{i}}_{\text{non-Doppler related}}+\underbrace{\tilde{\boldsymbol{b}}_{k}^{i}(\sigma_{\text{d},k}^{i})^{-2}(\tilde{\boldsymbol{b}}_{k}^{i})^{\mathsf{T}}}_{\text{Doppler related}}.\label{eq:cramer_rao_bound_elements_data_decomposeB}
\end{align}
Clearly, Doppler provides non-negative information to both UE and landmark states. Due to the block-diagonal structure of
$\boldsymbol{J}_{\text{data},k}$, \eqref{eq:cramer_rao_bound_elements_data_decomposeA} and \eqref{eq:cramer_rao_bound_elements_data_decomposeB} can be interpreted as the equivalent FIMs of the UE and landmark states.

To gain further insights, we expand
$\tilde{\boldsymbol{a}}^i_{k}$ and $\tilde{\boldsymbol{b}}^i_{k}$ as
\begin{align}
    & \tilde{\boldsymbol{a}}_{k}^{i}=\big[-\frac{\left([\boldsymbol{v}_{k}]_{1:2}-d_{k}^{i}[\boldsymbol{q}_{\text{AOA},k}^{i}]_{1:2}\right)^{\mathsf{T}}}{\Vert\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{UE},k}\Vert}, \boldsymbol{v}_{\perp,k}^{\textsf{T}}\boldsymbol{q}_{\text{AOA},k}^{i},0 \big]^{\mathsf{T}}, \\
    & \tilde{\boldsymbol{b}}_{k}^{i} =\frac{\partial \boldsymbol{x}_{\text{inc},k}^{i}}{\partial \boldsymbol{x}_{\text{LM}}^{i}}\frac{\boldsymbol{v}_{k}-d_{k}^{i}\boldsymbol{q}_{\text{AOA},k}^{i}}{\Vert\boldsymbol{x}_{\text{inc},k}^{i}-\boldsymbol{x}_{\text{UE},k}\Vert},
\end{align}
where $\boldsymbol{v}_{\perp,k}=[ -v\sin \varpi_{k}, v\cos \varpi_{k}, 0]^{\mathsf{T}}$ (i.e., $\boldsymbol{v}^\top_{\perp,k}\boldsymbol{v}_{k}=0$). Then,  the additive Doppler-related parts in \eqref{eq:cramer_rao_bound_elements_data_decomposeA} and \eqref{eq:cramer_rao_bound_elements_data_decomposeB} can be summarized as follows:
\begin{itemize}
\item \emph{\ac{ue} position:} as $\boldsymbol{v}_{k}-d_{k}^{i}\boldsymbol{q}_{\text{AOA},k}^{i}=\boldsymbol{v}_{k}-\boldsymbol{v}_{k}^{\textsf{T}}\boldsymbol{q}_{\text{AOA},k}^{i}\boldsymbol{q}_{\text{AOA},k}^{i}$, the  information brought by Doppler is along the direction of the velocity, when the velocity is orthogonal to a path's $\boldsymbol{q}_{\text{AOA},k}^{i}$ direction, i.e., when $\boldsymbol{v}_{k}^{\textsf{T}}\boldsymbol{q}_{\text{AOA},k}^{i}=0$. On the other hand, when the velocity is parallel to a path's $\boldsymbol{q}_{\text{AOA},k}^{i}$ direction, the Doppler of that path does not provide direct Fisher information.  A similar argument holds for landmark locations.

\item \emph{Heading $\varpi_{k}$:} the direct contribution is from the projection of $\boldsymbol{v}_{\perp,k}$ on the direction of $\boldsymbol{q}_{\text{AOA},k}^{i}$ of each path. Hence, if the \ac{ue} velocity is parallel to a path's $\boldsymbol{q}_{\text{AOA},k}^{i}$ direction,
the direct Fisher information  is 0.

\item \emph{Clock bias $b_{k}$:} there is no direct Fisher information contribution. However, the direct contributions to the other dimensions still benefit estimate of the clock bias.
\end{itemize}
Therefore, if the \ac{ue} moves alongside the same or the opposite direction of the  $\boldsymbol{q}_{\text{AOA},k}^{i}$, the corresponding group of measurement  does not  has any Doppler-related information. Moreover, a smaller $(\sigma_{\text{d},k}^{i})^{2}$ results in larger additional Doppler-related information.

\begin{comment}

As $\boldsymbol{J}_{\text{data},k}$ is block-diagonal, the equivalent \ac{FIM} of $\boldsymbol{s}_{k}$ and the $i$-th landmark is the same as \eqref{eq:cramer_rao_bound_elements_data_decompose} and \eqref{eq:cramer_rao_bound_elements_data_decompose2}
\begin{align}
&\boldsymbol{J}_{\text{UE},k}^{\text{E}}=[\boldsymbol{J}_{\text{data},k}]_{1:4,1:4},\\
    &\boldsymbol{J}_{\text{LM},k}^{i,\text{E}}=[\boldsymbol{J}_{\text{data},k}]_{3i+2:3i+4,3i+2:3i+4} \qquad \quad  i=1,\cdots,I.
\end{align}
Clearly, the equivalent \acp{FIM} are non-decreasing, compared to the case without Doppler, as positive semi-definitive parts were added in \eqref{eq:cramer_rao_bound_elements_data_decompose} and \eqref{eq:cramer_rao_bound_elements_data_decompose2}. This is because the complete state remains the same, and an additional dimension of measurement is added. Moreover, a smaller $(\sigma_{\text{d},k}^{i})^{2}$ results in larger $\boldsymbol{J}_{\text{UE},k}^{\text{E}}$ and $\boldsymbol{J}_{\text{LM},k}^{i,\text{E}}$. As the \ac{PIM} is recursively updated by \eqref{eq:fim_recursive}, if $\boldsymbol{J}_{\text{data},k}$ non-decreases, $\boldsymbol{J}_{k}$ also non-decreases, so that the bounds non-increases.

\end{comment}

\section{EK-PMB(M) SLAM filter}
The \ac{slam} framework in this paper follows the EK-PMB \ac{slam} filter proposed in \cite{ge2022computationally}. In this section, the basics of the PMB(M) density and the EK-PMB(M) SLAM filter are summarized in the following.

\subsection{Basics of PMB(M) Density}
The map of the environment $\mathcal{X}$ is formulated as a \ac{PMBM} \ac{rfs}, which can be viewed as the union of two disjoint \acp{rfs}: $\mathcal{X}_{\mathrm{U}}$ for the set of undetected landmarks, which are the landmarks that have never been detected, and $\mathcal{X}_{\mathrm{D}}$ for the set of detected landmarks, which are the the landmarks that have been detected at least once before. Its density follows \cite{garcia2018poisson}
\begin{equation}
    f(\mathcal{X})=\sum_{\mathcal{X}_{\mathrm{U}}\biguplus\mathcal{X}_{\mathrm{D}}=\mathcal{X}}f_{\mathrm{P}}(\mathcal{X}_{\mathrm{U}})f_{\mathrm{MBM}}(\mathcal{X}_{\mathrm{D}}),\label{PMBM}
\end{equation}
where $\mathcal{X}_{\mathrm{U}}$ and  $\mathcal{X}_{\mathrm{D}}$ are modeled as a \ac{PPP} and a \ac{MBM}, respectively, with densities following
\begin{align}
    f_{\mathrm{P}}(\mathcal{X}_{\mathrm{U}})&=e^{-\int\lambda(\boldsymbol{x})\mathrm{d}\boldsymbol{x}}\prod_{\boldsymbol{x} \in \mathcal{X}_{\mathrm{U}} }\lambda(\boldsymbol{x}),\label{PPP} \\
     f_{\mathrm{MBM}}(\mathcal{X}_{\mathrm{D}})&\propto \sum_{j \in \mathbb{I}}w^{j}\sum_{\mathcal{X}^{1}\biguplus \dots \biguplus \mathcal{X}^{|\mathcal{X}_{\mathrm{D}}|}=\mathcal{X}_{\mathrm{D}}}\prod_{i=1}^{|\mathcal{X}_{\mathrm{D}}|}f^{j,i}_{\mathrm{B}}(\mathcal{X}^{i}).\label{MBM}
\end{align}
Here, $\lambda(\cdot)$ denotes the intensity function of the \ac{PPP} density, and  $\mathbb{I}$ is the index set of all global hypotheses with weights satisfying $\sum_{j\in\mathbb{I}}w^{j}=1,w^{j}\ge 0$, and the global hypotheses in \ac{slam} correspond to different \acp{DA} \cite{williams2015marginal}. Each individual component in \eqref{MBM} is termed as a Bernoulli process, and $f_{\mathrm{B}}^{j,i}(\cdot)$ denotes the Bernoulli density of the $i$-th landmark under the $j$-th global hypothesis, following
\begin{equation}
f^{j,i}_{\mathrm{B}}(\mathcal{X}^{i})=
\begin{cases}
1-r^{j,i} \quad& \mathcal{X}^{i}=\emptyset \\ r^{j,i}f^{j,i}(\boldsymbol{x}) \quad & \mathcal{X}^{i}=\{\boldsymbol{x}\} \\ 0 \quad & \mathrm{otherwise},
\end{cases}
\end{equation}
where $f^{j,i}(\cdot)$ denotes the state density, and $r^{j,i} \in [0,1]$ denotes the probability that the corresponding landmark exists.
We note that if there is only one mixture component in the \ac{MBM}, then \eqref{MBM} reduces to a \ac{MB}, and \eqref{PMBM} reduces to a \ac{PMB}.

\subsection{ EK-PMB(M) SLAM Filter Recursion} \label{EK-PMB}
The EK-PMB(M) SLAM filter follows the Bayesian filtering recursion with \acp{rfs}, and it computes the joint posterior $f(\boldsymbol{s}_{k+1},\mathcal{X}|\mathcal{Z}_{1:k+1})$  by \cite{ge2022computationally}
\begin{align}
     f(\boldsymbol{s}_{k+1},\mathcal{X}|\mathcal{Z}_{1:k+1}&) \propto   \ell(\mathcal{Z}_{k+1}|\boldsymbol{s}_{k+1},\mathcal{X}) f(\mathcal{X}|\mathcal{Z}_{1:k}) \nonumber \\ &\times f(\boldsymbol{s}_{k+1}|\mathcal{Z}_{1:k}), \label{joint_posterior}
\end{align}
with $f(\boldsymbol{s}_{k+1}|\mathcal{Z}_{1:k})=\int f(\boldsymbol{s}_{k}|\mathcal{Z}_{1:k})f(\boldsymbol{s}_{k+1}|\boldsymbol{s}_{k}) \text{d} \boldsymbol{s}_{k}$ and $\ell(\mathcal{Z}_{k+1}|\boldsymbol{s}_{k+1},\mathcal{X})$ denoting the RFS likelihood function, given by  \cite[eqs.\,(5)--(6)]{garcia2018poisson}. Instead of tracking the joint density, the filter keeps track of the marginal  \ac{ue} $f(\boldsymbol{s}_{k}|\mathcal{Z}_{1:k})$ and map $f(\mathcal{X}|\mathcal{Z}_{1:k})$ posteriors by marginalizing out the map and the \ac{ue}  state from the joint posterior, respectively. The marginal posteriors are given by
\begin{align}
     f(\boldsymbol{s}_{k+1}|\mathcal{Z}_{1:k+1})  &= \int  f(\boldsymbol{s}_{k+1},\mathcal{X}|\mathcal{Z}_{1:k+1})\delta \mathcal{X},  \label{eq:UpMarObjMarg1}\\
    f(\mathcal{X}|\mathcal{Z}_{1:k+1}) &= \int f(\boldsymbol{s}_{k+1},\mathcal{X}|\mathcal{Z}_{1:k+1}) \mathrm{d}\boldsymbol{s}_{k+1},\label{eq:UpMarObjMarg2}
\end{align}
where $\int \psi(\mathcal{X})\delta \mathcal{X}$ refers to the set integral \cite[eq.~(4)]{williams2015marginal}.

The EK-PMB(M) \ac{slam} filter proposed in \cite{ge2022computationally}  determines $\gamma\ge 1$ most likely \acp{DA} for each prior global hypothesis with corresponding weights, and the joint posterior of the \ac{ue} state and the map conditioned on each \ac{DA} is computed by utilizing an \ac{ekf}-like update step. Then the joint posterior in \eqref{joint_posterior} can be acquired by the weighted summation of the densities for \acp{DA}, followed by \eqref{eq:UpMarObjMarg1} and \eqref{eq:UpMarObjMarg2} to compute the \ac{ue} state and the map. To avoid the exponential increase of \acp{DA}, we can approximate the PMBM density to a PMB density at the end of each update step, which reduces the complexity significantly.

It is important to mention that the weights of $\gamma\ge 1$ most likely \acp{DA} for each prior global hypothesis are computed by using measurements \cite[eq.\,(29)]{ge2022computationally}. Since more information is provided with the involvement of Doppler, finding the correct \acp{DA} becomes more likely,
as correct local hypothesis weights computed in \cite[eqs.\,(22)--(24)]{ge2022computationally} become more prominent.

\section{Results}

\subsection{Simulation Environment}
We consider a 5G downlink scenario with a single \ac{bs} located at $[0 \, \text{m},0 \, \text{m},40\, \text{m}]^{\mathsf{T}} $, 4 \acp{va} located at  $[\pm 200 \, \text{m},0 \, \text{m},40 \, \text{m}]^{\mathsf{T}}$, $[0 \, \text{m}, \pm 200 \, \text{m},40 \, \text{m}]^{\mathsf{T}}$, and 4 \acp{sp} located at  $[\pm 99 \, \text{m},0 \, \text{m},10 \, \text{m}]^{\mathsf{T}}$, $[0 \, \text{m}, \pm 99 \, \text{m},10 \, \text{m}]^{\mathsf{T}}$. The \ac{ue} does a counterclockwise  constant turn-rate movement around the \ac{bs} according to \eqref{eq:dynamic_model}, with $v = 22.22 \text{ m/s}$,  $\omega = \pi/10 \text{ rad/s}$, $T=0.5 \text{ s}$, and $\mathbf{Q}= \text{diag}(0.2^2 \text{ m}^2, \,0.2^2 \text{ m}^2, \,0.001^2 \text{ rad}^2, \,0.2^2 \text{ m}^2)$, and it takes $K = 40$ samples for the \ac{ue} to circle the road once. The measurement covariance matrix is set as  $\boldsymbol{R}_{k}^{i} = \text{blkdiag}(10^{-2} \text{ m}^2, 2.5\times 10^{-3} \cdot \mathbf{I}_4  \text{ rad}^2, \sigma_{\textsf{d}}^{2})$. We initialize the \ac{ue} at $[70.7285 \, \text{m},0 \, \text{m},\pi/2 \, \text{rad}, 300 \, \text{m}]^{\mathsf{T}}$ with covariance as $\text{diag}[ 0.3 \, \text{m}^{2},0.3 \, \text{m}^{2} \, 0.0052 \, \text{rad}^{2},0.3 \, \text{m}^{2}]$. The \ac{bs} is a priori known to the \ac{ue}. We implemented the EK-PMB SLAM filter with considering $\gamma=10$ best \acp{DA} every time step, and compared the results of cases using different levels of $\sigma_{\textsf{d}}$ and without considering Doppler. The mapping  performance is quantified by the \ac{gospa} distance \cite{rahmathullah2017generalized} for both \acp{va} and \acp{sp}, separately, the positioning performance is evaluated by the root mean squared error (RMSE), and we also compare bounds of different cases. More details on parametric settings of the filter can be found in \cite{ge2022computationally}.  The results are averaged over 100 Monte Carlo  simulations.

\subsection{Results and Discussion}
We firstly analyze how different levels of $\sigma_{\textsf{d}}$ affect the bounds. Fig.\,\ref{Fig.bounds} shows how the bounds
change with $\sigma_{\textsf{d}}$, compared with the bounds without considering Doppler. We observe that involving Doppler as a dimension of the measurement improves the positioning and mapping performances of the \ac{slam} system, as the error bounds are lower than in the case where the Doppler measurement is ignored. This effect is stronger when Doppler measurements are more accurate, i.e., when $\sigma_{\textsf{d}}$ is smaller.

\begin{figure}[t]
\center
\input{bounds}\vspace{-5mm}
\caption{The \ac{PEB},  the \ac{HEB}, the \ac{CEB} and the average \ac{LEB} change with $\sigma_{\textsf{d}}$. The benchmarks are the bounds of the case without considering Doppler.}\vspace{-4mm}
\label{Fig.bounds}
\end{figure}

\begin{figure}[t]
\center
\input{mapping}\vspace{-5mm}
\caption{Comparison of mapping performances for VAs between two cases: with and without considering Doppler.}\vspace{-5mm}
\label{Fig.mapping}
\end{figure}

To validate the benefits of Doppler in mmWave radio SLAM,  we then implement the EK-PMB SLAM filter on two cases: 1) involving Doppler into the measurement with $\sigma_{\textsf{d}}=0.1 \text{ m/s}$ and 2) without considering Doppler, and compare the mapping performance for both \acp{va} and \acp{sp}, and the positioning performance between two cases, as illustrated in Fig.\,\ref{Fig.mapping}, Fig.\,\ref{Fig.mapping_SP} and Fig.\,\ref{Fig.bar}, respectively. Fig.\,\ref{Fig.mapping} and Fig.\,\ref{Fig.mapping_SP} demonstrate that the SLAM filter can map the environment for both cases, and the mapping accuracy improves with more measurements being received, as overall, all GOSPA distances gradually decrease over time.  Clearly, when Doppler is involved, the SLAM filter has better mapping performance, as the solid lines are lower than the dashed lines in both Fig.\,\ref{Fig.mapping} and Fig.\,\ref{Fig.mapping_SP}. Fig.\,\ref{Fig.bar} indicates that considering Doppler results in better \ac{ue} state estimates, as lower RMSEs can be acquired. The reasons are that considering Doppler as a dimension of the measurement brings lower bounds. In addition, it also helps the SLAM filter to solve the \ac{DA} problem, which improves the average weight for the correct \ac{DA} from 0.6763 to 0.8622, with weights of all selected $\gamma=10$ \acp{DA} summed to 1 every time step.

\begin{figure}
\center
\input{mapping_SP}
\caption{Comparison of mapping performances for SPs between two cases: with and without considering Doppler.}
\label{Fig.mapping_SP}
\end{figure}

\begin{figure}
\center
\input{Doppler for 5G SLAM/positioning}
\caption{Comparison of \ac{ue} state estimation between two cases: with and without considering Doppler.}
\label{Fig.bar}
\end{figure}

\section{Conclusions}
In this paper, we exploited Doppler as a part of measurement in bistatic mmWave radio SLAM, formulated the EK-PMB SLAM filter on the new measurement model, and provided the \ac{PCRB} for the model. Our results theoretically indicate that the involvement of Doppler helps the SLAM filter to acquire better mapping and positioning performance than the case without considering Doppler. The better Doppler observations are, the lower \acp{PCRB} and the more accurate  map and \ac{ue} state estimates will be. The implementation of the EK-PMB SLAM filter validates the theoretical benefits of involving  Doppler. Future work will include extending the \ac{ue} model to a model with unknown speed, the inclusion of high-dimensional channel estimation, as well as extending the \ac{slam} problem to a \ac{slat} problem.

\scriptsize{
\section*{Acknowledgment}
This work was partially supported by the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by Knut and Alice Wallenberg Foundation, and the Vinnova 5GPOS project under grant 2019-03085, by the Swedish Research Council under grant 2018-03705, and by the Academy of Finland under the grants \#315858, \#328214, \#319994, \#323244, \#346622.
}
\balance

redpajama_set_name: RedPajamaArXiv


Several sites have tried to persuade readers that this is a hoax. It is not. The academic departments referenced in the story have deep ties to the Military/Industrial complex and have issued denials about any observed near earth object possibly striking earth. However they will not explain the specifics of their research in Antarctica. Only recently did they admit to using an military underwater robot and hot water drill that drills through 1km of ice per day. It is the vehemence of the denials and lack of transparency that is becoming very suspicious. If you caught Starfire Tor on CoastToCoastAM last week you may have heard her mention that 'they', the PowersThatBe know that an incoming near earth object is going hit earth sometime in the next 2 years. There is evidence to suggest the object is going to strike one of the poles… most likely the south pole. Special Scientific teams have been down in Antarctica mapping the ice shelf for probable weak points. The object is rumored to be 800 meters wide and when it hits the south pole the entire ice shelf will collapse within months. The PowerThatBe and the super rich are making preparations. They have seed vaults and underground bunkers prepared. This is the reason they do not care about the world's environment or the financial system. It won't matter. A billion will die. The oceans will rise 70 meters. Mankind will lose half the food supply. You have 1 year to prepare. The Mayans, the Sumerians, the Babylonian civilization, native Americans and countless other cultures have foretold this event. The earth is going to get hit by a small object that is in tow of a much larger object. The larger object revisits the earth on a predictable cycle. It will flip the earth upside down. Americans Military forces now control all of the high mountain ranges and have been building bunkers for years, this is the reason for staying in Afghanistan. Only the super rich elite will be invited escape the apocalypse. 90% of the death and destruction will come from the panic and wars after the event. Trillions of Dollars have been sucked out of the economy to build news space shuttles and space projects. Most of the UFO sightings are these black-project vehicles. A University of British Columbia Professor published an online article that projected an 800m asteroid would hit Antarctica in the fall of 2012. His article was on the www.phas.ubc.ca website for 2 days before it abruptly disappeared. The initial data was gathered by The Balloon-borne Large-Aperture Submillimeter Telescope (BLAST) at McMurdo Station, Antarctica. The theorized asteroid was then tracked by Canada-France-Hawaii Telescope on Mauna Kea, which (with the Adaptive Optics Bonnette) supplies probably the sharpest images currently obtainable from the ground. One week after this mysterious article disappeared, "Canadian and American astronauts say the world needs to prepare for the big one — the asteroid impact that could one day devastate the Earth.Veteran Canadian astronaut Chris Hadfield is president of the Association of Space Explorers, which has prepared a detailed report on the asteroid threat. The Canadian Space Agency intends to launch NEOSSat next March to look for asteroids that may be hiding near the sun. The $15-million suitcase-sized satellite will circle about 700 kilometres above the Earth. A Canadian Space Agency official says NEOSSat is expected to detect hundreds of new asteroids during its first year of operation. It will also monitor the heavy traffic of satellites now orbiting the Earth to try to prevent possible collisions." Link The article hypothesized a 94% probability that the asteroid would impact on the Filchner-Ronne ice shelf and cause 2/3 of ice on Antarctica to crumble into sea. The article surmised the total collapse of all Antarctic ice within 2-4 months after impact. "A deep enough impact would crack the ice shelf like a window and total structural collapse would be inevitable, a few months at the outside." If the ice on Antarctica was added the the world's ocean it would raise them by 70 meters.


Furniture — This furniture, is one of the 16 designs published under the gallery: Design Home Office Space Easily. In the next page, we have the design home office space inspiring nifty home office space design home office space photos [ ... ]. While in	the previous page, we have the design home office space photo of fine office space designers modern home office space amazing [ ... ].


Alger Barrett, 78 of Barbourville, Kentucky, formally of Leslie County departed this life on Thursday, May 17, 2018 at the Hazard Appalachian Regional Medical Center in Hazard, Kentucky. He was the son of the late Chester and Chelsea Collett Barrett. He enjoyed playing bluegrass music, fishing, watching westerns, and baseball games. Alger served our country in the United States Army.
Hes was preceded in death by his parents; three brothers Mckinnley, Paul, Wendal Barrett; three sisters Hazel, Virgie and Goldie. Those surviving are sons Bill Barrett and wife Tonya , of Barbourville, Kentucky, Brett Barrett of Barbourville, Ky and step-son James Brock. Daughters are Chantel Bays and Gary, Algina Sizemore, Anna Scott, and step daughter Stephanie Brock Cornett. Five brothers Chester Barrett Jr and Bonnie of Stinnett, Elmer Barrett of Stinnett, Caleb Barrett of Essie, Lynn Barrett and Wilma Barrett of Stinnett, Ky, James and Ethel Barrett of London, Ky.. Fourteen grandchildren and ninteen great-grandchildren A host of other family and friends.
Funeral Services will be conducted at 2:00 P.M. Saturday May 19th, 2018 in the chapel of the Wolfe and Sons Funeral Home with Taylor Adams and James Bowling officiating. Burial will follow in the Bowens Creek Cemetery, Essie, Kentucky. The family will receive friends at the funeral home on Saturday at 11:00 A.M. until the funeral hour.
Wolfe and Sons Funeral Home is in charge of the arrangements.


Tony Pinville & Charles Ollion founded Heuritech in 2013 after they met whilst undergoing their PhD degrees in Artificial Intelligence. At Heuritech, we believe that image recognition is the future of fashion & retail industries, where billions of images are created and shared each day.
Heuritech is building a new kind of visual technology, tailored for the fashion and luxury industry - an empowerment tool for its design, business and marketing teams.
2019 will be a year of rapid scaling for our tools and infrastructure, we are looking for ambitious individuals to join us on this technical challenge !


text:
Firms left baffled by endless new initiatives as Federation of Small Businesses calls for simplification
Published: 17:02 EST, 28 February 2015 | Updated: 05:13 EST, 1 March 2015
John Allan, FSB national chairman
The Coalition has come under fire for creating too many schemes to support businesses.
Firms are baffled by the 'constant churn' of new initiatives, according to the Federation of Small Businesses, which wants them simplified.
On Thursday the Business, Innovation and Skills Select Committee found the provision of business support too 'vague and complex'.
Commenting on the findings, John Allan, FSB national chairman, said: 'While we have seen progress on a more coherent offer nationally through the Business Growth Service and an improved website, there remains some way to go before we see this fully integrated at all levels.
'This is important, as well designed business support can improve firms' productivity. Many small businesses don't understand their options or the support available.
Use the £200m sat on by energy giants to fund advice for... Boss of state-owned Business Bank fights back with £100m...
'Prepare for parental leave or risk being sued': Businesses...
'This is not helped by the constant churn in schemes, which leads to confusion. What's needed is a far more rigorous, focused approach with a clear understanding of businesses' needs.'
Queen's Speech debate: The highlights

redpajama_set_name: RedPajamaCommonCrawl


Thankfully all of these problems can be resolved by turning the video acceleration under Windows 98 down two notches (Select 'Properties' on 'My Computer' –> 'Performance' –> 'Graphics' –> Turn hardware acceleration to one above 'None'). This doesn't negatively affect performance at all.
On a side note – this is one of the reasons why I appreciate the Microsoft WHQL Program (WHQL stands for Windows Hardware Qualification – and is what you need to do to get a driver to load without warnings under Windows 2000 and XP). The WHQL program gives us a complete test suite to run on our Windows 2000 and XP video drivers – that ensures that we never have problems like this – on those platforms.
Unfortunately for Windows 98 we have to test the driver by just running a number of programs and confirming that there are no display problems.
Does reverting back the the pre-SP1 additions also work?
It would address the Hebrew and GDI+ issues – but not the DirectX issue.
Thanks so much for this input. I just installed VM on Thursday nite. This could be the ultimate solution for the company I just started working for. They have a legacy order entry system that runs on a mainframe. It can't work in W2k or higher operating system. I've had some video weirdness during the close of the appplication and the VM window went haywire. But by Right-Alt to go full screen a few times refreshed the VM window and I was able to go on my testing way. We have over 80 computers that will now be able to be upgraded to Windows XP and run their legacy application. I am in awe with VM. Thanks for your input.