# Quantum simulation of hadronic states with Rydberg-dressed atoms

Zihan Wang<EMAIL_ADDRESS>Blackett Laboratory, Imperial College
London, London SW7 2AZ, United Kingdom  Feiyang Wang  Blackett Laboratory,
Imperial College London, London SW7 2AZ, United Kingdom  Joseph Vovrosh
Blackett Laboratory, Imperial College London, London SW7 2AZ, United Kingdom
PASQAL SAS, 2 av. Augustin Fresnel, 91120 Palaiseau, France  Johannes Knolle
Blackett Laboratory, Imperial College London, London SW7 2AZ, United Kingdom
Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4,
D-80799 München, Germany  Department of Physics TQM, Technical University of
Munich, 85748 Garching, Germany  Florian Mintert  Blackett Laboratory,
Imperial College London, London SW7 2AZ, United Kingdom  Helmholtz-Zentrum
Dresden-Rossendorf, Bautzner Landstraße 400, 01328 Dresden, Germany  Rick
Mukherjee<EMAIL_ADDRESS>Blackett Laboratory, Imperial
College London, London SW7 2AZ, United Kingdom  Center for Optical Quantum
Technologies, Department of Physics, University of Hamburg, Luruper Chaussee
149, 22761 Hamburg, Germany

###### Abstract

The phenomenon of confinement is well known in high-energy physics and can
also be realized for low-energy domain-wall excitations in one-dimensional
quantum spin chains. A bound state consisting of two domain-walls can behave
like a meson, and in a recent work of Vovrosh et al. [PRX Quantum 3, 040309
(2022)] , it was demonstrated that a pair of mesons could dynamically form a
meta-stable confinement-induced bound state (consisting of four domain-walls)
akin to a hadronic state. However, the protocol discussed in Vovrosh et al.
[PRX Quantum 3, 040309 (2022)] involving the use of interactions with
characteristically non-monotonic distance dependence is not easy to come by in
nature, thus, posing a challenge for its experimental realization. In this
regard, Rydberg atoms can provide the required platform for simulating
confinement-related physics. We exploit the flexibility offered by interacting
Rydberg-dressed atoms to engineering modified spin-spin interactions for the
one-dimensional transverse field Ising model. Our numerical simulations show
how Rydberg-dressed interactions can give rise to a variety of effective
potentials that are suitable for hadron formation, which opens the possibility
of simulating confinement physics with Rydberg platforms as a viable
alternative to current trapped-ion experiments.

Introduction.— Confinement is one of the cornerstones of high-energy physics.
In quantum chromodynamics (QCD), the asymptotic freedom of the strong
interaction implies that the interaction energy between a pair of elementary
particles such as quarks grows with increasing separation Greensite (2011). As
a result of this, quarks remain confined within a certain distance forming
composite objects such as hadrons. Since high-energy experiments are extremely
costly and resource intensive, there is increasing interest in quantum
simulations that can guide experimental activities Bañuls  _et al._ (2020);
Surace  _et al._ (2020); Yang  _et al._ (2020); Zhou  _et al._ (2022).

Interestingly, various aspects of confinement physics are realized in
condensed matter systems, in particular, one-dimensional quantum spin chains
McCoy and Wu (1978); Delfino  _et al._ (1996); Kormos  _et al._ (2017); Liu
_et al._ (2019). In such models, the initial state acts as a source of
quasiparticle excitations that are well described by doman-walls which are
regions of up (down) spins followed by a region of down (up) spins as shown
schematically in Fig. 1(b). It has been shown that bound states of domain-
walls which behave like mesons (two-domain-wall state) can appear in the
transverse field Ising model (TFIM) with a longitudinal field Kormos  _et al._
(2017) or with long-range interactions Liu  _et al._ (2019). Although
confinement in one-dimensional spin chains is much simpler than confinement
physics in higher dimensions studied in QCD, it is nevertheless a noteworthy
step in quantum simulation.

Figure 1: (a) Setup involving an array of trapped Rydberg-dressed atoms
simulating the spin-1/2 system: Each atom has a four-level structure with two
ground states ($\ket{g_{1}},\ket{g_{2}}$) coupled to Rydberg states
($\ket{r_{1}},\ket{r_{2}}$) with a pair of excitation lasers such that
$\ket{g_{i}}$ only couples to $\ket{r_{j\neq i}}$ with Rabi frequency
$\Omega_{j}$ and detuning $\Delta_{j}$. Atoms interact via the Rydberg-Rydberg
interaction, which is highlighted in the top orange box denoted by
$\hat{H}_{\mathrm{vdw}}$. (b) Dynamical formation of a tetraquark during a
collision event between mesons: Numerical results showing dynamics of two
$1$-meson wavepackets (meson of width 1 site as shown schematically) where the
majority of the wavepackets are reflected elastically, but with a small and
finite probability to form a tetraquark at the center. Here $P_{i}$ measures
the probability that a domain-wall is located at site $i$. Parameters used for
the numerical simulation are discussed in Sec. 5 of the Appendix.

The formation of meson bound states in one dimensional quantum spin models has
interesting consequences for non-equilibrium physics Kormos  _et al._ (2017),
the propagation of energy Mazza  _et al._ (2019), and quantum information
Vovrosh and Knolle (2021); Vovrosh  _et al._ (2021). Although scattering
processes between mesons are predominantly elastic, the formation of a stable
multi-meson bound state such as a tetra-quark (four-domain wall state)
requires deep inelastic collisions Milsted  _et al._ (2022); Karpov  _et al._
(2022); Surace and Lerose (2021). One way to enhance the overlap of the
initial asymptotic scattering free meson states with the hadronic bound states
is by quenching the kinetic energy of the individual mesons. This approach
requires a high degree of precision and control in the temporal variation of
the external field. An alternative approach to obtain hadronic bound states is
to fuse a pair of mesons by inducing non-monotonic meson-meson interaction
Vovrosh  _et al._ (2022). However, realizing non-monotonic spin interactions
in conventional simulators, such as trapped ions or cold atoms, is challenging
in general.

In this article, we detail the microscopic implementation of dynamical hadron
formation on a Rydberg quantum simulator and identify parameters suitable for
the experimental realization of the protocol discussed in Ref. Vovrosh  _et
al._ (2022). Rydberg dressing Che  _et al._ (2015); Zeiher  _et al._ (2016);
Balewski  _et al._ (2014); Young  _et al._ (2018); Madjarov  _et al._ (2020);
Henkel  _et al._ (2010); Keating  _et al._ (2013); Petrosyan and Mølmer
(2014); Glaetzle  _et al._ (2015); Li  _et al._ (2018); Wu  _et al._ (2021);
Gil  _et al._ (2014); van Bijnen and Pohl (2015); Mukherjee (2019); Keleş  _et
al._ (2020) is a well-known and useful concept that allows for a small
admixture of Rydberg state to the ground state atoms, which in turn provides
unique flexibility to manipulate the strength and shape of the interactions
using parameters that are easily tunable in experiments. The goal is to
exploit these highly tunable anisotropic Rydberg-dressed interactions to
provide the necessary non-monotonic potential that will facilitate the
formation of the hadron during the mesonic dynamics.

Figure 2: Schematic diagrams outlining the different possible channels during
adiabatic elimination of the different Rydberg states belonging to the same
Zeeman manifold in the weak coupling limit: Second-order processes contain
only singly excited states as intermediate channels as shown in (a), while
fourth-order processes also include eigenstates $\ket{E_{i=0,1,2,3}}$ which
are the result of the Rydberg interaction induced Zeeman mixing when both
atoms are in Rydberg states. $\\{i,j,k,l\\}\in\\{1,2\\}$ with bar representing
the compliment of value.

Setup.— Consider a one-dimensional array of trapped atoms with each atom
possessing a pair of ground states $\ket{g_{1}}$ and $\ket{g_{2}}$ that are
weakly coupled to Rydberg states $\ket{r_{2}}$ and $\ket{r_{1}}$ respectively
using excitation lasers that are characterized by Rabi frequencies
$\Omega_{i=1,2}$ and detunings $\Delta_{i=1,2}$ as shown in Fig. 1(a). A
magnetic field along the direction of the laser is applied shown as $B$ in
Fig. 1(a) which defines the quantization axis and by choosing (left/right)
circularly polarized excitation lasers, one obtains the level scheme described
in Fig. 1(a). The indices $\\{i,j,l,k\\}$ label the internal levels for the
individual atoms and not the position of the atoms. The ground state energy
$E_{g}$ is set to zero without loss of generality. A pair of atoms separated
by distance $\mathbf{r}$ can interact only when they are both excited to their
Rydberg states. Focusing on this two-atom picture, the Hamiltonian that
describes the laser excitation process along with the atom-atom interaction is
given as

$\displaystyle\hat{H}=\hat{H}_{0}(\mathbf{r})+\hat{H}^{\prime},$ (1)
$\displaystyle\hat{H}_{0}(\mathbf{r})=\sum_{i=1}^{2}-\Delta_{i}\ket{r_{i}}\bra{r_{i}}+\hat{H}_{\mathrm{vdw}}(\mathbf{r}),$
$\displaystyle\hat{H}^{\prime}=\sum_{i,j=1;i\neq
j}^{2}\left(\frac{\Omega_{j}}{2}\ket{g_{i}}\bra{r_{j}}+h.c.\right).$

Here $\hat{H}_{\mathrm{vdw}}(\mathbf{r})$ is the Hamiltonian that describes
the van der Waals interaction between two Rydberg atoms assuming that the
Rydberg states are far away from any Förster resonance and is given as

$\hat{H}_{\mathrm{vdw}}(\mathbf{r})=\sum_{\eta}\frac{C_{6}(\eta)}{r^{6}}\mathcal{\hat{D}}_{\eta}(\theta,\phi).$
(2)

The van der Waals interactions have two effects: a shift in Rydberg levels and
a mixing of different Rydberg states within the same Zeeman manifold Walker
and Saffman (2008). Depending on the choice of the Rydberg states
($\ket{r_{i=1,2}}$), there are different channels of intermediate states that
couple to $\ket{r_{i=1,2}}$ which are denoted by $\eta$. The overall strength
of the Rydberg interaction is determined by summing over all the relevant
channels as defined in Eq. 2. For each channel, there is a van der Waals
coefficient $C_{6}(\eta)$ which depends on the radial component of the dipole
matrix element, while $\mathcal{\hat{D}}_{\eta}(\theta,\phi)$ is the operator
that projects onto the subspace of magnetic quantum numbers allowed by the
selection rules. Thus $\mathcal{\hat{D}}_{\eta}(\theta,\phi)$ has the
information about the anisotropy in the Rydberg interactions. Explicit
expressions for $C_{6}(\eta)$ and $\mathcal{\hat{D}}_{\eta}(\theta,\phi)$ are
provided in Sec. 1 of the Appendix.

Rydberg-dressing.— Although bare ground state atoms do not interact, effective
interactions between ground state atoms can be induced as a result of Rydberg-
dressing Glaetzle  _et al._ (2015). Rydberg dressing is a consequence of weak
optical coupling of the ground states to their corresponding Rydberg states,
which is satisfied when $\Delta_{1,2}\gg\Omega_{1,2}$ Johnson and Rolston
(2010); Honer  _et al._ (2010); Balewski  _et al._ (2014); Mukherjee  _et al._
(2016). In this limit, perturbation theory is used to determine
$\tilde{V}^{ij}_{lk}(\mathbf{r})$ which denote the effective interactions
between a pair of Rydberg-dressed atoms with ground states
$\ket{\tilde{g}_{i,j=1,2}}$ and $\ket{\tilde{g}_{l,k=1,2}}$ that are separated
by a distance vector $\mathbf{r}$ Glaetzle  _et al._ (2015). The final
expressions $\tilde{V}^{ij}_{lk}$ depend on the bare Rydberg interactions as
well as the laser parameters $(\Omega,\Delta)$ and have the following form,

$H_{\text{eff}}=\sum_{\begin{subarray}{c}i,j\\\
l,k\end{subarray}}\tilde{V}^{ij}_{lk}(\mathbf{r})\ket{\tilde{g}_{i}\tilde{g}_{j}}\bra{\tilde{g}_{l}\tilde{g}_{k}},$
(3)

where

$\displaystyle\tilde{V}^{ij}_{lk}(\mathbf{r})$
$\displaystyle=U^{ij(2)}_{lk}(\mathbf{r})+U^{ij(4)}_{lk}(\mathbf{r}),$ (4)
$\displaystyle U^{ij(2)}_{lk}(\mathbf{r})$
$\displaystyle=\sum_{a}\frac{\bra{g_{i}g_{j}}\hat{H}^{\prime}\ket{a}\bra{a}\hat{H}^{\prime}\ket{g_{l}g_{k}}}{E_{g}-E_{a}},$
(5) $\displaystyle U^{ij(4)}_{lk}(\mathbf{r})$
$\displaystyle=\sum_{abc}\frac{\langle
g_{i}g_{j}|\hat{H}^{\prime}|a\rangle\langle a|\hat{H}^{\prime}|b\rangle\langle
b|\hat{H}^{\prime}|c\rangle\langle
c|\hat{H}^{\prime}|g_{l}g_{k}\rangle}{\left(E_{g}-E_{a}\right)\left(E_{g}-E_{b}\right)\left(E_{g}-E_{c}\right)}$
$\displaystyle-U^{ij(2)}_{lk}\sum_{a}\frac{\bra{g_{i}g_{j}}\hat{H}^{\prime}\ket{a}\bra{a}\hat{H}^{\prime}\ket{g_{l}g_{k}}}{\left(E_{g}-E_{a}\right)^{2}}.$
(6)

The effective interaction between the dressed ground states has a non-zero
contribution from the even order terms in the perturbation theory indicated by
$U^{ij(2,4)}_{lk}$. Higher than fourth-order processes are negligible for our
purposes. The different channels that contribute to the second and fourth-
order processes have been illustrated in Fig. 2(a,b). The two-atom states with
single Rydberg excitation are labeled as
$\ket{a},\ket{c}=\\{\ket{r_{1}g_{2}},\ket{r_{2}g_{1}},\ket{g_{1}r_{2}},\ket{g_{2}r_{1}}\\}$
with energy $E_{a,c}$ and $\ket{b}=\ket{E_{i=0,1,2,3}}$ correspond to both
atoms excited to Rydberg states with energy $E_{b}=E_{i=0,1,2,3}$. The energy
eigenstates $\ket{E_{i=0,1,2,3}}$ are obtained by diagonalizing $\hat{H}_{0}$,
which already take into account the Zeeman mixing of relevant states. The
strength of Rydberg-dressed interactions is controlled by laser parameters,
which are encapsulated in Eq. 5-6. The most generalized form for $\tilde{V}$
and $U$ for arbitrary $(\theta,\phi)$ is complex but can be evaluated
numerically. Nevertheless, for a specific case of $(\theta=0,\phi=0)$, one
obtains relatively simple analytical expressions, which are provided in Sec. 2
of the Appendix.

