Multi-RIS-Enabled 3D Sidelink Positioning

Hui Chen, Member, IEEE,
Pinjun Zheng,
Musa Furkan Keskin, Member, IEEE,
Tareq Al-Naffouri, Senior Member, IEEE
and Henk Wymeersch, Senior Member, IEEE

H. Chen, M. F. Keskin and H. Wymeersch are with the Department of Electrical Engineering, Chalmers University of Technology, 412 58 Gothenburg, Sweden (Email: hui.chen; furkan; henkw@chalmers.se). P. Zheng and  T. Y. Al-Naffouri are with the Division of Computer, Electrical and Mathematical Science & Engineering, King Abdullah University of Science and Technology (KAUST), Thuwal, 23955-6900, KSA. (Email: {pinjun.zheng; tareq.alnaffouri}@kaust.edu.sa). H. Chen and P. Zheng are co-first authors; they contributed equally to this paper.

This work was supported, in part, by the European Commission through the EU H2020 RISE-6G project under grant 101017011, by the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research (OSR) under Award No. ORACRG2021-4695, and by the 6G-Cities project from Chalmers.

Positioning is expected to support communication and location-based services in the fifth/sixth generation (5G/6G) networks. With the advent of reflective reconfigurable intelligent surfaces (RISs), radio propagation channels can be controlled, making high-accuracy positioning and extended service coverage possible.
However, the passive nature of the RIS requires a signal source such as a base station (BS), which limits the positioning service in extreme situations, such as tunnels, dense urban areas, or complicated indoor scenarios where 5G/6G BSs are not accessible. In this work, we show that with the assistance of (at least) two RISs and sidelink communication between two user equipments (UEs), the absolute positions of these UEs can be estimated in the absence of BSs. A two-stage 3D sidelink positioning algorithm is proposed, benchmarked by the derived Cramér-Rao bounds. The effects of multipath and RIS profile designs on positioning performance are evaluated, and localization analyses are performed for various scenarios. Simulation results demonstrate the promising positioning accuracy of the proposed BS-free sidelink communication system. Additionally, we propose and evaluate several solutions to eliminate potential blind areas where positioning performance is poor, such as removing clock offset via round-trip communication, adding geometrical prior or constraints, as well as introducing more RISs.

3D positioning, reconfigurable intelligent surface, 5G/6G, sidelink communication, Cramér-Rao bound.

§ INTRODUCTION

Knowing the position of a mobile device is crucial in various scenarios, starting from emergency services (e.g., disaster rescue) and then extended to daily location-based applications (e.g., navigation and augmented reality) <cit.>.
gnss is one of the most successful positioning systems, which is highly effective in outdoor scenarios but has persistent challenges in dense urban areas and indoor environments.
With the increased frequency and bandwidth of the mmwave/thz systems in the fifth/sixth generation (5G/6G) systems, promising positioning solutions are emerging to complement gnss, and also making communication and positioning integrated and beneficial to each other <cit.>.
Position information can be extracted from channel estimation using radio signals, which can further assist communication with handover <cit.>, and re-establishment of communication links <cit.>. Such integration provides versatile 5G/6G radio systems that can support both communications and positioning functions without introducing extra infrastructure deployments.

Positioning in the 5G new radio (NR) has been studied in TR38.855 <cit.> and more recently in TR38.895 <cit.>, with initial efforts being carried out in both academia and industry. Based on existing mmwave positioning research works, huge potential has been shown in angle-based positioning <cit.>, multipath resolvability <cit.>, positioning under mobility <cit.>, and 6D positioning scenarios <cit.>. Verification and evaluation of the onsite positioning systems have also been carried out with 5G base stations in  indoor <cit.> and outdoor scenarios <cit.>. However, no existing works have reported 5G positioning with comparable performance reported in theoretical analysis or expected in future use cases. The model mismatch (e.g., caused by hardware impairment <cit.>, the effect of multipath <cit.>, erroneous motion model <cit.>) and harsh propagation channels <cit.> constitute major factors that prevent the radio positioning system from achieving high-accuracy performance.
These factors cause errors in channel parameter estimation and further affect the positioning performance, especially for UEs located far away from the BSs, where snr is low and angle estimation error propagates with distance.
Laying out a denser network with more active anchors can mitigate the above-mentioned positioning errors. However, network deployment cost increases with densification (especially for positioning where usually multiple bs are needed at the same time), requiring new enablers to accomplish positioning tasks.

One of the most promising enablers is sidelink communication (or device-to-device communication),  introduced in 3GPP Release 12 <cit.>, and more recently standardized in Release 16 <cit.> to support FR1 and the mmwave range FR2. With direct communication between devices, cooperative positioning is possible, which reduces the requirement for densely deployed BSs. In general, relative position information between each device/vehicle can be obtained in a cooperative positioning network given a sufficient number of vehicles <cit.>. With an anchor provided in a global coordinate system, the true positions of all the devices can be obtained. Moreover, the sidelink can also be implemented in partial coverage and out-of-coverage areas for positioning, where the relative location will be beneficial to vehicles in various applications such as platooning, collision avoidance, and so on <cit.>.

Another promising technology that has been studied extensively for positioning (yet not standardized) is ris <cit.>.[Throughout this work, we consider passive RISs due to their advantage of low power consumption and low deployment cost. However, other types of RIS also exist, such as active RIS <cit.>, hybrid RIS <cit.>, and simultaneous transmissive and receiving (STAR) RIS <cit.>, which are left for future work.] RISs consist of configurable elements with the ability to reshape the channel by changing the phase of the incident signals. For communication, RISs are able to provide improved snr, reduced interference, and extended coverage under blockage. From the positioning point of view, RISs can work as additional passive anchors and provide high-resolution angular information by virtue of a large number of RIS elements. With the assistance of RIS, various positioning scenarios are created, from the simplest scenario where a UE can be localized in a siso system <cit.> to joint localization and mapping with multiple ris <cit.>. In bi-static and multi-static sensing scenarios, when an object is equipped with a RIS, the object can be passively localized with transmitter and receiver anchors <cit.>. More recent works show joint RIS calibration and UE positioning can be performed simultaneously within a simo system, providing a practical solution for RIS calibration <cit.>. All of these works have shown a huge potential for RIS in 5G/6G positioning.

Both sidelink communication and RIS are promising enablers for 5G/6G systems, which have been separately studied in most of the works. The discussion on combining these two technologies for low-latency and high-reliability communication has been discussed in <cit.>, while positioning works appear quite recently <cit.> with a major feature of no bs being involved, resulting in the tx and rx both with unknown positions. In <cit.>, sidelink positioning with RISs is discussed at a high level with localization and sensing scenarios, architectures, and protocols being discussed. Regarding the technical contributions, the work in <cit.> requires the cooperation of multiple UEs with only time-of-arrival information being considered, and self-localization has been studied in <cit.> where the ue is equipped with a full-duplex array. However, the former case did not benefit from the high angular resolution of RISs, and the latter scenario requires extra hardware cost to upgrade half-duplex arrays into (radar-like) full-duplex arrays. More recent work shows that cooperative localization can benefit from spatial frequency estimation with a single RIS. However, at least three UEs are required, and the coordination of multi-UE transmission increases positioning overhead <cit.>.

To the best of our knowledge, this is the first technical work that discusses the multi-RIS-enabled 3D sidelink positioning problem, which aims to locate both the transmitter and the receiver using one-way pilot signal transmission without full-duplex hardware requirements. We will show that with a sufficient number of ris (at least two) involved, the 3D absolute positions of two single-antenna UEs can be estimated using sidelink communication in the absence of BSs, making ubiquitous positioning possible.
Such an application can be used in partial coverage and out-of-coverage scenarios (such as tunnels or rural areas) without any BSs involved <cit.>.
By relying on the deployment of quasi-passive RISs and the cooperation of UEs/vehicles via sidelink communication, localization tasks can be accomplished with low-cost systems.

In this work, we consider a 3D siso sidelink communication scenario with two UEs and several passive RIS anchors. The contributions of this work can be summarized as follows:

* We formulate the problem of multi-RIS enabled 3D siso sidelink positioning. In this scenario, the RISs (at least two) work as passive anchors with known positions and orientations. With sidelink communication, the 3D positions (absolute positions with respect to the global coordinate system) of both UEs and the clock offset between them can be estimated. Note that this problem is fundamentally different from previous positioning works as we do not assume the state of one of the tx or rx to be known as in <cit.>, and we do not assume a full-duplex transceiver is implemented as in <cit.>.

* We derive the crb for both channel parameter estimation and positioning, which serve several purposes: a) to benchmark the proposed positioning algorithms; b) to evaluate the different designs of RIS profiles; c) to provide guidelines on blind areas (where positioning task cannot be completed) evaluation and anchor deployment optimization.

* We adopt a time-orthogonal RIS profile design scheme to assist channel estimation by differentiating the LOS path, and each of the RIS paths from each other. With this scheme, we design positioning-oriented RIS profiles based on directional and derivative codebooks from prior UE information, which can be further improved with power control.

* We develop a low-complexity channel parameter estimation to obtain the delays and spatial frequencies (separate estimation of aoa and aod is not possible in this scenario due to inherent ambiguity, which will be described in Section <ref>). Based on the delay and spatial frequency estimates from multiple RISs, a 3D-search positioning algorithm is developed to estimate the 3D positions of both UEs and the clock offset between them. In addition, maximum likelihood estimators for channel parameter estimation and positioning are also formulated for refining results.

* Extensive simulations are carried out to show the effectiveness of the derived performance analysis and the proposed algorithm. The effects of multipath and RIS profile designs on positioning performance are evaluated. Several RIS deployment strategies (e.g., placed on one side or both sides of the road), and further sidelink positioning system designs are suggested.

The structure of this paper is organized as follows. Section <ref> discusses the system model, based on which problem formulation will be described.
The performance analysis, including the lower bounds for channel parameters and position estimation, is provided in Section <ref>. Section <ref> details the methodology of the RIS profile design and positioning algorithm. Simulation results are presented in Section <ref>, followed by the conclusion of this work in Section <ref>.

Notations and Symbols: Italic letters denote scalars (e.g. $a$), bold lower-case letters denote vectors (e.g. $\av$), and bold upper-case letters denote matrices (e.g. $\Am$). $(\cdot)^\top$, $(\cdot)^\herm$, $(\cdot)^*$, $(\cdot)^{-1}$, $\trace(\cdot)$, and $\Vert{\cdot}\Vert$ represent the transpose, Hermitian transpose, complex conjugate, inverse, trace, and $\ell$-2 norm operations, respectively; $\Am \odot \Bm$, $\Am \otimes \Bm$, $\av \circ \bv$ are the Hadamard product, Kronecker product, and outer product, respectively; $[\cdot,\ \cdot,\ \cdots, \cdot]^\top$ denotes a column vector; $[\cdot]_{i,j}$ is the element in the $i$-th row, $j$-th column of a matrix, and $[\cdot]_{a:b,c:d}$ is the submatrix constructed from the $a$-th to the $b$-th row, and the $c$-th to $d$-th column of a matrix; $\mathrm{Re}\{a\}$ extracts the real part of a complex variable, $\angle(a)$ returns the phase of a complex number $a$; $\bm{1}_{N}$ denotes an $N\times 1$ all ones vector, and $\mathbf{I}_{N}$ denotes a size-$N$ identity matrix.

(image) [anchor=south west];
(image.north east);

at (0.37*,0.5*)[rotate=0,anchor=north]$x$;
at (0.405*,0.685*)[rotate=0,anchor=north]$y$;
at (0.295*,0.77*)[rotate=0,anchor=north]$z$;

at (0.15*,0.8*)[rotate=0,anchor=north]Origin of the Global
Coordinate System:
$[0,0,0,]^\top \unit[]{m}$;
at (0.25*,0.62*)[rotate=0,anchor=north]$\pv_1$;
at (0.52*,0.88*)[rotate=0,anchor=north]$\pv_2$;
at (0.21*,0.13*)[rotate=0,anchor=north]$\pv_\text{T}$ ?;
at (0.80*,0.55*)[rotate=0,anchor=north]$\pv_\text{R}$ ?;
at (0.09*,0.4*)[rotate=0,anchor=north]$\theta_{\text{A},1}$;
at (0.1*,0.58*)[rotate=0,anchor=north]$\phi_{\text{A},1}$;
at (0.68*,0.97*)[rotate=0,anchor=north]$\theta_{\text{D},2}$;
at (0.47*,0.98*)[rotate=0,anchor=north]$\phi_{\text{D},2}$;
at (0.54*,0.42*)[rotate=0,anchor=north]$\Hm_\text{U}$;
at (0.34*,0.47*)[rotate=0,anchor=north]$\Hm_\text{R,1}$;
at (0.57*,0.68*)[rotate=0,anchor=north]$\Hm_\text{R,2}$;
at (0.13*,0.28*)[rotate=0,anchor=north]$\Hm_\text{R,MP,1}$;
at (0.87*,0.72*)[rotate=0,anchor=north]$\Hm_\text{R,MP,2}$;
at (0.57*,0.17*)[rotate=0,anchor=north]$\Hm_\text{U,MP}$;

Illustration of multi-RIS-enabled 3D sidelink positioning. With the help of multiple (at least two) RIS anchors, the positions of both UEs (with respect to the global coordinate system) and the clock offset between them can be estimated through a one-way sidelink communication.

§ SYSTEM MODEL

In this section, we describe the geometry model, signal model, and problem statement of the considered multi-RIS-enabled 3D sidelink positioning.

§.§ Geometry Model

We consider a 3D siso scenario with $L>1$ RISs and two unsynchronized single-antenna ue, where the 3D positions of both UEs need to be estimated via sidelink communication, as shown in Fig. <ref>. The tx and the rx ue are located at $\pv_{\text{T}}$, $\pv_{\text{R}}$ $\in \mathbb{R}^3$, respectively. The positions (array centers) and orientations of $L$ ris are denoted by $\pv_1, \ldots, \pv_L$ $\in \mathbb{R}^{3}$,
and Euler angle vectors $\ov_1, \ldots, \ov_L \in \mathbb{R}^{3}$ (which can be mapped into rotation matrices $\Rm_1, \ldots, \Rm_L$ that belong to the special orthogonal group in three dimensions denoted $\text{SO(3)}$ <cit.>), respectively.

For simplicity, we assume all the ris consist of $N = N_1\times N_2$ RIS elements with $N_1$ and $N_2$ as the number of rows and columns, respectively. In addition, without loss of generality, all the RIS elements are located on the Y-Z plane of each RIS's local coordinate system with the $n$-th element located at $\zv_{n} = [0, z_{n,2}, z_{n,3}]$. The aoa $\vpv_{\text{A},\ell}$ at the $\ell$-th RIS from the tx UE and the aod $\vpv_{\text{D},\ell}$ from the same RIS to the rx UE  can then be expressed as

\begin{align}
\vpv_{\text{A},\ell} & =
\begin{bmatrix}
\phi_{\text{A},\ell}\\
\theta_{\text{A},\ell}
\end{bmatrix}=
\begin{bmatrix}
\arctan2( t_{\text{T},\ell,2}, t_{\text{T},\ell,1})\\
\arcsin( t_{\text{T},\ell,3})
\end{bmatrix},
\label{eq:aoa}
\\
\vpv_{\text{D},\ell} & =
\begin{bmatrix}
\phi_{\text{D},\ell}\\
\theta_{\text{D},\ell}
\end{bmatrix}=
\begin{bmatrix}
\arctan2( t_{\text{R},\ell,2}, t_{\text{R},\ell,1})\\
\arcsin( t_{\text{R},\ell,3})
\end{bmatrix},
\label{eq:aod}
\end{align}

where $\phi$ and $\theta$ are the azimuth and elevation angles, respectively. Let $\tv_{\text{T},\ell} = [ t_{\text{T},\ell,1},  t_{\text{T},\ell,2},  t_{\text{T},\ell,3}]^\top$ and $\tv_{\text{R},\ell} = [t_{\text{R},\ell,1},  t_{\text{R},\ell,2},  t_{\text{R},\ell,3}]^\top$ denote the direction vectors in the local coordinate system of the $\ell$-th RIS to the TX and RX, respectively.  These vectors can be expressed using global positions $\pv_\text{T}$, $\pv_\text{R}$, $\pv_\ell$ and rotation matrix $\Rm_\ell$ as

\begin{align}
\tv_{\text{T},\ell} & =  \Rm_\ell^{-1}\frac{\pv_\text{T} - \pv_\ell}{\Vert \pv_\text{T} - \pv_\ell \Vert}=
\begin{bmatrix}        \cos(\phi_{\text{A},\ell})\cos(\theta_{\text{A},\ell})
\\
\sin(\phi_{\text{A},\ell})\cos(\theta_{\text{A},\ell})\\
\sin(\theta_{\text{A},\ell})\\
\end{bmatrix}
\label{eq:tv_T}
\\
\tv_{\text{R},\ell} & =  \Rm_\ell^{-1}\frac{\pv_\text{R} - \pv_\ell}{\Vert \pv_\text{R} - \pv_\ell \Vert}=
\begin{bmatrix}        \cos(\phi_{\text{D},\ell})\cos(\theta_{\text{D},\ell})
\\
\sin(\phi_{\text{D},\ell})\cos(\theta_{\text{D},\ell})\\
\sin(\theta_{\text{D},\ell})\\
\end{bmatrix}.
\label{eq:tv_R}
\end{align}

We also assume the existence of $M=\sum_{\ell=0}^L M_\ell$ mpc with $M_0$ TX-sp-RX mpc, and $M_\ell$ mpc in the $\ell$-th ris channel (in the forms of TX-SP-RIS-RX or TX-RIS-SP-RX). The corresponding aod, aoa, and direction vectors at the ris with respect to the sp $\pv_{\ell, m}$ can be defined similarly based on (<ref>) to (<ref>).

§.§ Signal Model

Assume $K$ subcarriers are adopted in a wideband system during the sidelink communication, and $G$ ofdm symbols are sent during the coherence time. The received signal block $\Ym\in\mathbb{C}^{K\times G}$ can be formulated as

\begin{equation}
\Ym  = {\Ym_\text{U}} + {\Ym_{\text{R}}} + \Nm,
\label{eq:signal_model}
\end{equation}

where $\Nm\in \mathbb{C}^{K\times G}$ is the additive white Gaussian noise matrix with each element $n_{k,g} \sim \mathcal{CN}(0, \sigma_n^2)$ and $\sigma_n^2 = WN_0$ depending on the bandwidth $W$ and the noise power spectral density (PSD) $N_0$, $\Ym_\text{U}$ and $\Ym_\text{R}$ are the received signal matrix of the uncontrollable paths and RIS paths.
We adopt a geometric channel model (e.g., the Saleh-Valenzuela model <cit.> used for mmWave and THz systems) to account for the multipath effect <cit.>, and the received signals can be modeled as

\begin{align}
\Ym_{\text{U}}
& = (\Hm_{\text{U}} + \Hm_{\text{U,MP}})\odot \Xm \label{eq:channel_LOS_MP}
\\
& = \Big(\underbrace{{\rho_{0}}\Dm(\tau_{0})}_{\text{LOS channel}}+ \underbrace{\sum_{m=1}^{M_0}\rho_{0,m}\Dm(\tau_{0,m})}_{\text{uncontrolled multipath channel}}\Big) \odot \Xm, \notag
\\
\Ym_{\text{R}}
& = \sum_{\ell=1}^{L}(\Hm_{\text{R},\ell} + \Hm_{\text{R,MP},\ell})\odot \Xm
\label{eq:channel_RIS_MP}
\\
& = \sum_{\ell=1}^L \Big(\underbrace{\rho_{\ell}\Dm(\tau_\ell)\odot\Am_\ell(\psiv_\ell)}_{\ell\text{-th RIS channel}} \notag
\\
& + \underbrace{\sum_{m=1}^{M_{\ell}}{\rho_{\ell,m}}\Dm(\tau_{\ell,m})\odot\Am_\ell(\psiv_{\ell,m})}_{\ell \text{-th RIS multipath channel}}\Big)\odot \Xm .\notag
\end{align}

Here, $\Hm_\text{U}$, $\Hm_\text{U, MP}$, $\Hm_\text{R}$  and $\Hm_\text{R, MP}$ are the channel matrix of the LOS path, uncontrolled multipath RIS path, and RIS multipath, respectively. The pilot signal matrix $\Xm$ is defined as

\begin{equation}
\Xm = \sqrt{P}\xv \deltav^\top \in \mathbb{C}^{K\times G}, \ \ \deltav = [\delta_1, \ldots, \delta_G]^\top\ \ (\Vert \deltav \Vert= \sqrt{G}),
\label{eq_transmission_matrix}
\end{equation}

where $\xv \in \mathbb{C}^K$ ($|x_k| = 1$) represents the transmitted symbols for $K$ subcarriers, and the transmission power of the $g$-th transmission is $\delta_g^2 P$, where $\deltav = \bm{1}_G$ indicates a constant transmit power during $G$ transmissions. Here, we use the same $\xv$ for all the $G$ transmissions for simplicity, and the aim of introducing $\deltav$ is to implement a constrained power control for each transmission/RIS beam to enhance positioning accuracy with the same total transmission power, which will be detailed in Section <ref>.
The complex channel gains of the los path, $m$-th uncontrolled multipath, $\ell$-th RIS path and the $m$-th multipath involving the $\ell$-th RIS are denoted as $\rho_{0}$, $\rho_{0, m}$, $\rho_{\ell}$, and $\rho_{\ell, m}$, respectively. The subscripts also apply to the signal propagation delays of different paths such as $\tau_{0}$, $\tau_{0, m}$, $\tau_{\ell}$, $\tau_{\ell, m}$.

The delay matrix $\Dm(\tau) = \dv(\tau) \mathbf{1}_G^\top \in \mathbb{C}^{K\times G}$ contains the delay information of a specific signal propagation path across different subcarriers as[We assume the movement within the coherence time is negligible and hence the delay at the $k$-th subcarrier is identical across different transmissions.]

\begin{equation}
[\Dm(\tau)]_{k,g} = d_k(\tau) = e^{-j2\pi k \Delta_f \tau},
\label{eq:delay_matrix}
\end{equation}

with $\Delta_f=W/K$ as the subcarrier spacing. The delay for the LOS channel $\tau_0$ and for the $\tau_\ell$ RIS channel can be expressed as

\begin{align}
\tau_{0} & = \frac{d_{0} + B}{c}= \frac{{\Vert \pv_\text{T} - \pv_\text{R} \Vert} + B}{c},
\\
\tau_{\ell} & =\frac{d_{\text{T},\ell} + d_{\text{R},\ell} + B}{c}=\frac{{\Vert \pv_\text{T} - \pv_{\ell} \Vert} + {\Vert  \pv_\text{R} - \pv_{\ell} \Vert} + B}{c},
\end{align}

with $B$ indicating the clock offset (converted to meters) between the two UEs. Considering a sp located at $\pv_{0,m}$ affecting the LOS channel and a sp located at $\pv_{\ell,m}$ affecting the $\ell$-th RIS channel, the delay of multipath TX-SP-RX $\tau_{0,m}$ and TX-SP-RIS-RX $\tau_{\ell, m}$ can be expressed as

\begin{align}
\tau_{0,m} & = \frac{d_{\text{T},m} + d_{\text{R},m} + B}{c}
\\
& = \frac{\Vert \pv_\text{T} - \pv_{0, m} \Vert + \Vert \pv_\text{R} - \pv_{0, m} \Vert + B}{c},\notag
\\
\tau_{\ell,m} & = \frac{d_{\text{T},m} + d_{\ell,m} + d_{\text{R},\ell} + B}{c}
\\
& = \frac{\Vert \pv_\text{T} - \pv_{\ell, m} \Vert + \Vert \pv_\ell - \pv_{\ell, m} \Vert + \Vert \pv_\text{R} - \pv_\ell \Vert + B}{c}.\notag
\end{align}

The delay of the TX-SP-RIS-RX multipath can also be defined similarly[Note that the single-bounce MPCs (TX-SP-RX) can be harnessed in MIMO and MISO systems <cit.>. In the SISO system, the delay estimation of the multipath is insufficient to obtain the position of the incidence point, and hence, no extra information will be provided for the UE location. By considering double- or multi-bounce MPCs (e.g., TX-SP-RIS-RX path) in RIS-aided SISO positioning, the performance might be improved. However, it is shown that the improvement is limited due to the severe pathloss, and the detection/association problems are challenging due to the multipath observation of each object via both single-bounce and double-bounce paths <cit.>. As a consequence, we only include MPCs in channel modeling and simulations, but the exploitation of MPCs for positioning is not considered.]. Considering the high pathloss of the RIS channel and the complexity of problem formulation, multi-RIS-involved paths (e.g., TX-RIS1-RIS2-RX path) are not considered throughout this work.