Mapping to spin models.— The two-atom interaction picture in the previous
section is now generalized to the spin model. The Rydberg-dressed ground
states $\ket{\tilde{g}_{1}}=\ket{\downarrow}$ and
$\ket{\tilde{g}_{2}}=\ket{\uparrow}$ have small admixtures of Rydberg states
into the bare ground states and are used to define the spin-1/2 system. Using
the standard definitions of the Pauli spin operators for
$\hat{\sigma}^{z}=1/2(\ket{\uparrow}\bra{\uparrow}-\ket{\downarrow}\bra{\downarrow})$
and
$\hat{\sigma}^{x}=1/2(\ket{\uparrow}\bra{\downarrow}+\ket{\downarrow}\bra{\uparrow})$,
the one-dimensional TFIM for $N$ spins is written as

$\displaystyle\hat{H}=$
$\displaystyle\sum_{p<q}\left[J_{z}(\mathbf{r}_{pq})\hat{\sigma}^{z}_{p}\hat{\sigma}^{z}_{q}+h_{x}(\mathbf{r}_{pq})\hat{\sigma}^{x}_{q}(\mathbf{r}_{pq})+h_{z}(\mathbf{r}_{pq})\hat{\sigma}^{z}_{q}\right],$
(7) $\displaystyle J_{z}(\mathbf{r})=$
$\displaystyle\frac{1}{4}\left(\tilde{V}^{11}_{11}(\mathbf{r})-\tilde{V}^{12}_{12}(\mathbf{r})-\tilde{V}^{21}_{21}(\mathbf{r})+\tilde{V}^{22}_{22}(\mathbf{r})\right),$
(8) $\displaystyle h_{x}(\mathbf{r})=$
$\displaystyle\frac{1}{2}\left(\tilde{V}^{12}_{11}(\mathbf{r})+\tilde{V}^{21}_{11}(\mathbf{r})+\tilde{V}^{11}_{12}(\mathbf{r})+\tilde{V}^{11}_{21}(\mathbf{r})\right.$
$\displaystyle\left.+\tilde{V}^{22}_{12}(\mathbf{r})+\tilde{V}^{22}_{21}(\mathbf{r})+\tilde{V}^{12}_{22}(\mathbf{r})+\tilde{V}^{21}_{22}(\mathbf{r})\right),$
(9) $\displaystyle h_{z}(\mathbf{r})=$
$\displaystyle\frac{1}{4}\left(\tilde{V}^{22}_{22}(\mathbf{r})-\tilde{V}^{11}_{11}(\mathbf{r})\right),$
(10)

with $p,q$ labelling the lattice sites with position vectors
$\mathbf{r}_{p},\mathbf{r}_{q}$ respectively and
$\mathbf{r}_{pq}=\mathbf{r}_{p}-\mathbf{r}_{q}$ is the distance vector between
the sites. The $h_{x/z}$ are the transverse/longitudinal field strengths,
respectively, which in the Rydberg-dressed picture also depend on the distance
vector as shown for $h_{x}$. The term $J_{z}(\mathbf{r})=J_{z}(r,\theta,\phi)$
determines the strength of the Ising interaction between spins with distance
$r$ away from each other and depends on the choice of the principal quantum
number for the Rydberg state. The angular dependence of $J_{z}$, which is
responsible for the non-monotonic nature of the interactions, can be explored
by varying the angle between the laser propagation direction and the relative
vector connecting the two spins Reinhard  _et al._ (2007); Glaetzle  _et al._
(2015). Typically for the same quantum number, the effective
transverse/longitudinal field strengths can be orders of magnitude smaller
than the effective Ising interactions. The parameters chosen in this work for
the mesonic dynamics are such that the fields are kept constant over distance,
and only $J_{z}$ varies with distance.

In the semi-classical limit for which $h_{x}\rightarrow 0$, one obtains
fermionic excitations that are well described by free domain walls. By
increasing the transverse field and setting $h_{z}=0$, these domain walls gain
kinetic energy. Also, for a nonzero transverse field, the number of domain
walls is only quasi-conserved, but for small enough transverse fields
$h_{x}\ll J$, one can project the dynamics to subspaces with a conserved
number of domain walls Rutkevich (2008). When including the long-range
character in the Ising interactions, one obtains confinement potentials for
otherwise non-interacting free domain walls Liu  _et al._ (2019). For this
work, we choose parameters such that the effective Ising interaction is
attractive. One issue that arises in the context of Rydberg dressing with van
der Waals type of interaction is that the asymptotic decay is too short-ranged
($1/r^{6}$ scaling) for the individual mesons to be stable Kormos  _et al._
(2017). This can be resolved by either exploiting the long-ranged dipole-
dipole interactions of the Rydberg states Letscher and Petrosyan (2018) or by
switching on the longitudinal field, $h^{z}\neq 0$ which also ensures the
stable formation of individual meson states Kormos  _et al._ (2017); Vovrosh
and Knolle (2021). The latter scheme is used in this work. In either case, the
short-distance dynamics of the two mesons can still be governed by the non-
monotonic effective interactions obtained from Rydberg dressing.

Figure 3: The figure shows the dependence of the transverse field as a
function of distance on various parameters. (a) Varying principal quantum
number $n$, with fixed
$(\Omega_{1},\Omega_{2},\Delta_{1},\Delta_{2})=(5,-10,30,-70)~{}2\pi$ MHz,
$(\theta,\phi)=(\pi/3,\pi/9)$ (b) Varying the scaling factor $\gamma$, with
fixed $n=100$,
$(\Omega_{1},\Omega_{2},\Delta_{1},\Delta_{2})=(5,-10,30,-70)\gamma~{}2\pi$
MHz, $(\theta,\phi)=(\pi/3,\pi/9)$ (c) Varying $\theta$, with fixed $n=100$,
$(\Omega_{1},\Omega_{2},\Delta_{1},\Delta_{2})=(5,-10,30,-70)~{}2\pi$ MHz,
$\phi=\pi/9$.

Mapping to domain-wall physics.— The mapping of the Ising model to the mesonic
Hamiltonian is well known Kormos  _et al._ (2017); Liu  _et al._ (2019) and
more details are provided in Sec. 3-4 of the Appendix. We focus on the
simplest case of a pair of 1-mesons (mesons of width $1$ site) located at
sites $j_{1}$ and $j_{2}$ which is denoted by
$\ket{\uparrow\ldots\uparrow\downarrow_{j_{1}}\uparrow\ldots\uparrow\downarrow_{j_{2}}\uparrow\ldots\uparrow}$.
Using perturbation theory for weak $h^{x}$, the Ising model given in Eq. 7 is
mapped to an effective Hamiltonian describing the interaction between the
1-mesons in the reduced basis, which is given as Vovrosh  _et al._ (2022)

$\displaystyle H_{\text{mesons}}=$
$\displaystyle\sum_{j_{1},j_{2}}h_{j}[\ket{j_{1},j_{2}+1}+\ket{j_{1}-1,j_{2}}]\bra{j_{1},j_{2}}+h.c.$
$\displaystyle+\big{(}U_{j}-\frac{4}{j^{\alpha}}\big{)}\ket{j_{1},j_{2}}\bra{j_{1},j_{2}},$
(11)

where $j=j_{2}-j_{1}$, $h_{j}$ is the inhomogeneous hopping term and $U_{j}$
is the effective interaction. The hadron-like bound state of two mesons is not
coupled to the asymptotic states of free mesons and can have a significant
energy gap, especially for small meson sizes. Nevertheless, as already
mentioned in the introduction, the transition from free mesons to bound-state
mesons can be enhanced either by abruptly changing the transverse field or by
modifying the long-range interactions. Both these schemes can be applied to
Rydberg setups, although the latter scheme is more challenging to obtain in
other cold-atom quantum simulators and is the focus of this work. Eq. Quantum
simulation of hadronic states with Rydberg-dressed atoms is used for the
numerical dynamics shown in Fig. 1(b).

Figure 4: The figure shows the spatial variation of the resonant-type
effective spin interaction $J^{z}(\mathbf{r})$ for different parameters.
Angular variation of $J^{z}(\mathbf{r})$ with respect to $\theta$ for
$\phi=0$, $n=60$ and
$(\Omega_{1},\Omega_{2},\Delta_{1},\Delta_{2})=(5,-10,30,-70)~{}2\pi$ MHz in
(a) and with respect to $\phi$ for $\theta=\pi/2$ in (b). Radial variation of
$J^{z}(\mathbf{r})$ for different principal quantum number $n$ for
$(\theta,\phi)=(\pi/4,\pi/2)$ and
$(\Omega_{1},\Omega_{2},\Delta_{1},\Delta_{2})=(5,-10,30,-70)~{}2\pi$ MHz in
(c) and scaling factor $\gamma$ in
$(\Omega_{1},\Omega_{2},\Delta_{1},\Delta_{2})=(5,-10,30,-70)\gamma~{}2\pi$
MHz for $n=100$ and $(\theta,\phi)=(\pi/4,\pi/2)$ in (d).

Results.— Motivated by prior works on quantum simulation of spin models using
Rydberg-dressing Glaetzle  _et al._ (2015); van Bijnen and Pohl (2015), we
engineer effective spin interactions for the TFIM that are tailored to achieve
successful hadron formation based on Vovrosh  _et al._ (2022). The key
requirement is to have an effective Ising interaction that has a local minimum
with a decaying tail in distance. The strength and shape of the effective
Ising interaction are controlled by the laser parameters $(\Omega,\Delta)$,
the direction of the excitation laser with respect to the inter-nuclear
distance, the lattice spacing, and the principal quantum $n$. In this article,
$\text{Rb}^{87}$ atoms were chosen with the two hyperfine sublevels of
$\ket{5^{2}S_{1/2}}$ as the ground states and
$\ket{r_{1}}=\ket{60~{}^{2}P_{1/2},m_{j}=-1/2}$,
$\ket{r_{2}}=\ket{60~{}^{2}P_{1/2},m_{j}=1/2}$ as the Rydberg states. The van
der Waals interactions between Rydberg states for Rb atoms were calculated
using a Python package Alkali-Rydberg Calculator (ARC)Robertson  _et al._
(2021). The anisotropic nature of the dressed potentials, which stems from
using non-spherical ($nP$) Rydberg states, is scanned by varying the direction
of laser propagation with respect to the lattice axis. The energy scales of
the detuning and the van der Waals interactions are considered to be large
enough such that the hyperfine interactions in the Rydberg state can be
ignored.

The dependence of the transverse field on distance for chosen set of
parameters is shown in Fig. 3. The typical values of the transverse field is
few orders of magnitude smaller than that of the Ising interaction for the
same quantum number. Similarly we show the three different types of effective
Ising potentials obtained via the Rydberg-dressing of Rb atoms: “Resonance-
type”, “Step-type” and the “Hump-type” in Figs. 4 \- 6. In these figures, one
can see the how the non-monotonic nature of the potentials can be controlled
using different parameters.

Figure 5: The figure shows the spatial variation of the step-type effective
spin interaction $J^{z}(\mathbf{r})$ for different parameters. Angular
variation of $J^{z}(\mathbf{r})$ with respect to $\theta$ for $\phi=0$, $n=60$
and $(\Omega_{1},\Omega_{2},\Delta_{1},\Delta_{2})=(0.1,10,10,-100)~{}2\pi$
MHz in (a) and with respect to $\phi$ for $\theta=\pi/2$ in (b). Radial
variation of $J^{z}(\mathbf{r})$ for different principal quantum number $n$
for $(\theta,\phi)=(\pi/2,0)$ and
$(\Omega_{1},\Omega_{2},\Delta_{1},\Delta_{2})=(0.1,10,10,-100)~{}2\pi$ MHz in
(c) and scaling factor $\gamma$ in
$(\Omega_{1},\Omega_{2},\Delta_{1},\Delta_{2})=(0.01,1,1,-10)\gamma~{}2\pi$
MHz for $n=100$ and $(\theta,\phi)=(\pi/2,0)$ in (d).

Figures 4 \- 6(a,b) explicitly depict the much-needed spatial anisotropy in
the effective Ising interaction $J^{z}(\mathbf{r})$. Cases of how the
different parameters affect the shape of the interaction are shown in Fig. 4
\- 6 (c,d). The position of the resonance, the range of the step function, and
the peak of the hump all of which can be varied by simply varying the
principal quantum number, while the width of the resonance, the height of the
step function, or the hump are controlled by the laser properties. The origin
of these shapes can be understood from Eq. (6). The energies $E_{a}$ and
$E_{c}$ are fixed in inter-atomic distance and correspond to the energies of
the singly excited state in the rotating frame. The distance dependence in
energy exists only for $E_{b}(r)$ which, as explained before, arises from the
Rydberg induced Zeeman mixing. For example, the resonance-type potentials in
the dressed interaction picture as shown in Fig. 4 occur whenever $E_{b}(r)$
has a zero crossing in the rotating frame. Since $E_{b}(r)$ is a non-monotonic
function of distance $r$, it is possible to have multiple zero crossings which
result in more than one resonance. More non-trivial Rydberg-dressed potential
shapes such as the hump and step (also referred to as soft-core potentials)
are also a result of the non-monotonic nature of $E_{b}(r)$ and are well
studied in other works Henkel  _et al._ (2010); Glaetzle  _et al._ (2015); van
Bijnen and Pohl (2015).

Figure 6: The figure shows the spatial variation of the hump-type effective
spin interaction $J^{z}(\mathbf{r})$ for different parameters. Angular
variation of $J^{z}(\mathbf{r})$ with respect to $\theta$ for $\phi=0$, $n=60$
and $(\Omega_{1},\Omega_{2},\Delta_{1},\Delta_{2})=(10,10,50,-50)~{}2\pi$ MHz
in (a) and with respect to $\phi$ for $\theta=\pi/2$ in (b). Radial variation
of $J^{z}(\mathbf{r})$ for different principal quantum number $n$ for
$(\theta,\phi)=(\pi/2,0)$ and
$(\Omega_{1},\Omega_{2},\Delta_{1},\Delta_{2})=(10,10,50,-50)~{}2\pi$ MHz in
(c) and scaling factor $\gamma$ in
$(\Omega_{1},\Omega_{2},\Delta_{1},\Delta_{2})=(10,10,50,-50)\gamma~{}2\pi$
MHz for $n=100$ and $(\theta,\phi)=(\pi/2,0)$ in (d).