The matrix $\Am_\ell(\psiv) \in \mathbb{C}^{K\times G}$ captures the effect of the $\ell$-th RIS's phase modulation with each element expressed as

\begin{align}
[\Am_\ell(\psiv)]_{k, g} & = a_g(\psiv) = \av(\varphiv_{\text{D}})^\top \Omegam_{\ell, g} \av(\varphiv_{\text{A}})
\label{eq_RIS_equivalent_gain_matrix}
\\
& = \omegav_{\ell, g}^\top (\av(\varphiv_{\text{D}}) \odot \av(\varphiv_{\text{A}})),\notag
\end{align}

where $\Omegam_{\ell, g} = \text{diag}(\omegav_{\ell, g}) \in \mathbb{C}^{N \times N}$ is a diagonal matrix and $\omegav_{\ell, g} = [\omega_{\ell, g, 1}, \ldots, \omega_{\ell, g, N}]$ ($|\omega_{\ell, g, n}|=1$) is a vector containing all the RIS element coefficients.
The steering vectors $\av(\varphiv_\text{A})$ and $\av(\varphiv_\text{D})$ (based on the far-field assumption[Assume a $\unit[30]{GHz}$ SIMO system with a $10\times10$ array is adopted, the Fraunhofer distance can be calculated as $\frac{2D^2}{\lambda} = \frac{2*(\sqrt{2}*0.05)^2}{0.01} = \unit[1]{m}$. Based on <cit.>, the near-field condition for a SISO system with a $10\times 10$ RIS is $\frac{d_1 d_2}{d_1 + d_2} < \frac{2D^2}{\lambda}$ with $d_1$ and $d_2$ as the distances from the RIS to TX and RX, respectively. This indicates that the far-field condition is valid if $d_1>\unit[2]{m}$ and $d_2>\unit[2]{m}$. This work considers TX and RX positions satisfying the far-field condition, and leaves near-field for future work.]) can be expressed as

\begin{align}
\av(\varphiv) & =
e^{j\frac{2\pi f_c}{c} \Zm^\top \tv(\boldsymbol{\varphi})},
\label{eq_av_steer}
\end{align}

with $\Zm = [\zv_1, \ldots, \zv_n] \in \mathbb{R}^{3\times N}$ containing the positions of all the RIS elements (in the local coordinate system of the RIS), and $\tv(\boldsymbol{\varphi})$ can be obtained from (<ref>) and (<ref>).
The aoa and aod of the $l$-th RIS path are denoted as $\psiv_\ell = [\varphiv_{\text{A},\ell}^\top, \varphiv_{\text{D},\ell}^\top]^\top$, defined in (<ref>) and (<ref>), while the multipath aoa/aod vector $\psiv_{\ell, m}$ can be defined similarly based on the position of the incidence point $\pv_{\ell, m}$.

§.§ Spatial Frequency in the RIS Channel

In a RIS-aided positioning system, if ris are used as reflecting anchors for tdoa-based localization (e.g., two-step positioning in <cit.>), at least 6 ris are needed for 3D sidelink positioning due to the unknown tx and rx states. However, by exploiting the angular information at ris, the number of the required RIS anchors can be reduced to 2, which is the scenario considered in this work.
In most of existing ris-aided localization works, the state of bs is known, and hence either $\phiv_\text{D}$ or $\phiv_\text{A}$ is known (e.g., in <cit.>). When full-duplex transceivers are available, such as in self-localization or monostatic sensing scenarios <cit.>, $\phiv_\text{D}$ and $\phiv_\text{A}$ are identical, and hence the channel parameters estimation of the ris path can provide the direction of the ue directly in most existing works.

In this sidelink positioning scenario, however, since the aod and aoa are both unknown, the positioning task is more challenging due to the coupling of aoa and aod as shown in (<ref>).
To facilitate the positioning task, we define a steering vector <cit.>

\begin{equation}
\av_\text{R}(\varphiv_\text{D}, \varphiv_\text{A}) = \av(\varphiv_\text{D}) \odot \av(\varphiv_\text{A}) = e^{j\frac{2\pi f_c}{c} \Zm^\top \tv_\text{R}(\boldsymbol{\varphi}_\text{D}, \boldsymbol{\varphi}_\text{A})},
\end{equation}

\begin{equation}
\tv_\text{R}(\boldsymbol{\varphi}_\text{D}, \boldsymbol{\varphi}_\text{A}) = \tv(\boldsymbol{\varphi}_\text{D}) + \tv(\boldsymbol{\varphi}_\text{A}).
\label{eq_coupled_aoa_aod}
\end{equation}

Note that the first row of the matrix $\Zm$ contains all zeros (RIS elements are located on the local Y-Z plane), meaning the first element of the vector $\tv_\text{R}$ cannot be estimated.[In the scenarios where RIS is not planar (e.g., a cylinder array), all the three elements in $\tv_\text{R}$ can be estimated directly, which is not discussed in this work.] By further defining $\xi$ and $\zeta$ as the second and third entry of the vector $\tv_\text{R}(\boldsymbol{\varphi}_\text{D}, \boldsymbol{\varphi}_\text{A})$ given by

\begin{align}
\label{eq:RIS_aoa_xi}
% spatial frequency 1: \gamma
\xi & = \sin(\phi_{\text{A}})\cos(\theta_{\text{A}}) + \sin(\phi_{\text{D}})\cos(\theta_{\text{D}})
\\
% spatial frequency 2: \chi
\zeta & = \sin(\theta_{\text{A}}) + \sin(\theta_{\text{D}}),
% AOA at RIS: \varphiv_{I}
% \begin{cases}
% \arctan2(t_{\text{B},0,2}, t_{\text{B},0,1}), \tv_{\text{B}, 0} = \frac{\pv_0 - \pv_\text{B}}{d_0} & l = 0
% \\
% \arctan2(t_{\text{B},l,2}, t_{\text{B},l,1}), \tv_{\text{B}, l} = \frac{\pv_{\ell} - \pv_\text{B}}{d_{l,1}} & l > 0
% \end{cases},
% \label{eq:RIS_aoa_ele}
\end{align}

and the matrix $\Am_\ell(\psiv)$ in (<ref>) can also be expressed using spatial frequencies $\xi$ and $\zeta$ as

\begin{equation}
[\Am_\ell(\xi, \zeta)]_{g,k} = \omegav_{\ell, g}^\top \av_\text{R}(\xi, \zeta) = \omegav_{\ell, g}^\top e^{j\frac{2\pi f_c}{c} \Zm^\top [0, \xi, \zeta]^\top}.
\label{eq:RIS_equivalent_gain_xi_zeta}
\end{equation}

These spatial frequencies (i.e., $\xi$, $\zeta$) become channel parameters to be estimated at the ris channel, instead of the aoa and aod. We note from (<ref>) that inherent ambiguities exist in estimating aod and aoa from the spatial frequencies $\xi$ and $\zeta$ since many different aod-aoa pairs can be mapped to a single $\xi$-$\zeta$ pair, making absolute sidelink positioning task challenging.

§.§ Problem Statement

We first define a channel parameter vector as $\etav_{\text{All}} = [\etav^\top, \etav_{\text{MP}}^\top]^\top$, where $\etav = [\etav_{0}, \etav_{1}, \ldots, \etav_{L}]\in \mathbb{R}^{5L+3}$ contains all the parameters related to LOS and RIS channels with $\etav_{0} = [\tau_0, \alpha_0, \beta_0]^\top$ and $\etav_{\ell} = [\xi_\ell, \zeta_\ell, \tau_\ell, \alpha_\ell, \beta_\ell]^\top$, and $\etav_{\text{MP}}$ contains all the multipath-related channel parameters. In the vector $\etav_\text{All}$, $\alpha$ and $\beta$ are the amplitude and phase of the complex channel gain (i.e., $\rho=\alpha e^{-j \beta}$), $\tau$ is the delay, $\xi$ and $\zeta$ are spatial frequencies. Since the mpc do not contribute to positioning and we do not exploit the mapping functionality in this work, we focus on the analysis of $\etav$, and the corresponding state vector $\sv = [\pv_\text{T}^\top, \pv_\text{R}^\top, B, \alpha_0, \beta_0, \ldots, \alpha_L, \beta_L]^\top \in \mathbb{R}^{3L+7}$ containing the 3D positions of two UEs, a clock offset $B$, and the complex channel gains.
To assist analysis and algorithm development, we further define a nuisance-free channel parameter vector $\etav_\text{N} = [\eta_{\text{N},0}, \etav_{\text{N},1}, \ldots, \etav_{\text{N},L}]^\top = [\tau_0, \tau_1, \xi_1, \zeta_1, \ldots, \tau_L, \xi_L, \zeta_L]^\top$ and a nuisance-free state vector $\sv_\text{N} = [\pv_\text{T}^\top, \pv_\text{R}^\top, B]^\top \in \mathbb{R}^{7}$. The vector $\sv_\text{N}$ may have a posterior density function $\text{Prob}(\sv)$ to optimize RIS profile for a better positioning performance, as will be discussed in Section <ref>.

Based on the definitions above, we are able to formulate the sidelink positioning problem as the estimation of the geometric channel parameter vector $\etav$ from the observed signal $\Ym$, and then calculate the state vector $\sv$ based on $\etav$ (detailed in Section <ref> and <ref>).
To make sure the number of channel parameters (i.e., $\etav$ with $5L + 3$ elements) is larger than the number of state parameters (i.e., $\sv$ with $3L+7$ elements), the minimum number of RISs needed is $L=2$ for the los-available. For the los-blockage scenario, at least $L=3$ RISs are needed. Since the positioning task can be performed by a one-way positioning pilot signal transmission, the problem formulation can be extended to multiple UEs, as will be discussed in Section <ref>.
Considering sidelink positioning requires multiple RISs as positioning anchors, the calibration of RISs is needed, and possible solutions can be found in <cit.>, which will not be discussed in this work.

§ LOWER BOUND ANALYSIS

In this section, we derive the crb for the estimation of the channel parameter vector $\etav$ and state vector $\sv$, based on the fim analysis. [In this work, we only consider the fim provided by the data (i.e., received signal) as $\bm{\mathcal{I}} = \bm{\mathcal{I}}_\text{D}$ for performance analysis and positioning algorithm design. If the position prior information is available, the corresponding fim $\bm{\mathcal{I}}_\text{P}$ will only be used for ris profile design in Section <ref>. However, $\bm{\mathcal{I}}_\text{P}$ can also be easily included as $\bm{\mathcal{I}} = \bm{\mathcal{I}}_\text{D} + \bm{\mathcal{I}}_\text{P}$, for performance analysis (e.g., Bayesian crb) or estimator development (e.g., maximum a posteriori estimator).]

§.§ CRB of Channel Parameter Estimation

Based on the defined channel parameter vector $\etav_\text{All}$, state vector $\sv$, and the signal model in (<ref>), (<ref>), (<ref>), the fim of channel parameter estimation can be obtained as <cit.> (Sec. 3)

\begin{align}
\bm{\mathcal{I}}({\boldsymbol\eta_\text{All}}) &
= \frac{2}{\sigma_n^2} \sum^K_{k=1}\sum^{G}_{g=1}\mathrm{Re}\left\{
\left(\frac{\partial{\mu}_{k,g}}{\partial{\boldsymbol\eta_{\text{All}}}}\right)^{\mathsf{H}}
\left(\frac{\partial{\mu}_{k,g}}{\partial{\boldsymbol\eta_{\text{All}}}}\right)\right\}.
\label{eq:FIM_measurement}
\end{align}

Here, $\mu_{k,g} = [\Ym_{\text{U}}]_{k,g} + [\Ym_{\text{R}}]_{k,g}$ is the noise-free observation of the received signal. Note that for fixed system parameters and pilot signals, the $\bm{\mathcal{I}}({\boldsymbol\eta_\text{All}})$ depends on the system geometry (i.e., the position of tx and rx, the position/orientation of ris, and ris profiles.

In contrast to a MIMO system where multipath components can help in localization, sp in the SISO or SIMO/MISO scenarios (e.g., multipath in the LOS and RIS channels of this work) do not affect the localizability of the problem but may degrade the performance of the algorithm (i.e., resulting in unresolvable paths). More discussions on the bounds considering nlos paths can be found in <cit.>.
In this work, we focus on the non-multipath channel parameters $\etav$, whose fim can be easily obtained as the equivalent fim <cit.> of the top-left submatrix as
$\bm{\mathcal{I}}({\boldsymbol\eta}) \in \mathbb{C}^{(5L+3)\times (5L+3)}$.
We can further define delay error bound (DEB) and spatial error bounds (SEBs) for $\tau_\ell, \xi_\ell, \zeta_\ell$ as

\begin{align}
\mathrm{EB}_{\tau_\ell} & = \sqrt{[\bm{\mathcal{I}}({\boldsymbol\eta})^{-1}]_{1+5(\ell-1), 1+5(\ell-1)}}\ , \ \ \ (\ell \ge 0),
\label{eq:DEB}\\
\mathrm{EB}_{\xi_\ell} & = \sqrt{[\bm{\mathcal{I}}({\boldsymbol\eta})^{-1}]_{5\ell-1, 5\ell-1}}\ ,
\ \ \ (\ell>0),
\label{eq:SEB1}\\
\mathrm{EB}_{\zeta_\ell} & = \sqrt{[\bm{\mathcal{I}}({\boldsymbol\eta})^{-1}]_{5\ell, 5\ell}}\ , \ \ \ (\ell>0).
\label{eq:SEB2}
\end{align}

§.§ CRB for 3D Sidelink Positioning

Based on $\bm{\mathcal{I}}({\boldsymbol\eta})$, the CRB of the state parameters $\sv$ can be obtained as

\begin{equation} \label{eq_crb_sidelink}
\mathrm{CRB} \triangleq \left[\bm{\mathcal{I}}(\sv)\right]^{-1} = \left[\Jm_\mathrm{S} \bm{\mathcal{I}}({\etav}) \Jm_\mathrm{S}^\top\right]^{-1},
\end{equation}

where $\Jm_\mathrm{S} \triangleq \frac{\partial {\etav}}{\partial \sv} \in \mathbb{R}^{(3L+7)\times(5L+3)}$ is the Jacobian matrix using a denominator-layout notation from the channel parameter vector $\etav$ to the state vector $\sv$.
We can further define the position error bounds (PEBs), and clock offset error bound (CEB) as

\begin{align}
\mathrm{PEB_\text{T}} & = \sqrt{\trace([\bm{\mathcal{I}}({\sv})^{-1}]_{1:3, 1:3})}
\label{eq:PEBT},\\
\mathrm{PEB_\text{R}} & = \sqrt{\trace([\bm{\mathcal{I}}({\sv})^{-1}]_{4:6, 4:6})}
\label{eq:PEBR},\\
\mathrm{CEB} & = \sqrt{([\bm{\mathcal{I}}({\sv})^{-1}]_{7, 7})}
\label{eq:CEB}.
\end{align}

The derived CRB can be used to assist RIS profile design (see Section <ref>), select weight coefficients in positioning algorithm (see Section <ref>), and evaluate the proposed positioning algorithm as well as RIS profile design performance (see Section <ref>).

§ METHODOLOGY

In this section, we describe RIS profile design strategy, codebooks for the scenarios with and without prior information, channel parameter estimation algorithms, and positioning algorithms.[Note that the effect of multipath is not considered when designing the RIS profile and positioning algorithms. For the architectures/protocols of the RIS-aided sidelink positioning systems and the coordination between the involved devices (e.g., RISs and UEs), please check <cit.> for more details.]

§.§ RIS Profile Design

In the positioning phase, we assume the RIS profiles (or codebooks) are always known at the UE side via RIS-aided positioning protocols <cit.>.
To assist channel parameter estimation, we adopt time-orthogonal RIS profiles to differentiate independent RIS paths from the others <cit.>. We first divide the total transmission $G$ into $\Gamma \ge L+1$ blocks (each block with $\tilde G = G/\Gamma$ OFDM symbols) and define a matrix $\Bm \in \mathbb{C}^{\Gamma \times (L+1)}$  containing orthogonal columns (e.g., from a DFT matrix) as <cit.>

\begin{equation}
% \begin{split}
\Bm = [\bv_0, \bv_1, \ldots, \bv_{L}],\ \ \ \ \text{s.t.}\ \ \Bm^\herm\Bm = \mathbf{I}_{(L+1)\times(L+1)},
% \end{split}
\end{equation}

where each element inside $\Bm$ has a unit amplitude (i.e., $|[\Bm]_{i,j}|=1$). By selecting $\omegav_{\ell,\tilde g} \in \mathbb{C}^{N}$ for $1\le \ell \le L$, and $1\le \tilde g\le \tilde G$, the rest of the RIS profiles can be obtained as

\begin{equation}
\omegav_{\ell, (i-1)\tilde G+\tilde g} = b_{\ell,i}^*\omegav_{\ell, \tilde g},\ (i = 1,\ldots, \Gamma),
\label{eq:orthogonal_profiles}
\end{equation}

where $b_{\ell, i}$ is the $i$-th element of the vector $\bv_{\ell}$. We further define the received signal for the $i$-th block as $\Ym^{(i)}$ ($i = 1, \ldots, \Gamma$), and the LOS path and all the RIS paths can be separated as

\begin{align} \label{eq_ytilde}
\tilde \Ym_\ell & = \frac{1}{\Gamma}\sum_{i=1}^\Gamma b_{ \ell,i}\Ym^{(i)}
% = \frac{1}{\Gamma}\sum_{i}^\Gamma(b_{ \ell,i}\Ym^{(i)}_\ell + \Nm^{(i)})
= \Ym^{(1)}_\ell + \tilde \Nm,
\end{align}

where $\tilde \Ym_\ell \in \mathbb{R}^{K\times \frac{G}{\Gamma}}$ has a smaller size than the received signal block $\Ym \in \mathbb{R}^{K\times G}$ defined in (<ref>), and $\tilde \Nm \in \mathbb{R}^{K \times \tilde G}$ with each element $n_{k,\tilde g} \sim \mathcal{CN}(0, \sigma_n^2/\Gamma)$.

In the following RIS profile design, we will only discuss the design of the first block, and the rest of the blocks can be obtained based on (<ref>) to form orthogonal profiles that assist channel parameter estimation.
Without any prior information on the UE positions, random codebooks are adopted. If the prior information is available, directional or directional and derivative codebooks can be used, as detailed below.

§.§.§ Random Codebook

Without any UE position prior information, a random codebook can be adopted. In this case, each element in the coefficients vector of the $\ell$-th RIS $\omegav_{\ell, g}$ is chosen with unit amplitude and random phase following $\angle \omega_{{\ell, g, n}} \sim \mathcal{U}[0, 2\pi)$.

§.§.§ Directional Codebook

Assume the position prior information of both ue is available (e.g., from previous estimations or based on statistical information), RIS profiles can be designed to improve positioning performance <cit.>. A directional (DIR) codebook is one of the simplest codebooks with the main idea of maximizing the received signal strength of the receiver, given the prior position information of two UEs.
By dropping the time index $g$, for the UEs located at $\pv_\text{T}$ and $\pv_\text{R}$, the optimal RIS profile that maximizes the received energy can be obtained based on (<ref>) and (<ref>) as

\begin{equation}
\omegav^{(1)}_\ell = \av^*(\varphiv_\ell) = e^{-j\frac{2\pi f_c}{c}\Zm^\top (\tv_{\text{T},\ell} + \tv_{\text{R},\ell})}.
\label{eq:RIS_profile_sum_beam}
\end{equation}

Here, $\tv_{\text{T},\ell}$ and $\tv_{\text{R},\ell}$ can be obtained based on (<ref>) and (<ref>), which are the direction vectors obtained from
$\pv_{\text{T}}$ and $\pv_{\text{R}}$, respectively. In real scenarios, however, we cannot know the true location of both UEs, and the prior information may appear in the form of certain distributions
(e.g., the posterior distribution of the nuisance-free state vector $\text{Prob(\sv)}$ with $\sv_\text{N} \sim \mathcal{N}(\bar \sv_\text{N}, \bar \Sigmam_{\sv_\text{N}}$).

In this case, we can sample two candidate positions $\pv_{\text{T},\tilde g}$ and $\pv_{\text{R},\tilde g}$ ($1\le \tilde g \le \tilde G$) for $\tilde G$ times based on $\text{Prob}(\sv_\text{N})$. The DIR codebook of the $\ell$-th RIS (for the first block of transmissions) can be obtained as

\begin{equation}
\Xim_\ell^{\text{DIR}} = [\omegav_{\ell,1}^{(1)}, \ldots, \omegav_{\ell,\tilde G}^{(1)}]  \in \mathbb{C}^{N\times \tilde G},
\label{eq_DIR_codebook}
\end{equation}

with each column $\omegav_{\ell,\tilde g}^{(1)}$ corresponding to the DIR beam sampled UE positions $\pv_{\text{T},\tilde g}$ and $\pv_{\text{R},\tilde g}$ based on (<ref>). The RIS profiles for the rest $(\Gamma-1)$ blocks can be obtained based on (<ref>).

§.§.§ Directional and Derivative Codebook

As has been shown in previous works <cit.>, maximizing the received signal strength at the targeted direction does not necessarily lead to an optimal positioning performance. Hence, the SNR-maximizing directional codebook can be suboptimal from the perspective of positioning. For given $\xi$ and $\zeta$ (computed from the positions $\pv_{\text{T}}$ and $\pv_{\text{R}}$),

the optimal RIS phase profiles (in the sense of minimizing the PEBs in (<ref>) and (<ref>) should lie in the subspace spanned by the following vectors <cit.>,<cit.>:

\begin{align}
\omegav^{(1)} &= \avr^*(\xi, \zeta)
= e^{-j\frac{2\pi f_c}{c}\Zm^\top [0,  \xi, \zeta]^\top},
\label{eq:sum_diff_beam_1}
\\
\omegav^{(2)} &= \frac{\partial \avr^*( \xi,  \zeta) }{\partial \xi}
= \omegav^{(1)} \odot ({-j\frac{2\pi f_c}{c}\Zm^\top [0, 1, 0]^\top}),
\label{eq:sum_diff_beam_2}
\\
\omegav^{(3)} &= \frac{\partial \avr^*( \xi,  \zeta) }{\partial \zeta}
= \omegav^{(1)} \odot ({-j\frac{2\pi f_c}{c}\Zm^\top [0, 0, 1]^\top}),
\label{eq:sum_diff_beam_3}
\end{align}

where the first RIS profile $\omegav^{(1)}$ is identical to the directional (DIR) beam defined in (<ref>), and $\omegav^{(2)}$, $\omegav^{(3)}$ are the so-called derivative (DER) beams. Since the elements in $\omegav^{(2)}$ and $\omegav^{(3)}$ do not have unit amplitude, gradient projection <cit.> is adopted to find the closest unit-amplitude profiles to $\omegav^{(2)}$ and $\omegav^{(3)}$.
The intuition behind the derivative beams is to induce large amplitude changes in response to small perturbations in spatial frequency in the local neighborhood of the true value (see <cit.>). This enables accurate estimation of spatial frequency, along with the high SNR provided by the directional beam.

Similar to the formulation of the DIR codebook in (<ref>), we can sample TX and RX UE positions for $\tilde G/3$ times based on their distribution. The DIR+DER codebook of the $\ell$-th RIS $\Xim_\ell \in \mathbb{C}^{N\times \tilde G}$ can be formulated as

\begin{equation}
\Xim_\ell = [\omegav_{\ell,1}^{(1)}, \omegav_{\ell,1}^{(2)}, \omegav_{\ell,1}^{(3)}, \ldots, \omegav_{\ell,\tilde G/3}^{(1)}, \omegav_{\ell,\tilde G/3}^{(2)}, \omegav_{\ell,\tilde G/3}^{(3)}]  .
\label{eq_DER_codebook}
\end{equation}

Similarly, the RIS profiles for the rest $(\Gamma-1)$ blocks can be obtained based on (<ref>).

§.§.§ Power Control of the DIR+DER Codebook with Prior Information

To further improve the positioning performance, power control can be adopted for the implemented codebook by using different transmit powers for each beam in the codebook.

More specifically, the optimization problem can be formulated (take RX UE for example) as minimizing the expectation of the squared RX PEB (defined in (<ref>)) given by

\begin{align}
\min_{\bm\delta \in \mathbb{R}^G} \int \text{Prob}(\sv)
& \mathrm{PEB^2_\text{R}}(\sv, \Xim_1, \ldots,  \Xim_L|\deltav) \mathrm{d} \sv,     \label{eq_opt_formulation_1}
\\
& \mathrm{s.t.} ~~  \Vert \deltav \Vert ^2 = G ~,\notag
\end{align}

where $\deltav$ is the power control vector defined in (<ref>), and $\text{Prob}(\sv)$ is the posterior distribution of the state vector $\sv$. Considering the high complexity of solving the problem (<ref>), especially when the integral operation is involved, we can further simplify the problem formulation as

\begin{align}
\min_{\bm\delta \in \mathbb{R}^G}
\frac{3}{\tilde G}\sum_{\tilde g = 1}^{\tilde G/3} & \mathrm{PEB^2_\text{R}}(\pv_{\text{T}, \tilde g}, \pv_{\text{R}, \tilde g}, \Xim_1, \ldots,  \Xim_L|\deltav),
\label{eq_opt_formulation_2}
\\
~~ & \mathrm{s.t.} ~~  \Vert \deltav \Vert ^2 = G ~,\notag
\end{align}

where $\pv_{\text{T}, \tilde g}, \pv_{\text{R}, \tilde g}$ are sampled based on the prior information. (<ref>) can be equivalently reformulated as a convex problem
(see Appendix <ref>),
which can be solved using convex optimization tools to provide optimal positioning performance for a given codebook <cit.>. Based on the insights from simulation results, optimal performance can be achieved when DER beams $\omegav^{(2)}$ and $\omegav^{(3)}$ are assigned with the same amount of power. In order to further relieve the computational burden, we propose to use the same power control coefficient $\frac{\sqrt{3}}{\sqrt{1+2\gamma_\text{P}^2}}$ for all the DIR beams (i.e., $\omegav^{(1)}$), and the same coefficient $\frac{\sqrt{3}\gamma_\text{P}}{\sqrt{1+2\gamma_\text{P}^2}}$ for all the DER beams (i.e., $\omegav^{(2)}$, $\omegav^{(3)}$), where $\gamma_\text{P}$ is the ratio between of the DER beam power and DIR beam power (only DIR beams are kept when $\gamma_\text{P}=0$). The optimization problem in (<ref>) can be simplified as

\begin{align}
\min_{\gamma_\text{P} \in \mathbb{R}}
\sum_{\tilde g = 1}^{\tilde G/3} \mathrm{PEB^2_\text{R}}(\pv_{\text{T}, \tilde g}, \pv_{\text{R}, \tilde g}, \Xim_1, \ldots,  \Xim_L|\gamma_\text{P}),   ~~ \gamma_\text{P} \ge 0 ~.
\label{eq_opt_formulation_3}
\end{align}

§.§ Channel Parameter Estimation Algorithm

Once the RIS profiles are designed, the system can send positioning pilot signals and perform positioning algorithms. Here, we describe a two-stage positioning algorithm, including a channel parameter extraction step and a positioning step. For each stage, a coarse estimation algorithm and a refined mle are developed for different performance and complexity tradeoffs.

§.§.§ Low-complexity Channel Parameters Estimator

By implementing the orthogonal RIS profile as described in Section <ref>, the uncontrolled path and each RIS path can be well-separated.
For the LOS path observation $\tilde \Ym_0$ from (<ref>), we first obtain the estimated channel elements and sum across all the $G$ transmissions as

\begin{equation}
\underbrace{\hv_0}_{{\hv_0\in \mathbb{R}^K}} = \sum_{g=1}^{\tilde G} [\tilde\Ym_0]_{:,g}
\odot \xv^*.
\label{eq:coarse_los_channel}
\end{equation}

The delay of the LOS path $\hat \tau_0$ can be estimated based on (<ref>), (<ref>) and (<ref>) as

\begin{equation}
\hat \tau_0 = \arg\max_{\tau} |\dv^\hermit(\tau) \hv_0|,
\label{eq_coarse_los_delay}
\end{equation}

where $\dv(\tau)$ is defined in (<ref>), and (<ref>) can be solved using an $N_\text{F}$ point dft <cit.>.

For the observation of the $l$-th RIS path $\tilde \Ym_\ell$, since RIS profiles are different from one transmission to another, we need to modify (<ref>) and (<ref>) as

\begin{align}
\underbrace{\hat \Hm_\ell}_{\hat \Hm_\ell\in \mathbb{R}^{K\times \tilde G}} & = \tilde\Ym_\ell \odot (\xv^* \mathbf{1}_{\tilde G}^\top), \label{eq:coarse_ris_H} \\
\hat \tau_\ell & = \arg \max_{\tau} \Vert \dv^\hermit(\tau) \hat \Hm_\ell   \Vert. \label{eq_coarse_delay_RIS}
% \hat \kappa_\ell & = \arg \max_k \sum_{g=1}^{\tilde G} [\text{abs}(\fv(\hat \Hm_\ell))]_{k, g}. \label{eq:coarse_ris_k}
\end{align}

Once the delay of the $\ell$-th RIS path $\tau_\ell$ has been obtained, the estimation of spatial frequencies $\hat \xi_\ell$ and $\hat \zeta_\ell$ can be formulated as

\begin{equation}
\label{eq:coarse_ris_xizeta}
[\hat \xi_\ell, \hat \zeta_\ell] = \argmin_{\xi, \zeta} \sum_{g,k} |\omegav_{\ell, g}^\top e^{\frac{j 2\pi}{\lambda_c}\Zm^\top [0, \xi, \zeta]^\top}d_k(\hat \tau_\ell) x_k \tilde y_{k,g}^*|,
\end{equation}

where $\tilde y_{k,g}$ is the element of the matrix $\tilde \Ym$, and the problem can be solved via a 2D search.

§.§.§ MLE for Channel Parameter Estimation

From the low-complexity channel parameters estimator, we can estimate the nuisance-free channel parameter vector $\hat \etav_{\text{N}, \ell}$ ($\ell = 0, 1, \ldots, L$). The MLE aims to find the optimal channel parameters as

\begin{equation}
[\hat \rho_\ell, \hat \etav_{\text{N},\ell}] = \argmin_{\rho_\ell, \boldsymbol{\eta}_{\text{N},\ell}}
\Vert \tilde\muv_\ell - \rho_{\ell}\muv_\ell(\etav_{\text{N},\ell})\Vert,
\label{eq_mle_channel_estimation_0}
\end{equation}

where $\tilde\muv_\ell = \text{vec}(\tilde \Ym_\ell)$, $\muv^\herm(\etav_{\text{N},0}) = \text{vec}(\Dm(\tau_0)\odot\Xm)$, and $\muv^\herm(\etav_{\text{N},\ell}) = \text{vec}(\Dm(\tau_\ell)\odot \Am_{\ell}(\xi_\ell, \zeta_\ell)\odot\Xm)$ ($\ell\ge 1$) that can be obtained from (<ref>) and (<ref>).
Since the channel gain $\rho_\ell$ is a complex constant, by letting $\partial \Vert \tilde\muv_\ell - \rho_{\ell}\muv_\ell(\etav_{\text{N},\ell}) \Vert^2/\partial \rho_{\ell} = 0$, we can obtain the channel gain as $\hat \rho_\ell = \frac{\muv^\herm(\etav_{\text{N},\ell})\tilde\muv}{\Vert \muv(\etav_{\text{N},\ell}) \Vert^2}$.
And hence, the mle can be formulated from (<ref>) with nuisance-free channel parameters only as

\begin{equation}
\hat \etav_{\text{N},\ell} = \argmin_{\etav_{\text{N},\ell}}\left\Vert \tilde \muv_{\ell} - \frac{\muv^\herm(\etav_{\text{N},\ell})\tilde\muv}{\Vert \muv(\etav_{\text{N},\ell}) \Vert^2}\muv_{\ell}(\etav_{\text{N},\ell})\right\Vert.
\label{eq:mle_channel_estimation}
\end{equation}

Channel Parameter Estimation

[1]
— Coarse Estimation —
Input: $\tilde \Ym_0$, $\tilde \Ym_\ell,\ \ell=1,\dots,L.$
Estimate $\hat{\eta}_{\text{N}, 0}=\hat{\tau}_0$ using (<ref>), (<ref>).
$\ell = 1$ to $L$
Estimate $\hat{\tau}_\ell$ using (<ref>) and (<ref>).
Estimate $\hat \xi_\ell$ and $\hat \zeta_\ell$ using 2D grid search in (<ref>).
$\hat{\etav}_{\text{N}, \ell}$ $\leftarrow$ $[\hat{\xi}_\ell, \hat{\zeta}_\ell, \hat{\tau}_\ell]^\top$.
$\hat{\etav}_{\text{N}, \ell}, \ell=0,\dots,L$.
— Refinement —
Input: Coarse estimates $\hat{\etav}_{\text{N}, \ell}, \ell=0,\dots,L$.
$\ell = 0$ to $L$
Obtain the refined $\hat \etav_{\text{N},\ell}$ by solving (<ref>) initialized with $\hat{\etav}_{\text{N}, \ell}$.
refined $\hat{\etav}_\text{N} = [\hat{\eta}_{\text{N}, 0}, \hat{\etav}_{\text{N}, 1}, \ldots, \hat{\etav}_{\text{N}, L}]^\top$.

§.§ Positioning Algorithm

§.§.§ Coarse Position Estimation

Based on the estimated channel parameter vector $\hat\etav$, we propose a 3D-search positioning algorithm to estimate 7 unknowns (i.e., $\pv_\text{T}$, $\pv_\text{R}$, and $B$). For a position candidate $\check \pv_\text{T}$ of the transmitter UE, the candidate direction vector $\check \tv_{\text{T},\ell}$ can be obtained from (<ref>). Based on the estimated spatial frequency $\hat \xi_\ell$ and $\hat \zeta_\ell$, the candidate direction vector of the $\ell$-th RIS $\check \tv_{\text{R},\ell}$ can be calculated as

\begin{equation}
\label{eq_candidate_direction_vec_RX}
\begin{split}
\check t_{\text{R},\ell, 2} & = \hat\xi_\ell - \check t_{\text{T},\ell, 2}\ ,\\
\check t_{\text{R},\ell, 3} & = \hat\zeta_\ell - \check t_{\text{T},\ell, 3}\ ,\\
\check t_{\text{R},\ell, 1} & = \sqrt{1-\check t_{\text{R},\ell, 2}^2-\check t_{\text{R},\ell, 2}^2}\ .
\end{split}
\end{equation}

Note that ambiguities exist in the estimated spatial frequencies due to $\hat \xi, \hat \zeta \in [-1, 1)$, while the true spatial frequencies $\xi, \zeta \in [-2, 2]$. This issue can be solved with prior location information to limit the searching area, or with a reduced RIS inter-element spacing (e.g., to $\lambda_c/4$ instead of $\lambda_c/2$, see <cit.>).

Based on the candidate direction vector $\check \tv_{\text{R},\ell}$ ($\ell \ge 1$) and known RIS states, we are able to calculate the candidate receiver UE position $\check \pv_{\text{R}}$ by getting the closest point to both aod direction vectors <cit.>. Given two bearing lines $\lv_i = \pv_i+r \check \tv_{\text{G},i}$ and $\lv_j = \pv_j+r \check \tv_{\text{G},j}$ ($i,j \in {1, \ldots, L}$ and $\check \tv_{\text{G},\ell} = \Rm_\ell \check \tv_{\text{R},\ell}$), the following equations hold

\begin{align}\label{eq:a1}
\check  \pv_{ij} - \check  \pv_{ji} & = -\check d_{ji}(\check \tv_{\text{G}, j}\times \check \tv_{\text{G}, i}),\\
\check d_{ji} & = \frac{(\check \tv_{\text{G}, j}\times \check \tv_{\text{G}, i})(\pv_j - \pv_i)}{|\check \tv_{\text{G}, j}\times \check \tv_{\text{G}, i}|},
\end{align}

where $\check \pv_{ij}$ is the closest point on the bearing line $\lv_i$ to the bearing line $\lv_j$ denoted as

\begin{equation}
\check \pv_{ij} = \pv_i + \check r_{ij} \check \tv_{\text{G},i}.
\label{eq:closest_points}
\end{equation}

By using least squares, $r_{ij}$ and $r_{ji}$ can be obtained as

\begin{equation}
\begin{bmatrix}
\check r_{ij} \\ \check r_{ji}
\end{bmatrix}
= (\check \Qm^\top \check \Qm)^{-1}\check \Qm^\top [\pv_j - \pv_i -\check d_{ji}(\check \tv_{\text{G}, j}\times \check \tv_{\text{G}, i})],
\end{equation}

with $\check \Qm = [\check \tv_{\text{G},i}, \check \tv_{\text{G},j}]$, and the candidate receiver UE position can be obtained as

\begin{equation}\label{eq:candidate_UE}
\check \pv_\text{R} = \sum_{{i< j}} w_{ij}\check\pv_{ij}, \ \  (\sum_{{i< j}} w_{ij}=1),
\end{equation}

where $i, j\in \{1,\ldots, L\}$ and $w_{ij}$ is the weight coefficient that can be chosen equally as $\frac{2}{L(L-1)}$ if no prior information is provided or based on the quality of the channel parameters estimation (e.g., SNR or PEB using prior information). Since this work is an early-stage work for sidelink positioning, we only perform simulation results (see Section <ref>) to illustrate the effect of different weight coefficients on estimation results. The optimization of weight selection is left for future work.

Finally, we can obtain the estimated clock offset $\check B$ as

\begin{equation}
\label{eq_coarse_position_checkB}
\check B = c \check \tau_0 - \Vert \check \pv_{\text{R}} - \check \pv_{\text{T}} \Vert,
\end{equation}

and the cost function can be formulated as

\begin{equation}\label{eq:coarse_position_costFun}
J(\check \pv_\text{T}) = \sum_{\ell=1}^L w_\ell
|\check B + \Vert \check \pv_{\text{T}} - \check \pv_\ell \Vert + \Vert \check \pv_{\text{R}} - \check \pv_\ell \Vert - c \check \tau_\ell|,
\end{equation}

with $w_\ell$ as the weight coefficient depending on the quality of the delay measurement, which can be set as $1$ by default. Note that the weights mainly affect coarse estimation, but their impacts on the refined results are limited, as will be shown in Section <ref>.
Among all the transmitter UE position candidates, the one with the lowest cost will be the estimated position as

\begin{equation}\label{eq:coarse_position_hatpv}
\hat \pv_\text{T} = \argmin_{\check \pv_\text{T}} J(\check \pv_\text{T}),
\end{equation}

and the rest of the state parameter vector can be obtained based on (<ref>) to (<ref>).

The coarse position estimation algorithm can be easily extended into the los-blockage scenarios by estimating the clock offset from the shortest ris path (take $\ell = 1$ for example) as

\begin{equation}
\check B = c\check \tau_1 - \Vert \check\pv_\text{R} - \pv_\text{1} \Vert - \Vert \check\pv_\text{T} - \pv_\text{1} \Vert,
\end{equation}

and the summation in the cost function (<ref>) starts from $\ell = 2$.

§.§.§ MLE for Positioning

The MLE refinement for positioning can be formulated as

\begin{equation}
\hat \sv_{\text{N}} = \argmin_{\sv_{\text{N}}}\ \left(\hat \etav_{\text{N}} - \etav_{\text{N}}(\sv_{\text{N}})\right)^\mathsf{T}
\Sigmam_{\bm{\eta}_\text{N}}^{-1}\left(\hat \etav_{\text{N}} - \etav_{\text{N}}(\sv_{\text{N}})\right),
\label{eq:mle_position_estimation}
\end{equation}

where $\Sigmam_{\eta_\text{N}} = \bm{\mathcal{I}}({\boldsymbol\eta}_\text{N})^{-1}$ is the covariance matrix of the estimated channel parameters, and the optimization problem in (<ref>) can be solved by, e.g., the trust-region method and the gradient of the cost function in (<ref>) is $-(\frac{\partial \etav_{\text{N}}(\sv_{\text{N}})}{\sv_{\text{N}}})^\mathsf{T}\bm{\Sigma_{\eta_\text{N}}}^{-1}\left(\hat \etav_{\text{N}} - \etav_{\text{N}}(\sv_{\text{N}})\right)$.

For the scenarios covariance matrix in MLE formulation is not available, we can set $\Sigmam_{\bm{\eta}_\text{N}}=\mathbf{I}$, leading to a least squares solution. The pseudo-codes for channel parameter estimation and position estimation can be found in Algorithm <ref> and Algorithm <ref>, respectively.

Position Estimation

[1]
— Coarse Estimation —
Input: Channel parameters $\hat{\etav}_\text{N}$, searching area $\mathcal{A}$.
candidate transmitter $\check \pv_\text{T}\in\mathcal{A}$
Compute $\check \tv_{\text{T},\ell},\ \ell=1,\dots,L$ based on (<ref>).
Compute $\check \tv_{\text{R},\ell},\ \ell=1,\dots,L$ based on (<ref>).
Obtain the candidate UE position $\check \pv_{\text{R}}$ through (<ref>)–(<ref>).
Compute $\check{B}$ using (<ref>).
Calculate the cost $J(\check \pv_\text{T})$ in (<ref>).
Select the optimal $\hat \pv_\text{T}$ that minimizes $J(\check \pv_\text{T})$ as (<ref>) and the corresponding $\hat \pv_\text{R}$ and $\hat B$.
$\hat\pv_\text{T}, \hat\pv_\text{R}, \hat B$.
— Refinement —
Input: coarse estimates $\hat\pv_\text{T}, \hat\pv_\text{R}, \hat B$.
Obtain refined $\hat \sv_\text{N}$ by solving (<ref>) with the initialization $\{\hat\pv_\text{T}, \hat\pv_\text{R}, \hat B\}$.
refined $\hat\sv_\text{N} = [\hat\pv_\text{T}^\top, \hat\pv_\text{R}^\top, \hat B]^\top$.

§.§ Complexity Analysis

In this subsection, we perform complexity analysis on the proposed channel parameter estimation in Section <ref> and positioning algorithms in Section <ref>.
In channel parameter estimation, $L+1$ 1D $N_\text{F}$-point dft are needed for delay estimation, resulting in complexity on the order of $\mathcal{O}_{1}(LN_\text{F}\log N_\text{F})$. For each of the $L$ RISs, a 2D search for spatial frequency estimation is needed, resulting in a complexity of $\mathcal{O}_2(LQ_1 Q_2 GKN)$, where $Q_1=|\mathcal{G}_\xi|$ and $Q_2=|\mathcal{G}_\zeta|$ denote the searching dimension of $\xi$ and $\zeta$, respectively. To refine the channel parameter estimation, we have $\mathcal{O}_3(LQ_3 GKN)$, where $Q_3$ is the number of iterations.
Regarding the positioning algorithm, a 3D search is needed to estimate the positions of both UEs and clock offset, giving $\mathcal{O}_4(L^2Q_4 Q_5 Q_6)$, where $L^2$ indicates the number of beam pairs to be calculated (e.g., there are $L(L-1)$ pairs of beams from $L$ RISs to obtain the candidate receiver position via intersections) and $Q_4, Q_5, Q_6$ represent the searching dimension of the position on the $x$, $y$, and $z$ axis, respectively. For the refinement via MLE, the complexity is $\mathcal{O}_5(L^2 Q_7)$, where $L^2$ indicates the multiplication of matrices, and $Q_7$ is the number of iterations to refine the positioning results.
In summary, the overall complexity of the positioning problem is given by

\begin{align}
& \mathcal{O}_\text{P} = \notag
\\
& \underbrace{\mathcal{O}_1(LN_\text{F}\log N_\text{F}) + \mathcal{O}_2(LQ_1 Q_2 GKN) + \mathcal{O}_4(L^2Q_4 Q_5 Q_6)}_{\text{Coarse Estimation}}\notag
\\
& + \underbrace{\mathcal{O}_3(LQ_3 GKN) + \mathcal{O}_5(L^2 Q_7)}_{\text{Refinement}}.
\label{eq_complexity_analysis}
\end{align}

Based on (<ref>), we can see that the system parameters ($G$, $K$, $N$) affect the channel parameter estimation stage ($\mathcal{O}_2$, $\mathcal{O}_3$), the search grids ($Q_1$, $Q_2$, $Q_4$, $Q_5$, $Q_6$) affect the coarse estimations ($\mathcal{O}_2$, $\mathcal{O}_4$). The number of RISs $L$ increases the complexity for all the positioning steps, and the complexity increases quadratically in the positioning algorithm ($\mathcal{O}_4$, $\mathcal{O}_5$).
With more RISs introduced, the performance will improve based on the same system parameters (e.g., frequency and time resources, transmit power), owing to the extra geometrical information provided, as will be shown in Fig. <ref> and Fig. <ref>. However, the processing delay increases due to the high complexity of the positioning algorithms and possibly the increased $G$ due to the design of time-orthogonal beams.
Considering the dimension of received signal symbols remains the same regardless of the number of RISs (only the channel changes), we anticipate machine learning and deep learning solutions have a huge potential to solve the scalability issue of a larger number of RISs and reduce the processing time of search-based algorithms.

§.§ Extensions to More UEs

Till now, we have discussed the multi-RIS-enabled sidelink positioning with a pair of transmitter and receiver UEs. For the scenarios with more than two UEs, two possible extensions can be considered; both require a coordinator to allocate resources and perform localization tasks <cit.>:

* In the first scenario, one specific UE can be selected as the transmitter that broadcasts positioning reference signals to surrounding UEs. In this case, the structure and the algorithm mentioned in this work can be directly applied, and the positioning process can be performed at the receiver UEs in parallel without interference issues.

* In the second scenario, multiple UEs can play as transmitters and cooperate to obtain their positions (information exchange needed). To avoid interference between transmitting UEs, the time or frequency resources can be properly allocated for different UEs. Although with a higher overhead delay than the first scenario, the cooperation strategy can improve the system performance and can even work with a single RIS anchor with at least 3 UEs (see initial results reported in <cit.>).

Due to the existence of blind areas, the information gains for different UEs from the same RIS will be different, which will be shown in the simulation section in terms of position error bound. In addition, the hardware capabilities and power constraints of different UEs require optimized power allocation and coordination. As a result, dedicated RIS profile optimization, resource allocation, and algorithm hyperparameter selection for multi-UE scenarios are needed in future work.

§ SIMULATION RESULTS

§.§ Simulation Parameters

We consider a 3D scenario with two single-antenna UEs and two RISs. The pilot signal $x_{g,k}$ has a constant amplitude. The unknown channel gains for the LOS path and the $\ell$-th RIS path are set as
$|\rho_0| = \frac{\lambda_c}{4\pi d_0}$

$|\rho_\ell| = \frac{\lambda_c^2}{16\pi^2 d_{\text{T},\ell}d_{\text{R},\ell}}$, both with random phase.
When sp are introduced, the channel gains are set as $|\rho_{0, m}| = \frac{\sqrt{4\pi c_{0,m}} \lambda_c}{16\pi^2 d_{\text{T},m}d_{\text{R},m}}$ and $|\rho_{\ell, m}| = \frac{\sqrt{4\pi c_{\ell, m}} \lambda_c^2}{64\pi^3 d_{\text{T},m}d_{\ell,m}d_{\text{R},\ell}}$, indicating respectively the propagation path of TX-SP-RX and TX-SP-RIS-RX, with $c_{\ell, m}$ as the rcs coefficient of the $m$-th sp in the $\ell$-th path.

For channel parameter estimation, $N_{\text{FFT}} = 2^{10}$ is adopted, and the step size of the 2D grid search is $0.02$. In the positioning step, a step size $\unit[0.2]{m}$ is used to search around the true position (set as a $\unit[2\times2\times2]{m^2}$ area).

The weight coefficients $w_{ij}$ in (<ref>) and $w_\ell$ in (<ref>) are assigned with equal values, and $\Sigmam_{\eta_{\text{N}}}$ in (<ref>) is set as $\mathbf{I}$ by default.
The subcarrier spacing is chosen as $\unit[120]{kHz}$ ($\sim 3300$ subcarriers with a total $\unit[400]{MHz}$ bandwidth), and $K=512$ subcarriers are used for positioning pilot signals.
The rest of the default simulation parameters can be found in Table <ref>.

Default Simulation Parameters

Types     Simulation Parameters

TX Position     $\pv_\text{T} = [-2, -4, 0]^\top$

RX Position     $\pv_\text{R} = [2, 3, 0]^\top$

RIS Positions     $\pv_1 = [-4, 0, 2]^\top$, $\pv_2=[4, 0, 2]^\top$

RIS Orientation     $\ov_1 = [0, 0, 0]^\top$, $\ov_2=[\pi, 0, 0]^\top$

RIS Size     $N = 10\times 10$

Carrier Frequency     $f_c = \unit[30]{GHz}$

Bandwidth     $W = \unit[400]{MHz}$

Number of Transmissions     $G = 192$

Number of Subcarriers     $K = 512$ (out of $3300$)

Clock Offset     $B = \unit[5]{m}$

Average Transmission Power     $P = 30$ dBm

Noise PSD     $N_0 = \unit[-173.855]{dBm/Hz}$

Noise Figure     $\unit[10]{dB}$

(a) Channel parameter estimation

(b) Position estimation
RMSE of the estimation results vs. derived CRBs: (a) channel parameter estimation, (b) position estimation, both are benchmarked by the derived CRBs. Coarse estimation results saturate at high transmission powers, whereas the refined results are able to attach the bounds.

§.§ Channel Parameters and Position Estimation Results

§.§.§ Positioning with a Limited Number of mpc

We first evaluate the performance of the localization algorithm with a single sp located at $[0, 2, 3]^\top\unit[]{m}$ (providing $M=3$ mpc with $c_\text{RCS} = \unit[0.5]{m^2}$). Considering the multi-RIS-enabled 3D sidelink positioning discussed in this work is a novel scenario, no other benchmark algorithms are available. Hence, we use the derived CRB to evaluate the effectiveness of the proposed algorithms. The channel parameter estimation and positioning results are shown in Fig. <ref> (a) and (b), respectively. It can be seen from both figures that the coarse estimations saturate to a certain level with the increased transmit power. However, when refinement processes are applied, the CRBs of channel parameters and state parameters can be attached when the transmit power is higher than $\unit[16]{dBm}$.
Note that when transmit power is low, coarse position results may perform better due to the constrained searching area, while the refinement process does not have such constraints (usually treated as no-information regions).
The results show the effectiveness of the derived bounds and the estimators at the high transmit power. Since refinement processes are involved, the tradeoff between positioning performance and complexity (e.g., the number of iterations) can be performed. In order to reduce the transmit power of the asymptotic region (i.e., $\unit[16]{dBm}$) for specific power-limited UE devices, RIS sizes can be increased, and antenna arrays at the UE side can be utilized for higher beamforming gain. However, the nf channel model and orientation estimations for both UEs must be considered.

Evaluation of the multipath on the estimator in terms of CDF (default multipath parameters are set as $P = \unit[30]{dBm}$, $c_\text{RCS}$ = $\unit[0.5]{m^2}$, $\text{height} = \unit[3]{m}$). We can see that different sp properties (e.g., reflection coefficient, positions) affect positioning performance differently, and large transmit power could help to combat the multipath.

§.§.§ The Effect of Multipath

We further explore the effect of multipath on sidelink positioning by creating two clusters of sp. The first cluster has $5$ sp that are distributed uniformly inside a disk on the x-y plane with a radius of $\unit[1]{m}$ centered at $[0, -3, 3]^\top \unit[]{m}$ (closer to the TX), providing $3\times5=15$ (one LOS and two RIS paths) mpc for the LOS channel and TX-RIS channels. The second cluster is set similarly, centered at $[0, 2, 3]^\top \unit[]{m}$, to affect the LOS channel and RIS-RX channel at the RX side. A total number of $M = 30$ mpc are generated in each realization to affect both the LOS channel and RIS channel. By default, the rcs coefficient is set as $c_{\text{RCS}} = \unit[0.5]{m^2}$ for all the SPs, the phase of each mpc path gain is chosen randomly for different realizations, and the transmit power is set as $P = \unit[30]{dBm}$ to make sure the estimator is working in the asymptotic region (i.e., close to the derived CRB).
Since the position estimation performance of the TX and RX show a similar trend (which makes sense due to a large TX UE estimation error will affect the positioning of the RX UE), we focus on evaluating the estimation error of the RX.

The cdf of the estimation error less than a certain value and the CRB without considering multipath for 200 channel realizations are shown in Fig. <ref>.
When the mpc with a small value of rcs coefficients are introduced, the impact of multipath is limited (see the black curve with diamond markers and the blue curves with circular markers). However, the impact increases with an increased
rcs from $\unit[0.5]{m^2}$ to $\unit[2]{m^2}$ (see red curve with square markers). When changing the height of both cluster centers from $\unit[3]{m}$ to $\unit[2]{m}$, the mpc are likely to make nlos paths non-resolvable, resulting a worse performance under the same rcs coefficient of $\unit[0.5]{m^2}$ (see the green curve with triangle markers). We can also see that a higher transmit power is needed to combat the effect of multipath to achieve satisfactory positioning accuracy. However, power consumption is critical on the UE side, and hence the development of the algorithm taking into account the multipath or the implementation of active ris can be considered for future work.
Considering the randomness of multipath, for simplicity, the simulation results in the following sections do not consider the effect of multipath, which lower bound the performance in the scenarios with multipath.

Evaluation of coarse and fine positioning results of the RX under different weight coefficients. We can see that the equal weight allocation to $w_{ij}$ achieves the best coarse estimation performance. However, even initialized using the worst coarse estimation (i.e., $w_{23}=1$), the refined results are not largely affected. In the refined results, the positions initialized using $w_{ij}=\frac{1}{3}$, $\sigma_\tau/\sigma_{\xi\zeta}=1$ is close to the PEB and perform better than the coefficient selection relying more on delay (i.e., $\sigma_\tau/\sigma_{\xi\zeta}=0.2$ where the delay is assumed to have less error).

(a)-1  (a)-2  (a)-3

(a) DIR beam $\omegav_1$

(b)-1  (b)-2  (b)-3

(b) DER beam $\omegav_2$

The beam pattern of (a) DIR beam $\omegav_1$ and (b) DER beam $\omegav_2$. The figures in column-1 show the normalized radiation patterns of two beams by changing $\xi$ and $\zeta$, and the DIR beam in (a)-1 reaches the maximum at ($\xi = -0.4443$, $\zeta = -0.5039$). The figures in column-2 and column-3 visualize the radiation patterns (in 2D and 3D, respectively) by changing the azimuth and elevation of the receiver $\varphiv_{\text{A},1}$, with a fixed $\varphiv_{\text{A},1} = [-1.1071, -0.2200]^\top$ (rad), the DIR beam in (a)-2 reaches the maximum at $\varphiv_{\text{D},1} = [0.4636, -0.2898]^\top$ (rad).

§.§.§ Evaluation of Weight Coefficients Selection

In the scenarios with 2 RISs, there is only one position candidate $\check{\pv}_\text{R} = \check{\pv}_{12}$ ($w_{12}=1$ in (<ref>)). When extending to 3 RISs, the choice of $w_{12}$, $w_{13}$ and $w_{23}$ affects the positioning performance. We visualize four weight selection results, namely, $w_{12}=1$ (i.e., use $\check{\pv}_{12}$ directly), $w_{13}=1$, $w_{23}=1$, and $w_{12} = w_{23}=w_{13}=\frac{1}{3}$ (denoted as $w_{ij}=\frac{1}{3}$), as shown in Fig. <ref>. We notice that $w_{ij}=\frac{1}{3}$ performs better in the coarse position estimate than other coefficients mentioned before (see solid curves). However, the refined position RMSE with $w_{23}=1$ and $w_{ij}=\frac{1}{3}$ coefficient allocation perform almost the same, with the default $\Sigmam_{\eta_{\text{N}}}=\mathbf{I}$ (denoted by $\sigma_{\tau}=\sigma_{\xi\zeta}=1$). It is also shown that relying more on delay estimation (i.e., assigning $\sigma_{\tau}^2$ to the delay-related entries and $\sigma_{\tau}^2$ to the rest of the entries in the diagonal $\Sigmam_{\eta_{\text{N}}}$) yield worse performance. Note that the fixed weights used in Fig. <ref> are not optimal, and the weight optimization problem is left for future work.

§.§ Evaluation of RIS Profiles

§.§.§ Visualization of DIR and DER beams

Based on the simulation parameters in Table <ref>, we first visualize the radiation patterns (i.e., the equivalent RIS gain $|\omegav^\top \av_\text{R}(\xi, \zeta) |$) of DIR beam $\omegav_1$ and DER beam $\omegav_2$ obtained from (<ref>) and (<ref>) for the first RIS $\pv_{1}$.

By changing the spatial frequencies $\xi$ and $\zeta$, the radiation patterns of two beams are shown in Fig. <ref> (a)-1 and (b)-1. If we assume the position of the TX is known and fix the transmitter angles as $\varphiv_{\text{A},1}$, the 2D radiation patterns of the two beams are visualized in Fig. <ref> (a)-2, (b)-2, and the 3D radiation patterns are visualized in Fig. <ref> (a)-3, (b)-3, respectively. We can see from the figures that the DIR beam maximizes the snr of the TX-RX link, while the DER beams are split at the dimension of $\xi$ and $\zeta$ compared with the DIR beam.
The DER beam $\omegav_3$ (derivation with respect to $\zeta$) shows a similar pattern to Fig. <ref> (b)-1, by splitting the beam from the $\zeta$ axis, which is not plotted.

The evaluation of different prior error levels on positioning with different RIS profile designs. We can see that by using a DIR codebook, the positioning error bounds decrease and then increase when reducing the prior error level, and this effect can be mitigated by power control.

The evaluation of different power allocation coefficients for different prior error levels and RIS sizes. We notice that power control is important when the prior error level is small, and the selection of the optimal coefficient also depends on the RIS size.

§.§.§ The Effect of Prior Error Level on RIS Profile Design

To evaluate the effect of prior error level on the $\text{PEB}_\text{R}$ for different RIS profile designs, we assume the covariance matrices of the prior information are set as $\bar\Sigmam_{\sv_\text{N}} = \sigma^2_\text{pri}\mathbf{I}\in \mathbb{R}^{3\times 3}$, for simplicity. Benchmarked by the random RIS profile (black dashed curve), the PEBs for the DIR codebook, and DIR+DER codebooks with different power allocations are shown in Fig. <ref>. We can see from the figure that both DIR and DIR+DER RIS profiles do not help when the prior error level is high. With more accurate prior information, the DIR profile can largely reduce the PEB. However, when the prior error is too small, the RIS profiles based on the DIR beams are configured to beamforming to a small area and provide less spatial diversity. In an extreme case, the positioning task cannot be completed with all the beams pointing to a single point. When adopting the DIR+DER profiles, however, this phenomenon can be mitigated by choosing a proper power allocation coefficient $\gamma_\text{P}$.

We further evaluate the effect of power allocation coefficients on the positioning error bound, as shown in Fig. <ref>. For a fixed RIS size (i.e., $10 \times 10$), the performance improved by power allocation is limited when the prior error level is high (cyan curve with diamond markers), and becomes crucial with a small error level (blue curve with circle markers). We can also see that the optimal coefficient slightly shifts from left (red triangle) to right (red cross) with the increase of RIS sizes, which is due to the narrow beamwidth requiring accurate prior information.

(a)-1  (a)-2  (a)-3  (a)-4

(a) Scenario 1: Two RISs are parallel facing each other.

(b)-1  (b)-2  (b)-3  (b)-4

(b) Scenario 2: Two RISs are located on the same plane facing the +x axis.

(c)-1  (c)-2  (c)-3  (c)-4

(c) Scenario 3: Two RISs are perpendicular to each other.

Visualization of PEBs of RX at different positions while the TX is fixed (random RIS codebooks). Three scenarios are considered where two RISs (same height) are located $\unit[2]{m}$ above the UEs (same height): a) two RIS are parallel facing each other. b). Two RISs are located parallel to the y-axis. c). Two RISs are forming an L shape facing the positive of the X axis and the negative of the Y axis (e.g., a corner). For each benchmark scenario (column-1), PEBs with known clock offset (column-2), PEBs with known z-axis (column-3), and PEBs with known TX position (column-4) are also visualized.