Both the hump-type potential and the double-resonance type effective
potentials are most suitable for hardron formation. However, we take the
latter as an example, as shown in Fig. 7(a), where a certain fixed lattice
spacing is chosen, which is depicted by yellow dots corresponding to the
location of the atomic sites. As a result of this choice, one obtains an Ising
interaction profile for the dressed atoms as shown in Fig. 7(b). The non-
vanishing interaction energy at short distances in Fig. 7(b) is needed to
provide confinement as was shown in earlier works Kormos  _et al._ (2017); Liu
_et al._ (2019) and we go further than these results by introducing the
additional barrier at $|r_{p}-r_{q}|\sim$ 4 in order to provide a more stable
trap for the tetra-quark after its formation. For future works, it would be
interesting to explore the collisions of magnon excitations in the presence of
humps (Fig. 6) where we anticipate to have very similar results although will
lack the analogy with confinement physics. Using the effective Rydberg-dressed
potential shown in Fig. 7(b), the dynamics of two $1$-mesons are performed,
which is depicted in Fig. 1(b). The individual mesons are initially placed far
apart from each other at site number 0 and 30, respectively, where they
experience the asymptotic tail of the potential. A non-zero longitudinal field
($|h^{z}/J^{z}|=1$) is also added in order to keep the individual mesons
stable. By switching on the transverse field to a value $|h^{x}/J^{z}|=0.3$,
the static mesons start moving towards each other. Numerically, the mesonic
dynamics is treated by defining a Gaussian wavepacket for each meson which
evolves under the Hamiltonian $H_{\text{mesons}}$, and the parameters are
provided in Sec. 5 of the Appendix. Thus, two meson wavepackets with opposite
momenta approach each other until they see the potential barrier of the
effective interaction, which is located at $|r_{p}-r_{q}|\sim 4$ in Fig. 7(b)
or at site 15 in Fig. 1(b). While a large portion of the individual meson
wavepacket gets reflected due to the potential barrier, some of it tunnels
through the barrier to form a tetraquark. As a result of the quantum
tunneling, the lifetime of the meta-stable bound state can be calculated
Vovrosh  _et al._ (2022) which is estimated to be $\tau\sim 384.6/J_{z}\sim
3-15$ ms for $J_{z}=4-20~{}2\pi$ kHz. The interaction-induced fusion of mesons
occurs at around $tJ_{z}\sim 100$ (see Fig. 1(b)) which is around 795 $\mu$s.
The bare Rydberg lifetimes for the states chosen in this work is around 100
$\mu$s which is further enhanced to 10-100 ms due to dressing. This clearly
exemplifies the merit of Rydberg-dressing, which not only allows the
flexibility to vary the width and location of the double-resonances but also
provides an enhanced lifetime due to the fact that the interactions occur
effectively between ground state atoms via the weak coupling to Rydberg
states.

Considerations of Experimental feasibility.— The timescales provided in the
previous section to observe coherent dynamics take into account only the
spontaneous emission of the Rydberg state. However, in a real experiment, we
have other processes that negatively impact the coherence of the dynamics such
as (i) black-body effects, (ii) interaction-induced dephasing that arises from
fluctuations in atomic positions, and (iii) the presence of accidental
resonances of near-degenerate pair states, especially at small inter-atomic
distances which may drastically limit the lifetime of the dressed states. Our
setup is envisaged to be in a cryogenic environment ($\sim$10-100 nK) that
would suppress the blackbody-radiation-induced transitions to other Rydberg
states as well as Doppler broadening similar to most of the current
experimental schemes for Rydberg-dressing Zeiher  _et al._ (2016); Steinert
_et al._ (2023); Hollerith  _et al._ (2022). At few $\mu$m inter-atomic
distances, the energy mismatch between two pairs of Rydberg states greatly
exceeds their coupling strength resulting in a van der Waals potential. Thus,
at such distances, the simplified description of the Rydberg-dressed
interactions shown in Figs. 4 \- 6 are verifiable to a very good
approximation. Naturally, this description becomes increasingly inaccurate for
smaller distances which are plagued with a dense region of resonances due to
multiple energy crossings from neighbouring pairs of Rydberg states. The
density of the spaghetti region experiences a notable decrease beyond a
certain distance threshold depending on the principal quantum number. For our
case, a distance of around 2 $\mu$m is at the borderline of where effects of
the deep spaghetti region qualitatively affect the shape of the potential. The
thorough numerical treatment of molecular interaction potentials at such small
distances was done in van Bijnen and Pohl (2015) and is beyond the scope of
our work. However, it was confirmed in van Bijnen and Pohl (2015) and later
seen in experiments Steinert  _et al._ (2023); Hollerith  _et al._ (2022),
that despite the “spaghetti region of interaction potentials at such small
distances”, it is possible to obtain potential shapes qualitatively of the
type Figs. 4 \- 6 provided we choose lattice spacings larger than 1-2 $\mu$m
for principal quantum numbers of about $n=35-70$ for Rb. Interestingly
Förster-resonant interactions can also be exploited to obtain a more favorable
ratio of interactions and decoherence rates which is experimentally
achievable. The dynamical behavior of the mesons presented in this work should
be observable with current techniques of state-selective imaging of atoms
using quantum gas microscopes Schauß  _et al._ (2012); Schauß (2018); Browaeys
and Lahaye (2020).

Figure 7: (a) Rydberg-dressed resonance-type potential chosen for the meson
dynamics with parameters $n=70$, $(\theta,\phi)=(0,0)$,
$(\Omega_{1},\Omega_{2},\Delta_{1},\Delta_{2})=(5,10,47,-53)~{}2\pi$ MHz. (b)
By choosing a lattice spacing $a=2.1~{}\mu$m shown as yellow dots in (a), the
effective normalized Ising interaction is obtained for which a tetraquark
potential well is produced around site three and is used to obtain the
dynamics in Fig. 1(b). Here $|r_{p}-r_{q}|$ denote the difference in lattice
sites between the pair of mesons.

In these setups, it is natural to have the many-body state correspond to the
case where all spins are down. The initial state preparation requires a
localized spin flip at specific sites to create the four-kink state shown in
Fig. 1 at $t=0$. This can be achieved using a $\pi$ pulse laser whose width is
around few hundred nm. In general, to create a meson of certain size will
depend on the lattice spacing and the width of the laser beam Labuhn  _et al._
(2014); Bernien  _et al._ (2017). As mentioned before, mesons of a certain
width with short-range interactions can be stabilized with small but finite
longitudinal field. In order to push the individual mesons in opposite
directions towards each other, one requires a controlled protocol of spin
exchange. For Rydberg dressed-ground states, this can be achieved using
protocols similar to Li  _et al._ (2022). Alternatively, by dressing to
excited states one can implement selective spin exchanges via resonant dipole-
dipole interaction or by flip-flop type of van der Waals interaction between
Rydberg states Schempp  _et al._ (2015); Yang  _et al._ (2019). The issue with
microwave dressing to Rydberg states is that one needs to have shorter
timescales for coherent dynamics compared to dressed-ground states. The ansatz
of describing the wavepacket with certain momenta as a Gaussian needs to be
verified experimentally and will depend on various controllable parameters of
the hopping term such as the duration and the strength of the microwave
dressing laser. Non-Gaussian wavepackets will give similar results depending
on the weight of the different velocities that contribute to the wavepacket
dynamics but for the benefit of having a better visualization of the formation
of hadrons, we chose Gaussian wavepackets that propagate effectively with a
single velocity. Moreover, unlike the parameters used in the theoretical
simulations, we believe that there is room to allow the collision between
mesons with a small uncertainty in momentum space. After the initial kick to
the mesons provided by the controlled nearest neighbour spin exchange, the
speed of the mesons is related to the transverse field strength $h_{x}$.
However, since our analysis is in the limit of $h_{x}\ll J$ which guarantees
the fixed number of domain wall approximation, we will have slow moving
mesons. Although the typical timescales required to observe collisions shown
in this work are long compared to what has been achieved in current
experiments, we believe that there is scope in optimizing the protocol
dynamics, either by minimising the initial meson separation or by maximising
$h_{x}$ while maintaining the domain wall approximation.

Conclusion and Discussion.— Quantum simulation of confinement physics and
other high-energy physics phenomena on table-top, controllable experiments is
an exciting and growing area of research Bañuls  _et al._ (2020); Surace  _et
al._ (2020); Yang  _et al._ (2020); Zhou  _et al._ (2022); Liu  _et al._
(2020). The formation of long-lived hadronic states by the fusion of
individual mesons during their dynamics requires sufficient tunability in the
spatial distribution of the interaction profile Vovrosh  _et al._ (2022).
Whilst there exist many ultra-cold platforms, such as trapped ions, dipolar
gases, and polar molecules, which have the ability to tune the strength of
their long-range interactions, it may be challenging to control the non-
monotonic interaction profiles in such setups. On the other hand, Rydberg-
dressed interactions have the unique flexibility to reproduce a variety of
non-monotonic potentials, some of which were identified in this work in the
context of the one-dimensional Ising model. The shape and strength of the
effective potential needed to facilitate the creation of hadronic states can
be modulated by laser parameters and the choice of the Rydberg state, all of
which are controlled in current experiments with sufficient precision. Apart
from inheriting the much-needed anisotropy from the Rydberg-Rydberg
interactions into the effective spin-spin interactions, the Rydberg-dressed
system has an enhanced lifetime compared to the bare Rydberg system making it
possible to observe the relevant coherent dynamics leading to hadron formation
potentially achievable within experimental timescales. For future work, it
would be interesting to see if it is possible to stabilize mesons with larger
widths, perhaps by dressing close to Förster resonances or with the use of
hump-like potentials. There is scope for optimizing the mesonic dynamics such
that the time of collision is reduced. Lastly, a more complete study of the
dynamics requires taking into account the dissipative and decoherent processes
which play a vital role in real experiments.

###### Acknowledgements.

We are indebted to the fruitful and inspiring discussions with Christian Gross
and Rick van Bijnen. JV acknowledges the Samsung Advanced Institute of
Technology Global Research Partner- ship and travel support via the Imperial-
TUM flagship part- nership. The research is part of the Munich Quantum Valley,
which is supported by the Bavarian state government with funds from the
Hightech Agenda Bayern Plus.

APPENDIX

##  1 Anisotropic Rydberg-Rydberg interaction between two atoms

A shorthand notation is used for the selected Rydberg states $\ket{nP_{\pm
1/2}m_{i};nP_{\pm 1/2}m_{j}}$ for a pair of atoms which is given as
$\ket{m_{i},m_{j}}$. As a result of the van der Waals interaction, the
original Rydberg states $\ket{m_{i},m_{j}}$ couple to a bunch of intermediate
states
$\ket{m_{\alpha}m_{\beta}}\equiv\ket{n_{\alpha},l_{\alpha},j_{\alpha},m_{\alpha};n_{\beta},l_{\beta},j_{\beta},m_{\beta}}$
as permitted by the selection rules. The dipole-dipole interaction operator is
defined as

$\hat{V}_{\mathrm{dd}}(\mathbf{r})=-\sqrt{\frac{24\pi}{5}}\frac{1}{r^{3}}\sum_{\mu,\nu}C_{\mu,\nu,\mu+\nu}^{1,1;2}Y_{2}^{\mu+\nu}(\vartheta,\varphi)^{*}\hat{d}_{\mu}^{1}\hat{d}_{\nu}^{2},$
(A1)

where $\hat{d}_{\mu,\nu\in\\{-1,0,1\\}}^{i}$ are the spherical basis
components of the dipole operator $\hat{d}^{i}$, $Y_{l}^{m}$ are the spherical
harmonic functions and $C_{m_{1}m_{2};M}^{j_{1},j_{2};J}$ as the Clebsch-
Gordan coefficients. The matrix element
$\bra{m_{i}m_{j}}\hat{V}_{\mathrm{dd}}\ket{m_{\alpha}m_{\beta}}$ can be
expanded into radial and angular parts as follows,

$\displaystyle\bra{m_{i}m_{j}}\hat{V}_{\mathrm{dd}}\ket{m_{\alpha}m_{\beta}}=(\bra{n,1,1/2,m_{i}}\otimes\bra{n,1,1/2,m_{j}})\hat{V}_{\mathrm{dd}}(\ket{n_{\alpha},l_{\alpha},j_{\alpha},m_{\alpha}}\otimes\ket{n_{\beta},l_{\beta},j_{\beta},m_{\beta}})$
(A2)
$\displaystyle=-\sqrt{\frac{24\pi}{5}}\frac{1}{r^{3}}\sum_{\mu,\nu}C_{\mu,\nu,\mu+\nu}^{1,1;2}Y_{2}^{\mu+\nu}(\vartheta,\varphi)^{*}\bra{n,1,1/2,m_{i}}d_{\mu}^{1}\ket{n_{\alpha},l_{\alpha},j_{\alpha},m_{\alpha}}\bra{n,1,1/2,m_{j}}d_{\nu}^{2}\ket{n_{\beta},l_{\beta},j_{\beta},m_{\beta}}$
$\displaystyle=(-1)^{s-m_{i}}\sqrt{\left(2\ell_{i}+1\right)\left(2j_{i}+1\right)\left(2\ell_{\alpha}+1\right)\left(2j_{\alpha}+1\right)}\left\\{\begin{array}[]{ccc}\ell_{i}&\ell_{\alpha}&1\\\
j_{\alpha}&j_{i}&s\end{array}\right\\}\left(\begin{array}[]{ccc}\ell_{\alpha}&1&\ell_{i}\\\
0&0&0\end{array}\right)$
$\displaystyle\times(-1)^{s-m_{j}}\sqrt{\left(2\ell_{j}+1\right)\left(2j_{j}+1\right)\left(2\ell_{\beta}+1\right)\left(2j_{\beta}+1\right)}\left\\{\begin{array}[]{ccc}\ell_{j}&\ell_{\beta}&1\\\
j_{\beta}&j_{j}&s\end{array}\right\\}\left(\begin{array}[]{ccc}\ell_{\beta}&1&\ell_{j}\\\
0&0&0\end{array}\right)$
$\displaystyle\times\left[-\sqrt{\frac{24\pi}{5}}\sum_{\mu,\nu}C_{\mu,\nu;\mu+\nu}^{1,1;2}\left(\begin{array}[]{ccc}j_{\alpha}&1&j_{i}\\\
m_{\alpha}&\mu&-m_{i}\end{array}\right)\left(\begin{array}[]{ccc}j_{\beta}&1&j_{j}\\\
m_{\beta}&v&-m_{j}\end{array}\right)Y_{2}^{\mu+\nu}(\vartheta,\varphi)^{*}\right]\frac{\mathcal{R}^{\alpha}_{i}\mathcal{R}^{\beta}_{j}}{r^{3}}$
$\displaystyle=\mathcal{M}^{ij}_{\alpha\beta}\frac{\mathcal{R}^{\alpha}_{i}\mathcal{R}^{\beta}_{j}}{r^{3}}.$

In the third line, the angular part is expanded using the Wigner-Eckart
theorem and is denoted by $\mathcal{M}^{ij}_{\alpha\beta}$ in the last line.
The radial part is expressed as $\mathcal{R}_{\alpha}^{i}=\int
r^{2}dr\psi_{n_{i},\ell_{i},j_{i}}(r)^{*}r\psi_{n_{\alpha},\ell_{\alpha},j_{\alpha}}(r)$.
Using the above expressions, the interaction Hamiltonian is written as

$\displaystyle\hat{H}_{\mathrm{vdw}}(\mathbf{r})=$
$\displaystyle\sum_{\alpha,\beta}\frac{\bra{m_{i}m_{j}}\hat{V}_{\mathrm{dd}}\ket{m_{\alpha}m_{\beta}}\bra{m_{\alpha}m_{\beta}}\hat{V}_{\mathrm{dd}}\ket{m_{i^{\prime}}m_{j^{\prime}}}}{\delta_{\alpha\beta}}\ket{m_{i}m_{j}}\bra{m_{i^{\prime}}m_{j^{\prime}}}$
$\displaystyle=$
$\displaystyle\sum_{\begin{subarray}{c}\eta=l_{\alpha},l_{\beta}\\\
l_{\beta},j_{\beta}\end{subarray}}\sum_{n_{\alpha},n_{\beta}}\frac{\mathcal{R}_{i}^{\alpha}\mathcal{R}_{j}^{\beta}\mathcal{R}_{i^{\prime}}^{\alpha}\mathcal{R}_{j^{\prime}}^{\beta}}{\delta_{\alpha\beta}r^{6}}\left[\sum_{m_{\alpha},m_{\beta}}\mathcal{M}^{ij}_{\alpha\beta}\mathcal{M}_{i^{\prime}j^{\prime}}^{\alpha\beta}\right]\ket{m_{i}m_{j}}\bra{m_{i^{\prime}}m_{j^{\prime}}}$
$\displaystyle=$
$\displaystyle\sum_{\eta}\frac{C_{6}(\eta)}{r^{6}}\mathcal{\hat{D}}_{\eta}(\theta,\phi)$
(A3)

where $\delta_{\alpha\beta}=E_{i}+E_{j}-(E_{\alpha}+E_{\beta})$ is the energy
difference between the relevant states,
$C_{6}(\eta)=\sum_{n_{\alpha},n_{\beta}}(\mathcal{R}_{i}^{\alpha}\mathcal{R}_{j}^{\beta}\mathcal{R}_{i^{\prime}}^{\alpha}\mathcal{R}_{j^{\prime}}^{\beta})/\delta_{\alpha\beta}$
is the van der Waals coefficient which determines the strength of the
interactions and
$\mathcal{\hat{D}}_{\eta}(\theta,\phi)=\sum_{m_{\alpha},m_{\beta}}\mathcal{M}^{ij}_{\alpha\beta}\mathcal{M}_{i^{\prime}j^{\prime}}^{\alpha\beta}\ket{m_{i}m_{j}}\bra{m_{i^{\prime}}m_{j^{\prime}}}$
is the angular operator which describes the anisotropy in the interactions.
The different channels have been collectively denoted by $\eta$. Based on the
properties of the Wigner 3j and 6j symbols, we can deduce the selection rules,
which determine that the matrix elements will be non-zero provided the
conditions $m_{i}+\mu=m_{\alpha}$ and $m_{j}+\nu=m_{\beta}$ are satisfied.
Having chosen $\ket{P_{\pm 1/2};P_{\pm 1/2}}$ as our Rydberg states, all the
possible inter-mediate states that define the channels are (i)
$\ket{S_{1/2};S_{1/2}}$, (ii) $\ket{D_{3/2};D_{3/2}}$
(iii)$\ket{S_{1/2};D_{3/2}}$ and (iv) $\ket{D_{3/2};S_{1/2}}$.

##  2 Effective interactions between Rydberg-dressed ground states for a pair
atoms $(\theta,\phi)=(0,0)$

The unperturbed Hamiltonian for a pair of atoms for $\theta=0,\phi=0$ is given
as

$\hat{H}_{0}(r)=\left(\begin{array}[]{cccc}V_{22}^{22}(r)-2\Delta_{1}&0&0&0\\\
0&V_{21}^{21}(r)-\Delta_{1}-\Delta_{2}&V_{12}^{21}(r)&0\\\
0&V_{21}^{12}(r)&V_{12}^{12}(r)-\Delta_{1}-\Delta_{2}&0\\\
0&0&0&V_{11}^{11}(r)-2\Delta_{2}\end{array}\right),$ (A4)

where
$V^{ij}_{lk}(r)=\bra{r_{i}r_{j}}\hat{H}_{\mathrm{vdw}}(r)\ket{r_{l}r_{k}},\\{i,j,l,k\\}\in\\{1,2\\}$
and are explicitly given as

$\displaystyle
V^{11}_{11}(r)=V^{22}_{22}(r)=\frac{2}{81}(2C^{(i)}_{6}+11C^{(ii)}_{6}+14C^{(iii)}_{6})/r^{6}$
(A5) $\displaystyle
V^{12}_{12}(r)=V^{21}_{21}(r)=\frac{2}{81}(4C^{(i)}_{6}+13C^{(ii)}_{6}+10C^{(iii)}_{6})/r^{6}$
$\displaystyle
V^{12}_{21}(r)=V^{21}_{12}(r)=\frac{8}{81}(C^{(i)}_{6}+C^{(ii)}_{6}-2C^{(iii)}_{6})/r^{6}.$

The $C^{(i,ii,iii)}_{6}$ are the dispersion coefficients based on the
different channels described in the previous section. Diagonalizing the
Hamiltonian $\hat{H}_{0}$ gives the following eigenstates

$\displaystyle\ket{E_{1}}=\ket{r_{2}r_{2}},~{}\ket{E_{2}}=\ket{r_{1}r_{1}}~{}~{}\text{with
eigenvalues}~{}~{}E_{1}(r)=V_{11}^{11}(r),~{}E_{2}(r)=V_{11}^{11}(r),$ (A6)
$\displaystyle\ket{E_{3}}=\frac{1}{\sqrt{2}}\left(\ket{r_{1}r_{2}}+\ket{r_{2}r_{1}}\right)~{}~{}\text{with
eigenvalue}~{}~{}E_{3}(r)=V_{12}^{12}(r)+V_{12}^{21}(r)-\Delta_{1}-\Delta_{2},$
$\displaystyle\ket{E_{4}}=\frac{1}{\sqrt{2}}\left(\ket{r_{1}r_{2}}-\ket{r_{2}r_{1}}\right)~{}~{}\text{with
eigenvalue}~{}~{}E_{4}(r)=V_{12}^{12}(r)-V_{12}^{21}(r)-\Delta_{1}-\Delta_{2}.$

As described in the main text, the effective interaction matrix between the
two-atom ground states as a result of Rydberg dressing is derived by treating
the laser field as a perturbation up to its fourth-order term. In the two atom
picture, using the Rydberg-dressed atoms
$\ket{\tilde{g}_{i=1,2}\tilde{g}_{j=1,2}}$ as the basis, the effective
interactions take the form

$\displaystyle\hat{H}_{eff}(r)$
$\displaystyle=\left(\begin{array}[]{cccc}\tilde{V}^{22}_{22}(r)&0&0&0\\\
0&\tilde{V}^{21}_{21}(r)&\tilde{V}^{21}_{12}(r)&0\\\
0&\tilde{V}^{12}_{21}(r)&\tilde{V}^{12}_{12}(r)&0\\\
0&0&0&\tilde{V}^{11}_{11}(r)\end{array}\right),$ (A11) where
$\displaystyle~{}~{}~{}\tilde{V}^{ij}_{lk}(r)=U^{ij(2)}_{lk}+U^{ij(4a)}_{lk}(r)+U^{ij(4b)}_{lk},~{}~{}~{}\\{i,j,l,k\\}\in\\{1,2\\}.$
(A12)

In the above Hamiltonian, $\tilde{V}^{ij}_{lk}=\tilde{V}^{lk}_{ij}$ due to the
symmetry of the Hamiltonian $H$. The matrix elements are given as

$\displaystyle
U_{11}^{11(2)}=\frac{\Omega_{2}^{2}}{2\Delta_{2}},~{}U_{22}^{22(2)}=\frac{\Omega_{1}^{2}}{2\Delta_{1}},~{}U_{21}^{21(2)}=\frac{\Omega_{1}^{2}}{4\Delta_{1}}+\frac{\Omega_{2}^{2}}{4\Delta_{2}},$
(A13) $\displaystyle
U_{11}^{11(4a)}(r)=-\frac{\Omega_{2}^{4}}{4\Delta_{2}^{2}}\frac{V^{11}_{11}(r)-2\Delta_{1}}{\left(V^{11}_{11}(r)-2\Delta_{1}\right)\left(V^{11}_{11}(r)-2\Delta_{2}\right)},~{}U_{22}^{22(4a)}(r)=-\frac{\Omega_{1}^{4}}{4\Delta_{1}^{2}}\frac{V^{11}_{11}(r)-2\Delta_{2}}{\left(V^{11}_{11}(r)-2\Delta_{1}\right)\left(V^{11}_{11}(r)-2\Delta_{2}\right)}$
(A14) $\displaystyle
U_{12}^{21(4a)}(r)=\frac{\Omega_{1}^{2}\Omega_{2}^{2}}{16\Delta_{1}^{2}\Delta_{2}^{2}}\frac{\left(\Delta_{1}+\Delta_{2}\right)^{2}V^{12}_{21}(r)}{\left(\Delta_{1}+\Delta_{2}-V^{12}_{21}(r)\right)^{2}-(V^{12}_{21}(r))^{2}},~{}U_{21}^{21(4a)}(r)=\frac{\Omega_{1}^{2}\Omega_{2}^{2}}{16\Delta_{1}^{2}\Delta_{2}^{2}}\frac{\left(\Delta_{1}+\Delta_{2}\right)^{2}\left(\Delta_{1}+\Delta_{2}-V_{12}^{12}(r)\right)}{\left(\Delta_{1}+\Delta_{2}-V_{12}^{12}(r)\right)^{2}-\left(V_{21}^{12}(r)\right)^{2}},$
(A15) $\displaystyle
U_{22}^{22(4b)}=-\frac{\Omega_{1}^{4}}{4\Delta_{1}^{3}},~{}U_{11}^{11(4b)}=-\frac{\Omega_{2}^{4}}{4\Delta_{2}^{3}},~{}U_{21}^{21(4b)}=-\frac{\Omega_{1}^{4}}{16\Delta_{1}^{3}}-\frac{\Omega_{2}^{2}\Omega_{1}^{2}}{16\Delta_{1}^{2}\Delta_{2}}-\frac{\Omega_{2}^{2}\Omega_{1}^{2}}{16\Delta_{1}\Delta_{2}^{2}}-\frac{\Omega_{2}^{4}}{16\Delta_{2}^{3}}.$
(A16)

where $U_{21}^{21(2)}=U_{12}^{12(2)}$ and
$U_{21}^{21(4a,b)}(r)=U_{12}^{12(4a,b)}$. The most general case corresponding
to arbitrary angles will have all the off-diagonal terms in the matrix
$\hat{H}_{eff}$ to be non-zero.

##  3 Ising model in the basis of the two-domain wall states

This and the next section have been detailed in Vovrosh  _et al._ (2022),
which we briefly summarize here. The Hamiltonian for the one-dimensional Ising
spin chain is given as

$\hat{H}=\sum_{p<q}\left[J_{z}(\mathbf{r}_{pq})\sigma^{z}_{p}\sigma^{z}_{q}+h_{x}\sigma^{x}_{p}+h_{z}\sigma^{z}_{p}\right],$
(A17)

where all the terms are described in the main text of the article. $h_{x,z}$
have been chosen to be constant over distance. For $h_{x}=0$, the natural
excitations in the model are described by the non-interacting domain walls.
One can now consider the two domain-wall sectors with the basis
$\ket{j,n}=\ket{\uparrow...\uparrow\downarrow_{j}...\downarrow\uparrow_{j+n}...\uparrow}$.
This basis $|j,n\rangle$ labels the position of the first domain wall from the
left, and $n$ is the relative distance of the second domain wall from the
first. On this basis, the Ising model can be re-written as the following,

$H=\sum_{j,n}V(n)\ket{j,n}\bra{j,n}-h\big{[}\ket{j+1,n-1}\\\
+\ket{j-1,n+1}+\ket{j,n+1}+\ket{j,n-1}\big{]}\bra{j,n},$ (A18)

and can be understood as a kinetic term for each domain-wall with hopping
strength $h$ and a potential $V(n)$. For two domain-walls that are $n$ sites
apart feel an attractive potential $V(n)$ of the form Liu  _et al._ (2019)

$V(n)=4n\zeta(\alpha)-4\sum_{1\leq l<n}\sum_{1\leq r<l}\frac{1}{r^{\alpha}},$
(A19)

and $\zeta(\alpha)$ is the Riemann zeta function. The most interesting regime
is $1<\alpha<2$ where $V(n)$ is unbounded, and pulling the domain walls
infinitely far apart becomes energetically costly, leading to confinement.
Nonetheless, for larger $\alpha$, the two-kink subspace shows the presence of
deep bound states and asymptotic states of freely propagating domain-walls can
be very energetically high with very little overlap with the bound states.
However, in order to obtain stable mesons for $\alpha=6$, one requires a non-
zero longitudinal field $h_{z}$.

The study of hadron-like excitations requires considering interactions among
mesons. For this purpose, we focus on the two 1-meson subspaces by extending
the two domain-wall projections. For the $1-$meson states, we have $n=1$ and
can be dropped from the notation. For a pair of $1-$meson subspace, we use the
$|j_{1},j_{2}\rangle$ notation, which clearly indicates the position of the
two individual mesons. On this basis, the projected Hamiltonian takes the
shorthand form