§.§ Localizability Discussion

§.§.§ PEB visualization of different RIS layouts

Based on the analysis in Section <ref>, at least two RISs are needed to enable sidelink positioning, under the far-field assumption. However, this may not always work, and the localizability also depends on the state of the RIS anchors and the UEs. We visualize the PEBs of the RX UE (with TX UE position fixed) on a 2D x-y plane where both UEs have an unknown but fixed height, RISs are $\unit[2]{m}$ above both UEs, and the TX UE is located at $[-1, -1, 0]^\top$. Without loss of generality, random RIS codebooks and $P = \unit[30]{dBm}$ are adopted. Three different RIS layouts are considered: (a) two RISs at different locations ($[-4, 0, 2]^\top$ and $[4, 0, 2]^\top$) facing each other, (b) two RISs are located on the same y-z plane ($[-4, -3, 2]^\top$ and $[-4, 3, 2]^\top$) facing the same direction, and (c) two RISs are perpendicular to each other ($[-4, 0, 2]^\top$ and $[0, 4, 2]^\top$). The results are shown in Fig. <ref> with different assumptions. Benchmarked by the default setup (column-1), the PEB for perfect synchronization (column-2), UEs with known heights (column-3), and positioning of RX with known TX position (column-4) are visualized. The positions behind the RIS are also plotted, as there exist certain types of RIS that are able to refract the signals rather than reflect them.