$\displaystyle H=$
$\displaystyle\sum_{j_{1},j_{2}}\big{[}V(1)\ket{j_{1},j_{2}}\bra{j_{1},j_{2}}+V(2)\ket{j_{1},j_{2}}\bra{j_{1},j_{2}}+I(j_{1},j_{2})\ket{j_{1},j_{2}}\bra{j_{1},j_{2}}\big{]}$
$\displaystyle-h\big{[}\ket{j_{1}+1,j_{2}}+\ket{j_{1}-1,j_{2}}+\ket{j_{1},j_{2}+1}+\ket{j_{1},j_{2}-1}\big{]}\bra{j_{1},j_{2}}+h.c.$
(A20)

Here, $I$ can be seen as the meson interaction such that

$I(j_{1},j_{2})=-\frac{4}{(j_{1}-j_{2})^{\alpha}}$ (A21)

Note the form of the interaction, $I$, depends on the choice of boundary
conditions. Here we consider open boundaries. This effective interaction
between mesons falls off as $|j_{2}-j_{1}|^{-\alpha}$ and thus, similar to
domain-walls that are spatially distant from each other, individual mesons
which are far apart from each other interact only weakly.

##  4 Effective Hamiltonian for a pair of 1-mesons

The mesonic dynamics is constrained to the $1-$meson subspace, which ignores
any spin-flip processes (justified for a weak enough transverse field $h$) and
is well described with perturbation theory. The Ising model expressed using
the two-domain-wall basis states as shown in Eq. 3 is factored into a diagonal
part $H_{0}$ and an off-diagonal part $V$ proportional to what we will
consider the perturbative parameter $h$, that is $H=H_{0}+hV$. Using the
Schrieffer-Wolff transformation, a unitary transformation is generated by the
operator $S$ such that we have the expansion $H^{eff}=e^{S}He^{-S}\sim
H_{0}+h^{2}M+O(h^{3})$. The Baker-Campbell-Hausdorff expansion of
$e^{S}He^{-S}$ up to second order in $h$ gives

$H^{eff}=e^{S}He^{-S}=H_{0}+hV+[S,H_{0}]+h[S,V]+\frac{1}{2}[S,[S,V]]+O(h^{3})\
,$ (A22)

First-order terms are eliminated by choosing $S$ such that $[S,H_{0}]=-hV$.
The simplest way to do this is element-by-element, enforcing that

$\langle p|S|q\rangle=\frac{h\langle p|V|q\rangle}{E_{p}-E_{q}}\ ,$ (A23)

where $|p\rangle,|q\rangle$ are eigenvectors of $H_{0}$ with eigenvalues
$E_{p},E_{q}$. With this enforced and absorbing the higher order terms into
the $O(h^{3})$ factor, the effective Hamiltonian becomes

$H^{eff}=H_{0}+\frac{h}{2}[S,V]+O(h^{3})\ .$ (A24)

Thus, the resulting effective Hamiltonian contains only on-site potentials and
two-site hopping terms. The hopping coefficient $h_{n,n+1}$ determines the
moving of $1-$meson by one site. There is also a second order effect in which
a $1-$meson can ‘hop to itself’ through a process in which a neighboring spin
flips and then flips again to leave the initial state unchanged, $h_{n,n}$
which can be seen as an effective on-site potential. These contributions are
explicitly given by

$h_{j}=\frac{h^{2}}{2}\bigg{[}\frac{1}{V(1)-V(2)+\frac{4}{{j+1}^{\alpha}}}+\frac{1}{V(1)-V(2)+\frac{4}{{j}^{\alpha}}}\bigg{]},$
(A25)
$U_{j}=2h^{2}\bigg{[}\frac{1}{V(1)-V(2)+\frac{4}{{j+1}^{\alpha}}}+\frac{1}{V(1)-V(2)+\frac{4}{{j-1}^{\alpha}}}\bigg{]},$
(A26)

where $j=j_{2}-j_{1}$ and thus the effective Hamiltonian for a pair of
$1-$mesons can be written as shown in Eq. [9] in the main article.

##  5 Initial meson wavepacket preparation and dynamics

We theoretically choose Gaussian wavepackets of mesons as initial conditions
to allow for high control of meson momentum. In Rydberg atom setups, the many-
body product states corresponding to a pair of mesons are readily prepared
through local quenches induced by lasers pulses at the appropriate sites.
Here, the lattice spacing and the width of the laser beam determine the size
of the mesons realized Omran  _et al._ (2019); Labuhn  _et al._ (2014);
Bernien  _et al._ (2017). The momentum at which these mesons propagate along
the spin chain can be controlled by the transverse field. The initial states
used in the collision of mesons in Fig. 1(b) are given by

$\ket{\psi_{t=0}}\propto\sum_{j_{1},j_{2}}e^{-ikj_{1}-\frac{(j_{1}+n-2c_{1})^{2}}{4\sigma^{2}}}e^{ikj_{2}-\frac{(j_{2}+n-2c_{2})}{4\sigma^{2}}}\ket{j_{1},n,j_{2},n},$
(A27)

where $k$ is the initial momentum of each meson, $c_{i}$ is the center of the
wavepacket of the $i^{th}$ meson, and $\sigma$ is the standard deviation of
the wavepacket in real space. The parameters for the meson wavepackets in Fig.
1(b) are $q=1.5$, $\sigma_{k}=0.5$ and $\sigma_{c}=3$.

## References

  * Greensite (2011) J. Greensite,  _Lecture Notes in Physics: An introduction to the confinement problem_ , Vol. 821 (Springer-Verlag Berlin Heidelberg, 2011).
  * Bañuls  _et al._ (2020) M. C. Bañuls, R. Blatt, J. Catani, A. Celi, J. I. Cirac, M. Dalmonte, L. Fallani, K. Jansen, M. Lewenstein, S. Montangero, C. A. Muschik, B. Reznik, E. Rico, L. Tagliacozzo, K. Van Acoleyen, F. Verstraete, U.-J. Wiese, M. Wingate, J. Zakrzewski, and P. Zoller, The European Physical Journal D 74, 165 (2020).
  * Surace  _et al._ (2020) F. M. Surace, P. P. Mazza, G. Giudici, A. Lerose, A. Gambassi, and M. Dalmonte, Phys. Rev. X 10, 021041 (2020).
  * Yang  _et al._ (2020) B. Yang, H. Sun, R. Ott, H.-Y. Wang, T. V. Zache, J. C. Halimeh, Z.-S. Yuan, P. Hauke, and J.-W. Pan, Nature 587, 392 (2020).
  * Zhou  _et al._ (2022) Z.-Y. Zhou, G.-X. Su, J. C. Halimeh, R. Ott, H. Sun, P. Hauke, B. Yang, Z.-S. Yuan, J. Berges, and J.-W. Pan, Science 377, 311 (2022).
  * McCoy and Wu (1978) B. M. McCoy and T. T. Wu, Phys. Rev. D 18, 1259 (1978).
  * Delfino  _et al._ (1996) G. Delfino, G. Mussardo, and P. Simonetti, Nuclear Physics B 473, 469 (1996).
  * Kormos  _et al._ (2017) M. Kormos, M. Collura, G. Takács, and P. Calabrese, Nature Physics 13, 246 (2017).
  * Liu  _et al._ (2019) F. Liu, R. Lundgren, P. Titum, G. Pagano, J. Zhang, C. Monroe, and A. V. Gorshkov, Phys. Rev. Lett. 122, 150601 (2019).
  * Mazza  _et al._ (2019) P. P. Mazza, G. Perfetto, A. Lerose, M. Collura, and A. Gambassi, Phys. Rev. B 99, 180302 (2019).
  * Vovrosh and Knolle (2021) J. Vovrosh and J. Knolle, Scientific Reports 11, 11577 (2021).
  * Vovrosh  _et al._ (2021) J. Vovrosh, K. E. Khosla, S. Greenaway, C. Self, M. S. Kim, and J. Knolle, Phys. Rev. E 104, 035309 (2021).
  * Milsted  _et al._ (2022) A. Milsted, J. Liu, J. Preskill, and G. Vidal, PRX Quantum 3, 020316 (2022).
  * Karpov  _et al._ (2022) P. I. Karpov, G.-Y. Zhu, M. P. Heller, and M. Heyl, Phys. Rev. Res. 4, L032001 (2022).
  * Surace and Lerose (2021) F. M. Surace and A. Lerose, New Journal of Physics 23, 062001 (2021).
  * Vovrosh  _et al._ (2022) J. Vovrosh, R. Mukherjee, A. Bastianello, and J. Knolle, PRX Quantum 3, 040309 (2022).
  * Che  _et al._ (2015) J. Che, Z. Zhang, C. Li, Y. Tian, Y. Zhang, X. Wang, and Y. Zhang, Laser Physics Letters 12, 025403 (2015).
  * Zeiher  _et al._ (2016) J. Zeiher, R. V. Bijnen, P. Schauß, S. Hild, J. Y. Choi, T. Pohl, I. Bloch, and C. Gross, Nature Physics 12, 1095 (2016).
  * Balewski  _et al._ (2014) J. B. Balewski, A. T. Krupp, A. Gaj, S. Hofferberth, R. Löw, and T. Pfau, New Journal of Physics 16, 063012 (2014).
  * Young  _et al._ (2018) J. T. Young, T. Boulier, E. Magnan, E. A. Goldschmidt, R. M. Wilson, S. L. Rolston, J. V. Porto, and A. V. Gorshkov, Phys. Rev. A 97, 023424 (2018).
  * Madjarov  _et al._ (2020) I. S. Madjarov, J. P. Covey, A. L. Shaw, J. Choi, A. Kale, A. Cooper, H. Pichler, V. Schkolnik, J. R. Williams, and M. Endres, Nature Physics 16, 857 (2020).
  * Henkel  _et al._ (2010) N. Henkel, R. Nath, and T. Pohl, Phys. Rev. Lett. 104, 195302 (2010).
  * Keating  _et al._ (2013) T. Keating, K. Goyal, Y.-Y. Jau, G. W. Biedermann, A. J. Landahl, and I. H. Deutsch, Phys. Rev. A 87, 052314 (2013).
  * Petrosyan and Mølmer (2014) D. Petrosyan and K. Mølmer, Physical Review Letters 113, 123003 (2014).
  * Glaetzle  _et al._ (2015) A. W. Glaetzle, M. Dalmonte, R. Nath, C. Gross, I. Bloch, and P. Zoller, Physical Review Letters 114, 173002 (2015).
  * Li  _et al._ (2018) Y. Li, A. Geißler, W. Hofstetter, and W. Li, Phys. Rev. A 97, 023619 (2018).
  * Wu  _et al._ (2021) X. Wu, X. Liang, Y. Tian, F. Yang, C. Chen, Y.-C. Liu, M. K. Tey, and L. You, Chinese Physics B 30, 020305 (2021).
  * Gil  _et al._ (2014) L. I. R. Gil, R. Mukherjee, E. M. Bridge, M. P. A. Jones, and T. Pohl, Phys. Rev. Lett. 112, 103601 (2014).
  * van Bijnen and Pohl (2015) R. M. W. van Bijnen and T. Pohl, Phys. Rev. Lett. 114, 243002 (2015).
  * Mukherjee (2019) R. Mukherjee, Phys. Rev. A 100, 013403 (2019).
  * Keleş  _et al._ (2020) A. Keleş, E. Zhao, and X. Li, Phys. Rev. A 101, 023624 (2020).
  * Walker and Saffman (2008) T. G. Walker and M. Saffman, Phys. Rev. A 77, 032723 (2008).
  * Johnson and Rolston (2010) J. E. Johnson and S. L. Rolston, Physical Review A - Atomic, Molecular, and Optical Physics 82, 033412 (2010).
  * Honer  _et al._ (2010) J. Honer, H. Weimer, T. Pfau, and H. P. Büchler, Phys. Rev. Lett. 105, 160404 (2010).
  * Mukherjee  _et al._ (2016) R. Mukherjee, T. C. Killian, and K. R. A. Hazzard, Phys. Rev. A 94, 053422 (2016).
  * Reinhard  _et al._ (2007) A. Reinhard, T. C. Liebisch, B. Knuffman, and G. Raithel, Phys. Rev. A 75, 032712 (2007).
  * Rutkevich (2008) S. B. Rutkevich, Journal of Statistical Physics 131, 917 (2008).
  * Letscher and Petrosyan (2018) F. Letscher and D. Petrosyan, Phys. Rev. A 97, 043415 (2018).
  * Robertson  _et al._ (2021) E. Robertson, N. Sibalic, R. Potvliege, and M. Jones, Computer Physics Communications 261, 107814 (2021).
  * Steinert  _et al._ (2023) L.-M. Steinert, P. Osterholz, R. Eberhard, L. Festa, N. Lorenz, Z. Chen, A. Trautmann, and C. Gross, Phys. Rev. Lett. 130, 243001 (2023).
  * Hollerith  _et al._ (2022) S. Hollerith, K. Srakaew, D. Wei, A. Rubio-Abadal, D. Adler, P. Weckesser, A. Kruckenhauser, V. Walther, R. van Bijnen, J. Rui, C. Gross, I. Bloch, and J. Zeiher, Phys. Rev. Lett. 128, 113602 (2022).
  * Schauß  _et al._ (2012) P. Schauß, M. Cheneau, M. Endres, T. Fukuhara, S. Hild, A. Omran, T. Pohl, C. Gross, S. Kuhr, and I. Bloch, Nature 491, 87 (2012).
  * Schauß (2018) P. Schauß, Quantum Science and Technology 3, 023001 (2018).
  * Browaeys and Lahaye (2020) A. Browaeys and T. Lahaye, Nature Physics 16, 132 (2020).
  * Labuhn  _et al._ (2014) H. Labuhn, S. Ravets, D. Barredo, L. Béguin, F. Nogrette, T. Lahaye, and A. Browaeys, Phys. Rev. A 90, 023415 (2014).
  * Bernien  _et al._ (2017) H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Nature 551, 579 (2017).
  * Li  _et al._ (2022) X. X. Li, J. B. You, X. Q. Shao, and W. Li, Phys. Rev. A 105, 032417 (2022).
  * Schempp  _et al._ (2015) H. Schempp, G. Günter, S. Wüster, M. Weidemüller, and S. Whitlock, Phys. Rev. Lett. 115, 093002 (2015).
  * Yang  _et al._ (2019) F. Yang, S. Yang, and L. You, Phys. Rev. Lett. 123, 063001 (2019).
  * Liu  _et al._ (2020) F. Liu, S. Whitsitt, P. Bienias, R. Lundgren, and A. V. Gorshkov, arXiv:2007.07258 (2020).
  * Omran  _et al._ (2019) A. Omran, H. Levine, A. Keesling, G. Semeghini, T. T. Wang, S. Ebadi, H. Bernien, A. S. Zibrov, H. Pichler, S. Choi, J. Cui, M. Rossignolo, P. Rembold, S. Montangero, T. Calarco, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Science 365, 570 (2019).


# 3D Orientation Estimation with Multiple 5G mmWave Base Stations

Mohammad A. Nazari1  Gonzalo Seco-Granados2  Pontus Johannisson3  Henk
Wymeersch1
1Department of Electrical Engineering  Chalmers University of Technology
Sweden
2Department of Telecommunications and Systems Engineering  Universitat
Autonoma de Barcelona  Spain
3RISE Research Institutes of Sweden  Sweden
email<EMAIL_ADDRESS>
###### Abstract

We consider the problem of estimating the 3D orientation of a user, using the
downlink mmWave signals received from multiple base stations. We show that the
received signals from several base stations, having known positions, can be
used to estimate the unknown orientation of the user. We formulate the
estimation problem as a maximum likelihood estimation problem in the the
manifold of rotation matrices. In order to provide an initial estimate to
solve the non-linear non-convex optimization problem, we resort to a least
squares estimation problem that exploits the underlying geometry. Our
numerical results show that the problem of orientation estimation can be
solved when the signals from at least two base stations are received. We also
provide the orientation lower error bound, showing a narrow gap between the
performance of the proposed estimators and the bound.

##  I Introduction

5G mmWave signals can provide accurate location information by virtue of their
large bandwidth and large arrays at transmitter and receiver [1, 2]. In
particular, in contrast to 4G localization, which requires several
synchronized base stations (BSs), in 5G mmWave, a single BS can be sufficient
to obtain a location fix in 3D. In order to localize using a single BS,
related work on 5G mmWave localization assumed perfect synchronization between
the user equipment (UE) and the BS [2, 3, 4, 5, 6]. In practice, this
assumption is not valid, prompting work on joint localization and
synchronization by exploiting multipath information [7, 8, 9]. In these works,
the environment is mapped while at the same time the UE is localized and
synchronized to the BS.

In addition to localization and synchronization, determining the orientation
of the user is important for purposes of beamforming and beam tracking [10,
11]. Orientation information is obtained from angle measurements at the user
side (i.e, angle of arrival (AoA) in downlink and angle of departure (AoD) in
uplink). As was demonstrated in [2], with a single BS, the UE orientation in
2D can be determined when the UE is equipped with an array. In contrast, in 3D
a single BS does not suffice to determine the UE orientation, since only
angles in azimuth and elevation can be measured, unless additional signal
sources are available, e.g., reflectors or scatterers [3], reconfigurable
intelligent surfaces [12], or additional BSs. Joint localization and
orientation estimation was also considered in [13], for anchor-free swarm
navigation system, in [14] for visible light positioning, in [15] for
realistic channel realizations with hybrid array architectures, and in [16]
for relative localization of vehicles. Orientation estimation from range
measurements is also possible, as described in [17]. However, in the above
contributions, the reliance on range measurements requires tight
synchronization. To avoid this, pure angle-based localization or orientation
estimation methods were pursued in [18].

In this paper, we consider a mmWave MIMO scenario where the downlink signals
are used to estimate the orientation of a UE with known position. While
rotation estimation and tracking is a problem that has seen extensive
treatment in the robotics literature [19], generally relying on an inertial
measurement unit (accelerometers and rate gyros), our formulation is unique as
it provides _absolute_ 3D orientation information. The major contributions of
this paper are as follows:

  * •

We derive the orientation error bound for the estimation of the rotation
matrix of a UE, using a constrained Fisher information analysis, which gives
the lower bound for the performance of any unbiased estimator.

  * •

We pose a least squares (LS) optimization problem on the manifold of rotation
matrices, giving a solution when the downlink signals from at least two BSs
yield AoA measurements at the UE.

  * •

The obtained LS estimate can then be refined using a maximum likelihood
optimization on the manifold of rotation matrices, leading to a solution with
lower root mean squared error, and approaching the orientation error bound.

Our results show that the proposed estimation algorithms are efficient,
approaching the orientation error bound, with low complexity.

#### Notations

We denote vectors and matrices with bold lowercase and uppercase letters
($\mathbf{x}$ and $\mathbf{X}$), respectively. The matrix $\mathbf{X}$
transpose is represented as $\mathbf{X}^{\top}$. We write the Kronecker
product as $\otimes$, and the cardinality of a set $\mathcal{M}$ as
$|\mathcal{M}|$.

##  II Problem Statement

###  II-A System Model

We consider a scenario where there are $M$ base stations (BSs) with known
positions and known orientations, in a global coordinate system, to be used
for the estimation of the orientation a user equipment (UE), in the
environment. We denote by $\mathbf{p}_{m}=[p_{x,m},p_{y,m},p_{z,m}]^{\top}$,
$m=1,\cdots,M$ the positions of the base stations. We further assume that the
position of the UE is known and equal to
$\mathbf{p}=[p_{x},p_{y},p_{z}]^{\top}$, while the orientation
$\mathbf{o}=[\alpha,\beta,\gamma]^{\top}$ is unknown and to be estimated111The
angels $\alpha$, $\beta$, and $\gamma$ are called Euler, or Tait-Bryan
angles.. The UE orientation determines a local frame of reference,
conveniently described by a $3\times 3$ rotation matrix in the special
orthogonal group $\mathrm{SO}(3)$: $\mathbf{R}\in\mathrm{SO}(3)$ (i.e., an
orthogonal matrix, satisfying $\mathbf{R}^{\top}\mathbf{R}=\mathbf{I}_{3}$ and
$\det(\mathbf{R})=+1$). While angles $\alpha$, $\beta$, and $\gamma$ have
physical meaning, they can be related to the rotation matrix $\mathbf{R}$.
Accordingly, we mean _estimating the rotation matrix $\mathbf{R}$_ when
referring to _3D orientation estimation_.

The rotation order is important when mapping between the Euler angles and the
rotation matrix. We consider the sequence of rotations around $z$,
$y^{\prime}$, and $x^{\prime\prime}$, as the following:

$\displaystyle\mathbf{R}=\mathbf{R}_{z}(\alpha)\mathbf{R}_{y}(\beta)\mathbf{R}_{x}(\gamma),$
(1)

where $\mathbf{R}_{z}(\alpha)$ denotes a rotation of $\alpha$ radians around
the Z-axis

$\displaystyle\mathbf{R}_{z}(\alpha)=\left[\begin{array}[]{ccc}\cos(\alpha)&-\sin(\alpha)&0\\\
\sin(\alpha)&\cos(\alpha)&0\\\ 0&0&1\end{array}\right],$ (5)

$\mathbf{R}_{y}(\beta)$ shows a rotation of $\beta$ radians around the Y-axis

$\displaystyle\mathbf{R}_{y}(\beta)=\left[\begin{array}[]{ccc}\cos(\beta)&0&\sin(\beta)\\\
0&1&0\\\ -\sin(\beta)&0&\cos(\beta)\end{array}\right],$ (9)

and $\mathbf{R}_{x}(\gamma)$ indicates a rotation of $\gamma$ radians around
the X-axis

$\displaystyle\mathbf{R}_{x}(\gamma)=\left[\begin{array}[]{ccc}1&0&0\\\
0&\cos(\gamma)&-\sin(\gamma)\\\
0&\sin(\gamma)&\cos(\gamma)\end{array}\right].$ (13)

The system model is visualized in Fig. 1.

###  II-B Signal Model

The BSs and the UE are equipped with arrays, which are capable of measuring
the angles of departure (AoD) in the BS side, and angles of arrival (AoA) in
the UE side. As a reference, we consider $\mathbf{R}=\mathbf{I}_{3}$ to
correspond to the UE being parallel with the XY plane, with axis aligned with
X and Y axis respectively. The signal observed by the UE is of the form

$\displaystyle\mathbf{y}_{t}=\sum_{m=1}^{M}\alpha_{m}\mathbf{a}(\bm{\theta}_{m})\mathbf{a}^{\top}(\bm{\psi}_{m})\mathbf{s}_{m,t}+\mathbf{n}_{t},\,t=1,\ldots,T,$
(14)

where $\alpha_{m}$ is the complex channel gain from BS $m$ to the UE,
$\mathbf{a}(\bm{\theta}_{m})$ is the UE response vector corresponding to AoA
$\bm{\theta}_{m}=[\theta^{(\mathrm{el})}_{m},\theta^{(\mathrm{az})}_{m}]^{\top}$
for elevation angle $\theta^{(\mathrm{el})}_{m}$ and azimuth angle
$\theta^{(\mathrm{az})}_{m}$. Similarly, $\bm{\psi}_{m}$ denotes the AoD in
elevation and azimuth from BS $m$. The transmitted signal by BS $m$ is
$\mathbf{s}_{m,t}$ (with power
$P_{m}=\mathbb{E}\\{\|\mathbf{s}_{m,t}\|^{2}\\}$) and $\mathbf{n}_{t}$ is
spatially and temporally white complex Gaussian noise with variance $N_{0}/2$
per real dimension. Since the UE location is known, the AoDs are known as
well. Under the considered model, it follows immediately that

$\displaystyle\theta^{(\mathrm{el})}_{m}$
$\displaystyle=\arccos\left(q_{z,m}/\|\mathbf{q}_{m}\|\right),$ (15a)
$\displaystyle\theta^{(\mathrm{az})}_{m}$ $\displaystyle=\arctan
2(q_{y,m},q_{x,m}),$ (15b)

with

$\displaystyle\mathbf{q}_{m}=\mathbf{R}^{\top}(\mathbf{p}_{m}-\mathbf{p}).$
(16)

The AoAs in the local coordinate system of the UE are also shown in Fig. 1.

###  II-C Measurement Model

We assume an estimator exists that determines estimates of $\bm{\theta}_{m}$
from the observation $\mathbf{y}_{t}$, $t=1,\ldots,T$ (e.g., see [2]). For
simplicity of the exposition, we will model these AoA estimates as mutually
independent with von Mises distributions, i.e.,

$\displaystyle p(\hat{\bm{\theta}}|\bm{\theta})$
$\displaystyle=\prod_{m=1}^{M}\frac{1}{2\pi
I_{0}(\kappa^{(\mathrm{el})}_{m}))}\exp(\kappa^{(\mathrm{el})}_{m}\cos(\hat{{\theta}}^{(\mathrm{el})}_{m}-{{\theta}}^{(\mathrm{el})}_{m})),$
(17) $\displaystyle\times\frac{1}{2\pi
I_{0}(\kappa^{(\mathrm{az})}_{m}))}\exp(\kappa^{(\mathrm{az})}_{m}\cos(\hat{{\theta}}^{(\mathrm{az})}_{m}-{{\theta}}^{(\mathrm{az})}_{m}))$
(18)
$\displaystyle\propto\exp(\bm{\kappa}^{\top}\cos(\hat{\bm{\theta}}-\bm{\theta}))$
(19)

where $I_{0}(\cdot)$ is the modified Bessel function of order $0$,
$\kappa^{(\mathrm{el})}_{m}$ is the concentration parameter of the $m$-th AoA
in elevation, $\kappa^{(\mathrm{az})}_{m}$ is the concentration parameter of
the $m$-th AoA in azimuth. In (19), we have overloaded the notation for
cosines and aggregated the AoAs, their estimates, and the corresponding
concentrations in the $2M\times 1$ vectors $\bm{\theta}$, $\hat{\bm{\theta}}$
and $\bm{\kappa}$. The concentration parameters $\bm{\kappa}$ depend on the
quality of the estimator. It is important to note that the AoAs are obtained
in the local frame of reference of the UE, which depends on the UE orientation
$\mathbf{R}$.

Figure 1: Schematic of 3D orientation estimation using downlink mmWave MIMO
signals from 2 BSs.

##  III Fisher Information Analysis

###  III-A Background on (Constrained) Fisher Information

When estimating an unknown vector $\bm{\eta}\in\mathbb{R}^{N}$, constrained to
lie on a manifold $\mathbf{h}(\bm{\eta})=\mathbf{0}$ defined by $K\geq 0$ non-
redundant constraints, from an observation $\mathbf{y}$, the error covariance
(under certain technical conditions) is lower bounded as [20]

$\displaystyle\mathbb{E}\left\\{(\bm{\eta}-\hat{\bm{\eta}})(\bm{\eta}-\hat{\bm{\eta}})^{\top}\right\\}$
$\displaystyle\succeq\mathcal{I}_{\mathrm{const}}^{-1}(\bm{\eta}),$ (20)

where

$\displaystyle\mathcal{I}_{\mathrm{const}}^{-1}(\bm{\eta})$
$\displaystyle=\mathbf{M}(\mathbf{M}^{\top}\mathcal{I}(\bm{\eta})\mathbf{M})^{-1}\mathbf{M}^{\top},$
(21)

in which $\mathcal{I}(\bm{\eta})\in\mathbb{R}^{N\times N}$ is the
unconstrained Fisher information matrix

$\displaystyle[\mathcal{I}(\bm{\eta})]_{n,n^{\prime}}=-\mathbb{E}\left\\{\frac{\partial^{2}}{\partial\eta_{n}\partial\eta_{n^{\prime}}}\log
p(\mathbf{y}|\bm{\eta})\right\\},$ (22)

and $\mathbf{M}\in\mathbb{R}^{N\times(N-K)}$ with
$\mathbf{M}^{\top}\mathbf{M}=\mathbf{I}_{N-K}$, satisfying

$\displaystyle\frac{\partial\mathbf{h}(\bm{\eta})}{\partial\bm{\eta}^{\top}}\mathbf{M}=\mathbf{0}_{K\times(N-K)},$
(23)

is obtained by collecting the orthonormal basis vectors of null-space of the
gradient matrix
${\partial\mathbf{h}(\bm{\eta})}/{\partial\bm{\eta}^{\top}}\in\mathbb{R}^{K\times
N}$. Note that as a special case without constraints ($K=0$), we obtain the
standard Fisher information matrix inequality
$\mathbb{E}\left\\{(\bm{\eta}-\hat{\bm{\eta}})(\bm{\eta}-\hat{\bm{\eta}})^{\top}\right\\}\succeq\mathcal{I}^{-1}(\bm{\eta})$.

###  III-B Measurement FIM

Given a likelihood of the form (19), the Fisher information of $\bm{\theta}$
is

$\displaystyle\mathcal{I}(\bm{\theta})=\mathrm{diag}(\bm{\kappa}\odot
I_{1}(\bm{\kappa})\oslash I_{0}(\bm{\kappa})),$ (24)

where $\odot$ and $\oslash$ denote pointwise product and division,
respectively, and $I_{1}(\cdot)$ is the modified Bessel function of order $1$.
The proof is provided in Appendix A.

###  III-C Orientation FIM

To obtain the Fisher information of the rotation matrix $\mathbf{R}$, we
vectorize $\mathbf{R}=[\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{3}]$ as