We can see from Fig. <ref> that the blind areas exist (yellow area), where the positioning cannot be done or will yield poor performance. However, with extra information, such as clock offset (column-2) or known UE heights (column-3), the blind area can be largely reduced. We also notice that the blind areas in the 3D space are changing continuously (see column-1 and column-4). Since the derivation of PEBs involves a high dimension of parameters such as the RIS positions/orientations and UE positions, it would be challenging to derive a closed-form solution to avoid these areas.

(a)  (b)  (c)  (d)

Interpretation of the blind area. (a) PEB Heatmap of a 2D scenario where RISs and UEs are on the same X-Y plane and the UEs have known heights; (b) Cost function for the RX located at $[-0.61, 2.41, 0]^\top$ (blind area); (c) Cost function for the RX located at $[1.01, 1.01, 0]^\top$ (non-blind area); (d) Candidate TX/RX pairs at the locations around the local minima in (c).

(a)  (b)  (c)  (d)

Visualization of RX PEBs for different numbers of RISs (with $30^\circ$ tilting down). Subfigures (a) and (b) show that the RIS orientations affect the PEB (compared with the benchmarks in Fig. <ref> (a)-1, (c)-1). With more RISs, the areas with poor positioning performance can be largely reduced.

§.§.§ Interpretation of the Blind Area

For a fixed TX UE position, worse performance (yellow area) happens in the location where the surrounding RX positions can provide similar geometrical information. To visualize this, we plot the noise-free cost function, defined in (<ref>), for a 2D scenario. In the localizable location, the cost function shows a clear, sharp global optimal. In the blind area, the cost of the optimal point is similar to the surrounding positions, and the same level of noise will cause a larger estimation error compared with the first scenario. We further choose several surrounding candidate RX UE positions and find the corresponding optimal TX that minimizes the cost function, as shown in Fig. <ref> (d). All these pairs provide similar spatial frequencies and TDOA observations as the ground truth TX/RX pair, and erroneous estimations are likely to happen at these positions.
In order to avoid the effect of local optima, global optimization methods can be adopted. In addition, prior location information can also effectively eliminate the local minima.
As mentioned earlier, for a fixed setup, the blind areas can be reduced by adopting a round-trip estimation to remove the clock offset, using geometric constraints (e.g., known UE height). In the next, we discuss the scenarios with more than two RISs.

§.§.§ Evaluation of More RISs and TX/RX Pairs

Considering the localization performance highly depends on the location of RIS anchors, we can not have a general analytical solution for different numbers of RISs. Here, we numerically evaluate the effect of the number of RISs on positioning with four candidate anchor positions ($[-4, 0, 2]^\top, [0, 4, 2]^\top, [4, 0, 2]^\top, [0, -4, 2]^\top$) covering an $\unit[8\times 8]{m^2}$ area. The RIS tilt angles are set as $30^\circ$ pointing down, and we can see that RIS orientation affects the blind area (see Fig. <ref> (a) vs. Fig. <ref> (c)-1, and Fig. <ref> (b) vs. Fig. <ref> (a)-1). In general, more RISs can increase positioning coverage; however, the improvement depends on the position of UEs and RISs. In addition, if the same orthogonal strategy is implemented, more blocks are needed, increasing difficulties in coordinating between these RISs and channel parameter estimation. To support this statement in a more general scenario, we assume that the TX and RX UEs can be located at grids with $x\in \{-3, -2, \ldots, 3\}\ \text{m}$, $y\in \{-3, -2, \ldots, 3\}\ \text{m}$, and $z\in \{0, 0.5\}\ \text{m}$. We further assume both UEs cannot be located at the same place, and hence a total number of ${{7\times7\times 2}\choose{2}} = 4753$ TX-RX pairs can be evaluated. The cdf of the $\text{PEB}_\text{R}$ for all the scenarios in Fig. <ref> is shown in Fig. <ref>. The CDFs for the RIS (L) scenario with known height (black dashed curve with cross markers) and known clock offset (black dashed curve with triangle markers) are also plotted, which validates the suggested solutions in reducing the blind areas.