$\displaystyle\mathbf{r}=\mathrm{vec}(\mathbf{R})=[\mathbf{r}_{1}^{\top},\mathbf{r}_{2}^{\top},\mathbf{r}_{3}^{\top}]^{\top}.$
(25)

The unconstrained Fisher information matrix
$\mathcal{I}(\mathbf{r})\in\mathbb{R}^{9\times 9}$ of $\mathbf{r}$ is then
obtained using the transformation matrix relating the measurements to the
elements of the unknown rotation matrix $\mathbf{R}$ as in the following:

$\displaystyle\mathcal{I}(\mathbf{r})=\bm{\Upsilon}\mathcal{I}(\bm{\theta})\bm{\Upsilon}^{\top},$
(26)

where

$\displaystyle\bm{\Upsilon}=\frac{\partial\bm{\theta}}{\partial\mathbf{r}},$
(27)

with the obtained elements as in Appendix B. However, to account for the
orthogonality constraint of the rotation matrix, i.e.,
$\mathbf{R}^{\top}\mathbf{R}=\mathbf{I}_{3}$, we have $K=6$ constraints

$\displaystyle\mathbf{h}(\mathbf{r})=[$
$\displaystyle\|\mathbf{r}_{1}\|^{2}-1,\mathbf{r}^{\top}_{2}\mathbf{r}_{1},\mathbf{r}^{\top}_{3}\mathbf{r}_{1},$
$\displaystyle\|\mathbf{r}_{2}\|^{2}-1,\mathbf{r}^{\top}_{2}\mathbf{r}_{3},\|\mathbf{r}_{3}\|^{2}-1]^{\top}=\mathbf{0}_{6\times
1}.$ (28)

The following matrix $\mathbf{M}$ is an orthonormal basis for the null-space
of the gradient matrix
${\partial\mathbf{h}(\bm{\eta})}/{\partial\bm{\eta}^{\top}}$ [17]:

$\displaystyle\mathbf{M}=\begin{bmatrix}-\mathbf{r}_{3}&\mathbf{0}_{3\times
1}&\mathbf{r}_{2}\\\ \mathbf{0}_{3\times 1}&-\mathbf{r}_{3}&-\mathbf{r}_{1}\\\
\mathbf{r}_{1}&\mathbf{r}_{2}&\mathbf{0}_{3\times 1}\end{bmatrix}.$ (29)

6

Finally, we can define the orientation error bound (OEB) as

$\displaystyle\mathrm{OEB}$
$\displaystyle=\sqrt{\mathrm{trace}(\mathcal{I}_{\mathrm{const}}^{-1}(\mathbf{r}))}$
(30)
$\displaystyle\leq\sqrt{\mathbb{E}\left\\{\|\mathbf{r}-\hat{\mathbf{r}}\|\right\\}}=\sqrt{\mathbb{E}\left\\{\|\mathbf{R}-\hat{\mathbf{R}}\|^{2}_{F}\right\\}},$
(31)

where $\|.\|_{F}$ is the Frobenius norm.

##  IV Methodology

In this section, we first describe the general principle of optimization over
the $\mathrm{SO}(3)$ manifold. Then, we describe a method to obtain an initial
estimate of the UE rotation, based on a least squares criterion. This estimate
is then refined through the maximum likelihood criterion.

###  IV-A Optimization on the $\mathrm{SO}(3)$ Manifold

To solve problems of the form

$\displaystyle\hat{\mathbf{R}}=\arg\min_{\mathbf{R}\in\mathrm{SO}(3)}f(\mathbf{R}),$
(32)

where $f:\mathrm{SO}(3)\to\mathbb{R}$ is a smooth function, we rely on the
method provided in [21, Chapter 4]. Starting from an initial estimate
$\hat{\mathbf{R}}^{(0)}$, we compute

$\displaystyle\hat{\mathbf{R}}^{(k+1)}=\mathrm{Ret}_{\hat{\mathbf{R}}^{(k)}}\left(-\varepsilon_{k}\mathrm{Proj}_{\hat{\mathbf{R}}^{(k)}}\left.\frac{\partial
f(\mathbf{R})}{\partial\mathbf{R}}\right|_{\mathbf{R}=\hat{\mathbf{R}}^{(k)}}\right),$
(33)

where $\mathrm{Proj}_{\mathbf{X}}(\cdot)$ is a projection onto the tangent
space (the set of real, skew-symmetric $3\times 3$ matrices) at $\mathbf{X}$,
$\mathrm{Ret}(\cdot)$ is a retraction from the tangent space onto
$\mathrm{SO}(3)$, and $\varepsilon_{k}>0$ is a suitable step size.
Intuitively, the gradient is calculated, projected to the tangent space (to
follow the space of $\mathrm{SO}(3)$ as closely as possible), the initial
matrix is updated, and then the updated matrix is normalized back into the
$\mathrm{SO}(3)$ space. The projection and retractions operations are given by
[21, eqs. (7.36) and (7.22)]

$\displaystyle\mathrm{Proj}_{\mathbf{X}}(\mathbf{U})$
$\displaystyle=\mathbf{X}\mathrm{skew}(\mathbf{X}^{\top}\mathbf{U}),$ (34)
$\displaystyle\mathrm{Ret}_{\mathbf{X}}(\mathbf{U})$
$\displaystyle=(\mathbf{X}+\mathbf{U})(\mathbf{I}_{3}+\mathbf{U}^{\top}\mathbf{U})^{-1/2},$
(35)

where $\mathrm{skew}(\mathbf{Z})=(\mathbf{Z}-\mathbf{Z}^{\top})/2$. It can be
verified that $\mathrm{Ret}_{\mathbf{X}}(\mathbf{U})\in\mathrm{SO}(3)$ when
$\mathbf{U}$ belongs to the tangent space at $\mathbf{X}\in\mathrm{SO}(3)$.
Hence, optimization on the manifold requires definition of an initial
estimate, the cost function $f(\mathbf{R})$ and its unconstrained gradient
${\partial f(\mathbf{R})}/{\partial\mathbf{R}}$.

###  IV-B Least Squares (LS) Estimation

According to (15), we can obtain an estimate
$\hat{\mathbf{q}}_{m}(\hat{\bm{\theta}}_{m})$ of ${\mathbf{q}}_{m}$ from the
estimated AoA $\hat{\bm{\theta}}_{m}$ from BS $m$ in azimuth and elevation,
since all the distances to the different BSs are assumed known. Hence, when
considering several such estimates, we can solve for $\mathbf{R}$ using the
relation (16). Hence, we use the following procedure. First, we select a
subset $\mathcal{M}\subseteq\\{1,\ldots,M\\}$ of BSs. From
$\hat{\bm{\theta}}_{m}$, $m\in\mathcal{M}$, we compute
$\hat{\mathbf{q}}_{m}(\hat{\bm{\theta}}_{m})$ from the relations

$\displaystyle\hat{q}_{z,m}=\|\mathbf{p}_{m}-\mathbf{p}\|\cos(\hat{\theta}^{(\mathrm{el})}_{m}),$
(36)
$\displaystyle\hat{q}_{y,m}=\|\mathbf{p}_{m}-\mathbf{p}\|\sin(\hat{\theta}^{(\mathrm{el})}_{m})\sin(\hat{\theta}^{(\mathrm{az})}_{m}),$
(37)
$\displaystyle\hat{q}_{x,m}=\|\mathbf{p}_{m}-\mathbf{p}\|\sin(\hat{\theta}^{(\mathrm{el})}_{m})\cos(\hat{\theta}^{(\mathrm{az})}_{m}).$
(38)

We then create a matrix $\mathbf{Q}\in\mathbb{R}^{3\times|\mathcal{M}|}$ that
contains the estimates $\hat{\mathbf{q}}_{m}(\hat{\bm{\theta}}_{m})$ as
columns, so that $\mathbf{Q}\in\mathbb{R}^{3\times|\mathcal{M}|}$. Finally, we
solve the LS problem (initialized with the identity matrix)

$\displaystyle\hat{\mathbf{R}}_{\mathrm{LS}}=$
$\displaystyle\arg\min_{\mathbf{R}\in\mathrm{SO}(3)}\|\mathbf{U}-\mathbf{R}\mathbf{Q}\|^{2}_{F},$
(39)

where
$\mathbf{U}=\mathbf{P}-\mathbf{p}\otimes\mathbf{1}^{\top}_{|\mathcal{M}|}$, in
which $\mathbf{P}=[\mathbf{p}_{m\in\mathcal{M}}]$. The gradient of the cost
function is

$\displaystyle\frac{\partial\|\mathbf{U}-\mathbf{R}\mathbf{Q}\|^{2}_{F}}{\partial\mathbf{R}}=-2\|\mathbf{U}-\mathbf{R}\mathbf{Q}\|_{F}{\mathbf{Q}}^{\top}.$
(40)

######  Remark 1.

The complexity of the method above grows with the number of used BS
$\mathcal{M}$. To reduce the complexity, since $\mathbf{R}$ has only three
degrees of freedom, at least 2 BSs should be used, so $|\mathcal{M}|=2$
suffices. These 2 BSs can be chosen based on their relative geometry (e.g.,
not colinear with the UE) and the concentration values in the von Mises
likelihood.

###  IV-C Maximum Likelihood (ML) Estimation

The ML estimate of $\mathbf{R}$ is obtained by maximizing the log-likelihood
function.

$\displaystyle\hat{\mathbf{R}}_{\mathrm{ML}}=\arg\min_{\mathbf{R}\in\mathrm{SO}(3)}\underbrace{-\bm{\kappa}^{\top}\cos(\hat{\bm{\theta}}-\bm{\theta}(\mathbf{R}))}_{=f(\mathbf{R})},$
(41)

initialized with the LS estimate (39). The gradient of the log-likelihood
function is given by

$\displaystyle\frac{\partial f(\mathbf{R})}{\partial\mathbf{R}}=$
$\displaystyle-\sum_{m=1}^{M}{\kappa}^{(\mathrm{el})}_{m}\sin(\hat{{\theta}}^{(\mathrm{el})}_{m}-{\theta}^{(\mathrm{el})}_{m}(\mathbf{R}))\frac{\partial{\theta}^{(\mathrm{el})}_{m}(\mathbf{R})}{\partial\mathbf{R}}$
$\displaystyle-\sum_{m=1}^{M}{\kappa}^{(\mathrm{az})}_{m}\sin(\hat{{\theta}}^{(\mathrm{az})}_{m}-{\theta}^{(\mathrm{az})}_{m}(\mathbf{R}))\frac{\partial{\theta}^{(\mathrm{az})}_{m}(\mathbf{R})}{\partial\mathbf{R}}.$
(42)

The gradients
$\frac{\partial}{\partial\mathbf{R}}{\theta}^{(\mathrm{el})}_{i}(\mathbf{R})$
and
$\frac{\partial}{\partial\mathbf{R}}{\theta}^{(\mathrm{az})}_{i}(\mathbf{R})$
are provided in Appendix B.

##  V Numerical Results

In this section, we analyse performance of the proposed estimators, and
compare it with the OEB from (30). We also show the orientation error bound
for a range of UE orientations.

Figure 2: OEB in dB ($10\log_{10}(\mathrm{OEB})$) for $\beta=-\pi/4$ vs.
$\alpha,\gamma\in[0,\pi]$ with $M=2$ BSs. Figure 3: OEB in dB
($10\log_{10}(\mathrm{OEB})$) for $\beta=-\pi/4$ vs. $\alpha,\gamma\in[0,\pi]$
with $M=3$ BSs.

###  V-A Simulation Scenario

Unless otherwise noted, we consider a scenario with $M=2$ BSs located at
$\mathbf{p}_{1}=[0,0,0]^{\top}$ and $\mathbf{p}_{2}=[0,50,0]^{\top}$ sending
downlink signals to the UE located at $\mathbf{p}=[50,0,-5]^{\top}$. The UE is
equipped with a uniform planar array with $16\times 16$ antenna elements with
half-wavelength spacing. The carrier frequency is 28 GHz and the transmitted
signals are set as
$\mathbf{s}_{m,t}=\sqrt{P_{m}}\mathbf{a}^{*}(\bm{\psi}_{m})/\|\mathbf{a}(\bm{\psi}_{m})\|$,
so that we can define the SNR as

$\displaystyle\mathrm{SNR}_{m}=\frac{|\alpha_{m}|^{2}TN_{m,\mathrm{tx}}P_{m}}{N_{0}},$
(43)

where $N_{m,\mathrm{tx}}$ denotes the number of transmit antennas at BS $m$,
assuming coherent combining across transmissions. Without loss of generality
we set $\mathrm{SNR}_{1}=\mathrm{SNR}_{2}=\mathrm{SNR}$. The user orientation
is set via the angles $\mathbf{o}=[\alpha,\beta,\gamma]^{\top}$. To obtain the
measurements and their likelihoods, we proceed as follows:

  * •

We derive the FIM of $\mathcal{I}_{\bm{Y}}(\bm{\theta})$ from
$\mathbf{y}_{t}$, $t=1,\ldots,T$ under the considered SNR. We use the
subscript $\bm{Y}$ to express that the FIM is based on the received waveforms
$\mathbf{y}_{t}$.

  * •

We equate $\mathcal{I}^{-1}(\bm{\theta})$ to
$\mathrm{diag}(\mathcal{I}^{-1}_{\bm{Y}}(\bm{\theta}))$.

  * •

We equate the appropriate diagonal elements in $\mathcal{I}(\bm{\theta})$ to
${\kappa}^{(\mathrm{el})}_{m}I_{1}({\kappa}^{(\mathrm{el})}_{m})/I_{0}({\kappa}^{(\mathrm{el})}_{m})$
and
${\kappa}^{(\mathrm{az})}_{m}I_{1}({\kappa}^{(\mathrm{az})}_{m})/I_{0}({\kappa}^{(\mathrm{az})}_{m})$,
for $m=1,\ldots,M$, as derived in (24). We then solve for
${\kappa}^{(\mathrm{el})}_{m}$ and ${\kappa}^{(\mathrm{az})}_{m}$.

All manifold optimization problems were solved with the Manopt toolbox [22].

###  V-B Results and Discussion