CDF of the $\text{PEB}_\text{R}$ for both UEs at different locations inside a $\unit[7\times 7]{m^2}$ area. RIS (L) indicates RISs are located in an L-shape, as shown in Fig. <ref> (a), and RIS (P) indicates RISs are located in parallel, as shown in Fig. <ref> (b).

§.§.§ Summary

We have shown that blind areas exist in the problem of multi-RIS-enabled 3D sidelink positioning with a limited number of RISs (i.e., two), and we have provided several ways to mitigate the effect, namely, round-trip positioning to remove clock bias, consider geometric constraints, and adopt more RISs. These discussions also open new directions for offline and online system optimization. The offline deployment of the anchors needs to consider the TX and RX position probability (e.g., the vehicles can only drive on the road with a certain movement model), as well as the surrounding environment map (e.g., where RISs can be installed). The online optimization needs to take advantage of the prior information and consider when, and to which UE, to trigger a positioning process, as the positioning performance may not meet the positioning performance requirements in the blind areas.

§ CONCLUSION

In this work, we have formulated and solved the multi-RIS-enabled 3D sidelink positioning problem.
In this problem, with the assistance of at least two RISs, the absolute positions of two unsynchronized UEs can be estimated via a one-way sidelink communication in the absence of BSs. Channel parameter estimation and positioning algorithms are developed and benchmarked by the derived CRBs. We discussed the effect of multipath on positioning performance and found the impact on the RIS channels is more significant. We also evaluated the benefit of RIS profile design with prior information to boost positioning performance. Most importantly, we have shown that blind areas exist in RIS-enabled sidelink positioning problems with interpretations. Several solutions can be considered to reduce the effect of blind areas, such as utilizing round-trip communication to remove clock offset, adding geometric constraints to reduce the number of unknowns, and adopting more RISs to increase positioning coverage. However, this work is just the starting point for sidelink positioning with simplified scenarios and channel models. Further directions can consider the high-mobility scenario with the Doppler effect and more accurate channel models that account for more features, such as the near-field effect and beam squint effect.
Moreover, machine learning algorithms can be developed to reduce the algorithm complexity of RIS profile design and weighting coefficients optimization in the positioning algorithm when more RISs are involved.

§ EQUIVALENT CONVEX REFORMULATION OF (<REF>)

To obtain a convex optimization problem from (<ref>), we first express the PEBs in (<ref>) explicitly as a function of $\deltav$: $\mathrm{PEB^2_\text{R}}(\pv_{\text{T}, \tilde g}, \pv_{\text{R}, \tilde g}, \Xim_1, \ldots,  \Xim_L|\deltav) = \trace([\bm{\mathcal{I}}({\sv_{\tilde g}})^{-1}]_{4:6, 4:6})$,

where, using (<ref>) and (<ref>),

\begin{align}
\label{eq_Isvg_}
\bm{\mathcal{I}}(\sv_{\tilde g}) &= \Jm_\mathrm{S} \bm{\mathcal{I}}({\etav_{\tilde g}}) \Jm_\mathrm{S}^\top ~, \\
\bm{\mathcal{I}}({\boldsymbol\eta_{\tilde g}})
&= \frac{2}{\sigma_n^2}\sum^{G}_{g=1} \sum^K_{k=1}\mathrm{Re}\left\{
\left(\frac{\partial{\mu}_{k,g}(\tilde g)}{\partial{\boldsymbol\eta}_{\tilde g}}\right)^{\mathsf{H}}
\left(\frac{\partial{\mu}_{k,g}(\tilde g)}{\partial{\boldsymbol\eta}_{\tilde g}}\right)\right\} ~,
\label{eq:FIM_measurement_r1_} \\ \label{eq:FIM_measurement_nu_r1_}
\mu_{k,g}(\tilde g) &= [\Ym_{\text{U}}(\tilde g)]_{k,g} + [\Ym_{\text{R}}(\tilde g)]_{k,g} ~.
\end{align}

Plugging (6)–(8) into (<ref>), we obtain $\mu_{k,g}(\tilde g) = [\Tm(\tilde g)]_{k,g} \sqrt{P} x_k \delta_g$,

where $\Tm(\tilde g) \triangleq \Hm_{\text{U}} + \Hm_{\text{U,MP}} + \sum_{\ell=1}^{L}(\Hm_{\text{R},\ell} + \Hm_{\text{R,MP},\ell})$.

We now insert this into (<ref>) to obtain $\bm{\mathcal{I}}({\boldsymbol\eta_{\tilde g}})
= \sum^{G}_{g=1} \gamma_g \Om_g(\tilde g) $

for some known matrix $\Om_g(\tilde g) $,  where $\gamma_g \triangleq \delta_g^2$.

Substituting this into (<ref>) yields

\begin{align}
\bm{\mathcal{I}}({\sv_{\tilde g}})
&= \sum^{G}_{g=1} \gamma_g \Sm_g(\tilde g) ~,
\label{eq:FIM_measurement_r1_v3_}
\end{align}

where $\Sm_g(\tilde g) \triangleq \Jm_\mathrm{S} \Om_g(\tilde g) \Jm_\mathrm{S}^\top $. Performing change of variables $\bm\gamma = [\gamma_1 \, \ldots \, \gamma_G]^\top$, the problem (<ref>) can be expressed as

\begin{align}
\min_{\bm\gamma \in \mathbb{R}^G}
\frac{3}{\tilde G}\sum_{\tilde g = 1}^{\tilde G/3} & \trace([\bm{\mathcal{I}}({\sv_{\tilde g}})^{-1}]_{4:6, 4:6})
\label{eq_opt_formulation_3_r1_}
~ & \mathrm{s.t.} ~~  \boldone^\top \bm\gamma = G ~.%\notag ~.
\end{align}

Introducing auxiliary variables $u_{\tilde g}$, (<ref>) can be recast as

\begin{align}
\min_{\bm\gamma \in \mathbb{R}^G, \{u_{\tilde g}\}}
&~~ \frac{3}{\tilde G}\sum_{\tilde g = 1}^{\tilde G/3} u_{\tilde g}
\label{eq_opt_formulation_4_r1_}
\\ \label{eq_const_g_}
~~  \mathrm{s.t.} &~~  \boldone^\top \bm\gamma = G  \,, \,
\trace([\bm{\mathcal{I}}({\sv_{\tilde g}})^{-1}]_{4:6, 4:6}) \leq u_{\tilde g}, \forall \tilde g ~.
\end{align}

From <cit.>, each constraint in (<ref>) can be equivalently reformulated as a linear matrix inequality (LMI) constraint, leading to (using (<ref>))

\begin{align}
\min_{\bm\gamma \in \mathbb{R}^G, \{u_{\tilde g}\}, \{u_{\tilde g,\ell} \}}
&~~ \frac{3}{\tilde G}\sum_{\tilde g = 1}^{\tilde G/3} u_{\tilde g}
\label{eq_opt_formulation_5_r1_}
\\
~~  \mathrm{s.t.} &~~  \boldone^\top \bm\gamma = G \notag ~,
\\ \nonumber
& ~~ \begin{bmatrix}
\sum^{G}_{g=1} \gamma_g \Sm_g(\tilde g) & \ev_\ell \\
\ev_\ell^\top & u_{\tilde g,\ell}
\end{bmatrix} \succeq 0, \ell = 4,5,6,
\\ \label{eq_const_g2_} &~ \sum_{\ell = 4}^{6} u_{\tilde g,\ell} \leq  u_{\tilde g}, \forall \tilde g ~,
\end{align}

where $u_{\tilde g,\ell}$ denote the newly introduced auxiliary variables. The problem (<ref>), which is equivalent to (<ref>), is a convex semidefinite programming (SDP) problem <cit.>.