We first evaluate the impact of the orientation and the number of BSs on the
OEB. To visualize this, we fix $\beta$ to $-\pi/4$, and sweep $\alpha$ and
$\gamma$ in the range $[0,\pi]$. The SNR is set to $-10~{}\mathrm{dB}$. Fig. 2
shows the corresponding result, with OEB values larger than 1 truncated. We
observe low OEB for most UE orientations, but there are several peaks, where
the OEB tends to infinity. Specifically, when $\alpha\approx\pi/2$ and
$\gamma\approx\pi/4$, the received ray from one of the BS hits the UE antenna
array on the broadside, and does not provide a high quality orientation
estimation. However, once this specific orientation changes, the downlink
signal arrives at the UE array in a more suitable direction, facilitating a
more satisfactory estimation, and accordingly lower OEB. Fig. 3 depicts the
OEB, when a third BS at the position $\mathbf{p}_{3}=[50,50,0]^{\top}$ is
added. As observed, the OEB peaks are eliminated, and the error bound is
greatly reduced for all considered orientations. The third BS can help to
ensure that sufficient rays reach the UE with appropriate angles, not leading
to unidentifiable orientation estimation.

Figure 4: Performance of Estimators and Comparison with the OEB vs. SNR.

We now fix the orientation to $[\alpha,\beta,\gamma]=[0.6\pi,0,-0.8\pi]$ and
evaluate the performance of the proposed estimators as a function of the SNR.
From 200 Monte Carlo simulations, we obtain an estimate of
$({\mathbb{E}\\{\|\mathbf{R}-\hat{\mathbf{R}}\|^{2}_{F}\\}})^{1/2}$, for both
LS and ML estimators. We also plot the orientation error bound (21). As
observed, the performance of both estimators improves as SNR increases, which
is expected. The gap between the performance of ML and LS is due to the fact
that in LS, the estimation neglects the distribution of measurements (i.e.,
the concentration of the von Mises distribution). The gap between the ML and
LS estimator depends on the values of $\kappa$ for each AOA. When all $\kappa$
values are identical, the gap between ML and LS becomes smaller. When it comes
to comparison of the OEB with our proposed estimators, it is observed that the
performance of ML estimator achieves the OEB for
$\mathrm{SNR}>-20~{}\mathrm{dB}$. This shows the efficiency of the proposed
algorithms.

##  VI Conclusions

In this paper, we studied the problem of estimation of the orientation of a UE
in 3D using downlink mmWave MIMO signals from multiple base stations. The UE
has an array of antennas which allows measuring the AoAs from the BSs, and
solving for the unknown 3D orientation. Rotations are represented with
rotation matrices in $\mathrm{SO}(3)$, which allows us to formulate and solve
LS and ML estimation problems on the manifold of 3D rotations. The solution of
the LS problem was used as an initialization to the ML estimation problem,
which is a non-convex optimization on the manifold $\mathrm{SO}(3)$. The
performance of the resulting ML estimate coincides with the OEB, which is
obtained by using the FIM of the rotation matrix subject to its orthogonality
constraint. Future research would include the case where the position of the
UE is unknown, and the study of dynamic situations with time-varying UE
position and orientation.

## Acknowledgment

This work was supported under the Wallenberg AI, Autonomous Systems and
Software Program (WASP), the Swedish Research Council under grant 2018-03701,
the Spanish Ministry of Science, Innovation and Universities under Projects
TEC2017-89925-R and PRX18/00638 and by the ICREA Academia Programme.

##  Appendix A FIM of Von Mises Distribution

We focus on one AoA, allowing us to remove all indices. The log-likelihood
function is then given by

$\displaystyle\log p(\hat{\theta}|\theta)$
$\displaystyle=\kappa\cos(\hat{\theta}-\theta)-\ln 2\pi I_{0}(\kappa).$ (44)

Then we have

$\displaystyle\mathcal{I}(\theta)$
$\displaystyle=\mathbb{E}\left(-\frac{\partial^{2}}{\partial\theta^{2}}\log
p(\hat{\theta}|\theta)\right)$ (45)
$\displaystyle=\mathbb{E}\left(\kappa\cos(\hat{\theta}-\theta)\right)$ (46)
$\displaystyle=\int_{\theta-\pi}^{\theta+\pi}\kappa\cos(\hat{\theta}-\theta)\frac{\mathrm{e}^{\kappa\cos(\hat{\theta}-\theta)}}{2\pi
I_{0}(\kappa)}\mathrm{d}\hat{\theta}.$ (47)

By change of the variable $\hat{\theta}$, the equation (47) is simplified as

$\displaystyle\mathcal{I}(\theta)$
$\displaystyle=\int_{-\pi}^{\pi}\kappa\cos(\hat{\theta})\frac{\mathrm{e}^{\kappa\cos(\hat{\theta})}}{2\pi
I_{0}(\kappa)}\mathrm{d}\hat{\theta},$ (48)

which in turn is analogous to the first-order Bessel function

$\displaystyle
I_{1}(\kappa)=\int_{0}^{\pi}\frac{1}{\pi}\cos(\hat{\theta})\mathrm{e}^{\kappa\cos(\hat{\theta})}\mathrm{d}\hat{\theta}.$
(49)

Since $\cos(\hat{\theta})=\cos(-\hat{\theta})$, one can write

$\displaystyle\mathcal{I}(\theta)$
$\displaystyle=2\int_{0}^{\pi}\kappa\cos(\hat{\theta})\frac{\mathrm{e}^{\kappa\cos(\hat{\theta})}}{2\pi
I_{0}(\kappa)}\mathrm{d}\hat{\theta}$ (50)
$\displaystyle=\frac{\kappa}{I_{0}(\kappa)}\int_{0}^{\pi}\frac{1}{\pi}\cos(\hat{\theta})\mathrm{e}^{\kappa\cos(\hat{\theta})}\mathrm{d}\hat{\theta}$
(51) $\displaystyle=\kappa\frac{I_{1}(\kappa)}{I_{0}(\kappa)}.$ (52)

##  Appendix B Gradient of the AoA with respect to the rotation matrix

We introduce

$\displaystyle\mathbf{u}^{(m)}$
$\displaystyle=\frac{\mathbf{p}-\mathbf{p}_{m}}{\|\mathbf{p}-\mathbf{p}_{m}\|},$
(53)

and $\mathbf{u}_{1}=[1,0,0]^{\top}$, $\mathbf{u}_{2}=[0,1,0]^{\top}$,
$\mathbf{u}_{3}=[0,0,1]^{\top}$. This allows us to express (15a)–(15b) as

$\displaystyle\theta^{(\mathrm{el})}_{m}$
$\displaystyle=\arccos(-\mathbf{u}_{3}^{\top}\mathbf{R}^{\top}\mathbf{u}^{(m)}),$
(54) $\displaystyle\theta^{(\mathrm{az})}_{m}$ $\displaystyle=\arctan
2(-\mathbf{u}_{2}^{\top}\mathbf{R}^{\top}\mathbf{u}^{(m)},-\mathbf{u}_{1}^{\top}\mathbf{R}^{\top}\mathbf{u}^{(m)}),$
(55)

We now make use of the following identities

$\displaystyle\frac{\partial\mathbf{a}^{\top}\mathbf{X}^{\top}\mathbf{b}}{\partial\mathbf{X}}$
$\displaystyle=\mathbf{b}\mathbf{a}^{\top},$ (56)
$\displaystyle\frac{\partial}{\partial\mathbf{x}}\arccos({u}(\mathbf{x}))$
$\displaystyle=-\frac{1}{\sqrt{1-{u}^{2}(\mathbf{x})}}\frac{\partial}{\partial\mathbf{x}}{u}(\mathbf{x}),$
(57) $\displaystyle\frac{\partial}{\partial\mathbf{x}}\arctan
2({u}(\mathbf{x}),{v}(\mathbf{x}))$
$\displaystyle=\frac{{v}(\mathbf{x})\frac{\partial}{\partial\mathbf{x}}{u}(\mathbf{x})-{u}(\mathbf{x})\frac{\partial}{\partial\mathbf{x}}{v}(\mathbf{x})}{{u}^{2}(\mathbf{x})+{v}^{2}(\mathbf{x})},$
(58)

in order to write

$\displaystyle\frac{\partial\theta^{(\mathrm{el})}_{m}}{\partial\mathbf{R}}$
$\displaystyle=\frac{\mathbf{u}^{(m)}\mathbf{u}_{3}^{\top}}{\sqrt{1-(\mathbf{u}_{3}^{\top}\mathbf{R}^{\top}\mathbf{u}^{(m)}))^{2}}},$
(59)
$\displaystyle\frac{\partial\theta^{(\mathrm{az})}_{m}}{\partial\mathbf{R}}$
$\displaystyle=\frac{(\mathbf{u}_{1}^{\top}\mathbf{R}^{\top}\mathbf{u}^{(m)})\mathbf{u}^{(m)}\mathbf{u}_{2}^{\top}-(\mathbf{u}_{2}^{\top}\mathbf{R}^{\top}\mathbf{u}^{(m)})\mathbf{u}^{(m)}\mathbf{u}_{1}^{\top}}{(\mathbf{u}_{1}^{\top}\mathbf{R}^{\top}\mathbf{u}^{(m)})^{2}+(\mathbf{u}_{2}^{\top}\mathbf{R}^{\top}\mathbf{u}^{(m)})^{2}}.$
(60)

From this, we immediately obtain
${\partial\theta^{(\mathrm{el})}_{m}}/{\partial\mathbf{r}}=\mathrm{vec}({\partial\theta^{(\mathrm{el})}_{m}}/{\partial\mathbf{R}})$
and
${\partial\theta^{(\mathrm{az})}_{m}}/{\partial\mathbf{r}}=\mathrm{vec}({\partial\theta^{(\mathrm{az})}_{m}}/{\partial\mathbf{R}})$.

## References

  * [1] J. A. del Peral-Rosado, R. Raulefs, J. A. López-Salcedo, and G. Seco-Granados, “Survey of cellular mobile radio localization methods: From 1G to 5G,” vol. 20, no. 2, pp. 1124–1148, 2018.
  * [2] A. Shahmansoori, G. E. Garcia, G. Destino, G. Seco-Granados, and H. Wymeersch, “Position and orientation estimation through millimeter-wave MIMO in 5G systems,” _IEEE Transactions on Wireless Communications_ , vol. 17, no. 3, pp. 1822–1835, 2018.
  * [3] Z. Abu-Shaban, X. Zhou, T. Abhayapala, G. Seco-Granados, and H. Wymeersch, “Error bounds for uplink and downlink 3D localization in 5G millimeter wave systems,” _IEEE Transactions on Wireless Communications_ , vol. 17, no. 8, pp. 4939–4954, 2018.
  * [4] A. Kakkavas, M. H. Castañeda García, R. A. Stirling-Gallacher, and J. A. Nossek, “Performance limits of single-anchor millimeter-wave positioning,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 11, pp. 5196–5210, 2019.
  * [5] J. Talvitie, M. Valkama, G. Destino, and H. Wymeersch, “Novel algorithms for high-accuracy joint position and orientation estimation in 5G mmWave systems,” in _IEEE Global Communications Conference (GLOBECOM) Workshops_ , 2017.
  * [6] R. Mendrzik, H. Wymeersch, G. Bauch, and Z. Abu-Shaban, “Harnessing NLOS components for position and orientation estimation in 5G millimeter-wave MIMO,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 1, pp. 93–107, 2019.
  * [7] H. Wymeersch, N. Garcia, H. Kim, G. Seco-Granados, S. Kim, F. Wen, and M. Fröhle, “5G mmWave downlink vehicular positioning,” in _IEEE Global Communications Conference (GLOBECOM)_ , 2018.
  * [8] R. Mendrzik, F. Meyer, G. Bauch, and M. Z. Win, “Enabling situational awareness in millimeter-wave massive MIMO systems,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 13, no. 5, pp. 1196–1211, 2019.
  * [9] Y. Ge, F. Wen, H. Kim, M. Zhu, F. Jiang, S. Kim, L. Svensson, and H. Wymeersch, “5G SLAM using the clustering and assignment approach with diffuse multipath,” _Sensors_ , vol. 20, no. 16, p. 4656, 2020.
  * [10] J. Zhao, F. Gao, W. Jia, S. Zhang, S. Jin, and H. Lin, “Angle domain hybrid precoding and channel tracking for millimeter wave massive MIMO systems,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 10, pp. 6868–6880, 2017.
  * [11] J. Bao and H. Li, “Motion aware beam tracking in mobile millimeter wave communications: A data-driven approach,” in _IEEE International Conference on Communications (ICC)_ , 2019.
  * [12] A. Elzanaty, A. Guerra, F. Guidi, and M.-S. Alouini, “Reconfigurable intelligent surfaces for localization: Position and orientation error bounds,” _arXiv preprint arXiv:2009.02818_ , 2020.
  * [13] R. Poehlmann, S. Zhang, A. Dammann, and P. A. Hoeher, “Fundamental limits for joint relative position and orientation estimation,” in _IEEE International Conference on Communications (ICC) Workshops_ , 2018.
  * [14] S. Shen, S. Li, and H. Steendam, “Simultaneous position and orientation estimation for visible light systems with multiple leds and multiple pds,” _IEEE Journal on Selected Areas in Communications_ , vol. 38, no. 8, pp. 1866–1879, 2020.
  * [15] W. Zheng and N. González-Prelcic, “Joint position, orientation and channel estimation in hybrid mmWave MIMO systems,” in _2019 53rd Asilomar Conference on Signals, Systems, and Computers_ , 2019.
  * [16] A. Kakkavas, M. H. C. Garcia, R. A. Stirling-Gallacher, and J. A. Nossek, “Multi-array 5G V2V relative positioning: Performance bounds,” in _IEEE Global Communications Conference (GLOBECOM)_ , 2018.
  * [17] S. P. Chepuri, G. Leus, and A. van der Veen, “Rigid body localization using sensor networks,” _IEEE Transactions on Signal Processing_ , vol. 62, no. 18, pp. 4911–4924, 2014.
  * [18] Y. Wu, B. Peng, H. Wymeersch, G. Seco-Granados, A. Kakkavas, M. H. C. Garcia, and R. A. Stirling-Gallacher, “Cooperative localization with angular measurements and posterior linearization,” in _IEEE International Conference on Communications (ICC) Workshops_ , 2020.
  * [19] T. D. Barfoot, “State estimation for robotics,” 2019.
  * [20] P. Stoica and Boon Chong Ng, “On the Cramér-Rao bound under parametric constraints,” _IEEE Signal Processing Letters_ , vol. 5, no. 7, pp. 177–179, 1998.
  * [21] N. Boumal, “An introduction to optimization on smooth manifolds,” Available online, Aug 2020. [Online]. Available: http://www.nicolasboumal.net/book
  * [22] N. Boumal, B. Mishra, P.-A. Absil, and R. Sepulchre, “Manopt, a Matlab toolbox for optimization on manifolds,” _Journal of Machine Learning Research_ , vol. 15, no. 42, pp. 1455–1459, 2014. [Online]. Available: https://www.manopt.org