\end{align}

The temporal Fourier modes of the LDHO kernel are obtained from (<ref>) leading to (<ref>).
In light of (<ref>) and taking account  the dispersion relations, the fast and slow decay times transform as follows

\begin{align}
\tau_{s}(k) = & \frac{2\tilde{\tau}_{c}}{B(k) \left(1- 2\tilde{\tau}_{c} \tilde{\omega}_{d} \right)}, \;
\\
\tau_{f}(k) = & \frac{2\tilde{\tau}_{c}}{B(k) \left(1 + 2\tilde{\tau}_{c} \tilde{\omega}_{d} \right)}\,.
\end{align}

Based on (<ref>),  the functions $B(k)$ in (<ref>) cancel out, and  the temporal Fourier modes are given by
_-ω(k,τ) = c_0  ^-ϵk/4ω̃_d τ̃_c  [ β_f  ^ -β_s |τ|B(k)  /2τ̃_c  - β_s   ^ -β_f |τ|B(k) /2τ̃_c ],
where $\beta_{s}=1- 2\tilde{\tau}_{c} \tilde{\omega}_{d}$ and $\beta_{f}=1 + 2\tilde{\tau}_{c} \tilde{\omega}_{d}$.
Recalling (<ref>) for $B(k)$, the IFT expression (<ref>), and the linearity of the IFT, it follows that

\begin{align}
\label{eq:C-over-lin}
C(r,\tau)= & \frac{c_{0}}{4\tilde{\omega}_{d} \tilde{\tau}_{c} } \left[ \beta_{f}  C_{s}(r,\tau)
-  \beta_{s}  C_{f}(r,\tau)\,  \right],
\end{align}

C_j(r,τ) = ^ -β_j |τ|/2τ̃_c   [ ^ -ξ β_j |τ| k /2τ̃_c  -ϵk  ] =
^ -β_j |τ|/2τ̃_c/(2π)^d/2  r^ν ∫_0^∞k k^d/2 J_ν(kr) ^-a_jk ,   j=s, f ,
where $a_{j}=\frac{\beta_{j} \xi \lvert \tau \rvert}{2\tilde{\tau}_{c}} + \epsilon$.
The integral above can be calculated using <cit.>:

∫_0^∞k   k^d/2 J_ν(kr) ^-a_jk =2a_j(2r)^ν Γ(ν+3/2)/√(π)( a^2_j+ r^2)^(d+1)/2 ,   where   ν=d/2 -1 .

In view of the above, the functions $C_{j}(r,\tau)$ in (<ref>) are given by

C_j(r,τ) =  ^ -β_j |τ|/2τ̃_c
1/( a^2_j + r^2)^(d+1)/2

Finally, the overdamped LDHO covariance kernel (<ref>) is obtained from (<ref>) and (<ref>).

An example of the overdamped kernel $C(r,\tau)$  is shown in Fig. <ref>.

As in the underdamped case,  the spectral decay hyperparameter $\epsilon$ stabilizes the variance (i.e., the behavior at $\tau=0$),  and $\xi$ adjusts the space-time interaction since for $\xi=0$ the space and time dependence in (<ref>) decouple.

Zero-lag marginal covariances
The spatial and temporal marginal kernels  are  obtained from (<ref>) by setting $\tau=0$ and $r=0$ respectively, following simple algebraic calculations. Thus we obtain

\begin{align}
\label{eq:ldho-over-linear-marginal-r}
\Cms(r) =& \frac{c_{0}^\ast \,\left( \beta_{f} -\beta_{s} \right)\, }{ \left(  \epsilon^2 + r^2 \right)^{(d + 1)/2}} ,
\end{align}

\begin{align}
\label{eq:ldho-over-linear-marginal-t}
\Cmt(\tau) = &  \frac{c_{0}^\ast \,\beta_{f} \, \E^{-\frac{\beta_{s} \lvert \tau \rvert}{2\tilde{\tau}_{c}}}}{ \left(\frac{\xi \beta_{s} \lvert \tau \rvert}{2\tilde{\tau}_{c}} + \epsilon\right)^{d + 1}}   -   \frac{c_{0}^\ast \, \beta_{s} \, \E^{-\frac{\beta_{f} \lvert \tau \rvert}{2\tilde{\tau}_{c}}}}{ \left(\frac{\xi \beta_{f} \lvert \tau \rvert}{2\tilde{\tau}_{c}} + \epsilon\right)^{d + 1}}  \,.
\end{align}

The spatial marginal kernel (<ref>) is proportional to $(\epsilon^2 + r^2)^{-(d+1)/2}$ as in the underdamped case (<ref>).  The temporal marginal kernel (<ref>) comprises a combination of slow and fast  exponential kernels. This is analogous to the purely temporal case (<ref>), albeit  the coefficients of the exponentials are renormalized and include temporal dependence.

Normalized $C(r,\tau)$ and isolevel contour lines (red online) in the overdamped regime, obtained from (<ref>) using
$\tilde{\omega}_{d}=\pi/10$, $\tilde{\tau}_{c}=0.8$, $\xi=0.4$,  $\epsilon=8$, and $d=2$.

§.§ Critical damping

In this regime it holds that $\tilde{\omega}_{0} \tilde{\tau}_{c}=1/2$.  The temporal Fourier modes, based on (<ref>) and the dispersion relations, are given by

_-ω(k,τ) = c_0 ^ -|τ|B(k)/2τ̃_c -ϵk [ 1 + |τ| B(k)/2τ̃_c ].

In this case, the LDHO spatiotemporal kernel generated by the IFT of the temporal Fourier modes (<ref>).
The spatiotemporal LDHO kernel in the critical damping regime is given by (<ref>).

The LDHO temporal Fourier modes  are obtained  from (<ref>) by inserting the dispersion functions (<ref>) and (<ref>).
Recalling (<ref>) for $B(k)$, it follows that
_-ω(k,τ) =c_0 ^ -|τ|(1+ξk)/2τ̃_c -ξk  ( 1+ |τ|/2τ_c + ξ k|τ|/2τ_c).

Based on the linearity of the IFT we obtain

\begin{align}
\label{eq:C-crit-lin}
C(r,\tau)= & c_{0}\E^{ -\frac{\lvert \tau \rvert}{2\tilde{\tau}_{c}}}\, \left( 1+ \frac{\lvert \tau \rvert}{2\tau_c}\right) C_{1}(r,\tau)
+ \frac{c_{0} \xi  \,\lvert \tau \rvert}{2\tilde{\tau}_{c}} \, \E^{ -\frac{\lvert \tau \rvert}{2\tilde{\tau}_{c}}}\,  C_{2}(r,\tau),\end{align}

where $C_{1}(r,\tau)$ and $C_{2}(r,\tau)$ are given by

\begin{align}
\label{eq:C1-crit-ift-lin}
C_{1}(r,\tau)= & \ift_{\bfk}\left[ \, \E^{-k\left(\epsilon + \frac{\xi\lvert \tau \rvert}{2\tilde{\tau}_{c}} \right)} \right],
\\
\label{eq:C2-crit-ift-lin}
C_{2}(r,\tau)= & \ift_{\bfk}\left[ \, k \,\E^{-k\left(\epsilon + \frac{\xi\lvert \tau \rvert}{2\tilde{\tau}_{c}} \right)} \right].
\end{align}

Based on (<ref>), the $\ift_{\bfk}$ (<ref>) can be expressed as

\[
C_{1}(r,\tau)=\frac{1}{(2\pi)^{d/2}\, r^\nu} \int_{0}^{\infty}\D k k^{d/2}\,J_{\nu}(kr)\,\E^{-a_{\Re}k}, \; \textrm{where} \; a_{\Re}=\epsilon + \frac{\xi\lvert \tau \rvert}{2\tilde{\tau}_{c}} \,.
\]

The integral in the above equation has been evaluated in (<ref>), based on which we obtain

C_1(r,τ)= Γ(d+1/2)/π^(d+1)/2 ( r^2 + a_^2)^(d+1)/2 .

Similarly, the IFT in (<ref>) can be expressed as

\[
C_{2}(r,\tau)=\frac{1}{(2\pi)^{d/2}\, r^\nu} \int_{0}^{\infty}\D k k^{d/2 +1}\,J_{\nu}(kr)\,\E^{-a_{\Re}k}, \; \textrm{where} \; a_{\Re}=\epsilon + \frac{\xi\lvert \tau \rvert}{2\tilde{\tau}_{c}} \,.
\]

The latter can be expressed in terms of $C_{1}(r,\tau)$ by noticing that the only difference between the two is an extra factor of $k$ in the integral for $C_{2}(r,\tau)$. Therefore, by taking advantage that $C_{2}(r,\tau)$ depends on $a_{\Re}$ only via $\E^{-a_{\Re}k}$, it follows that

\[
C_{2}(r,\tau)= - \frac{\partial C_{1}(r,\tau)}{\partial a_{\Re}}\,.
\]

Thus, we obtain by differentiating $C_{1}(r,\tau)$ in (<ref>)

C_2(r,τ)= Γ(d+1/2)  (d+1)  a_/π^(d+1)/2 ( r^2 + a_^2)^(d+3)/2 .
This concludes the proof.

The critically damped LDHO kernel (<ref>)  includes four independent  hyperparameters: $c_{0}, \tilde{\tau}_{c}, \epsilon, \xi$;  $\tilde{\omega}_{d}=0$ at critical damping.  The critical-damping kernel (<ref>) can be viewed as the limit of the overdamped kernel (<ref>) for $\tilde{\omega}_{d} \to 0$, which implies $\beta_{s} \to 1$,  $\beta_{f} \to 1$.  The comments in Remark <ref> regarding the role of $\xi$ and $\epsilon$ also hold for the critically damped case.

Zero-lag marginal covariances
The spatial and temporal marginal kernels  are  obtained from (<ref>) by setting $\tau=0$ and $r=0$ respectively. We thus obtain

\begin{align}
\label{eq:ldho-crit-marginal-r-linear}
\Cms(r)  & =  c_{0}\, C_{1}(r, 0)=
\frac{c_{0} \Gamma(\frac{d+1}{2})}{\pi^{(d+1)/2} \left( r^{2} + \epsilon^{2}\right)^{(d+1)/2}}\,,
\end{align}

\begin{align}
\label{eq:ldho-crit-marginal-t-linear}
\Cmt(\tau) & = c_{0}\, \frac{\,\Gamma(\frac{d+1}{2})}{ \pi^{(d+1)/2}\, a_{\Re}^{d+1}}\, \E^{-\frac{\lvert \tau \rvert}{2\tilde{\tau}_{c}}}\, \left[ \left( 1 + \frac{\lvert \tau \rvert}{2\tilde{\tau}_{c}}\right)   + \frac{\xi\lvert \tau \rvert}{2\tilde{\tau}_{c}} \, \frac{(d+1) }{a_{\Re}}\right]\,.
\end{align}

§ S4. KERNELS BASED ON THE ORNSTEIN-UHLENBECK ODE

The covariance of the O-U process is given by
$C(\tau)=\sigma^2\, \exp(-\lvert \tau \rvert/\tau_c)$ where $\sigma^2 =\sigma_{\eta}^{2}\tau_{c}/{2}$ <cit.>. The radial dispersion relations are given by
$\sigma^{2} \to \sigma^{2}(k) = \sigma_{0}^{2} A(k), \, \tau_{c} \to \tau_{c}(k)= \tilde{\tau}_{c}/B(k)$.  The O-U temporal Fourier modes for radial dispersion functions are thus given by

\begin{align}
\label{eq:t-ft-modes-O-U}
\tC_{-\omega}(k,\tau) &  = \sigma^{2}_{0}  \, A(k) \, \exp\left[ -\lvert \tau\rvert\,  B(k)/\tilde{\tau}_{c}\right].
\end{align}

Since $A(k)$ is dimensionless, $[\sigma_{0}^{2}]=[\mathrm{X}]^{2}[\mathrm{L}]^d$, where $\mathrm{L}$ represents length, so that the FT (<ref>) be dimensionally correct.
Based on the IFT (<ref>), the O-U covariance kernel is given by the integral $(\nu=d/2-1)$:

\begin{equation}
\label{eq:invcovft-O-U}
C({r},\tau)= \frac{\sigma^{2}_{0}}{(2\pi)^{d/2} r^{\nu}} \int_{0}^{\infty}  k^{d/2}\,
{J_{\nu}(k  r)}  A(k) \, \E^{ -\lvert \tau\rvert\,  B(k)/\tilde{\tau}_{c}} \, \D k\,.
\end{equation}

In the following, we derive spatiotemporal  kernel expressions for two different choices of dispersion functions.

* $A(k)=\E^{-\beta k^2}, \; B(k)=a+bk^2$ where $a,b, \beta>0$ are hyperparameters with units $[b]=[\beta]=[L]^2$, $[a]=[L]^{0}$:

The spectral integral (<ref>) becomes

\[
C({r},\tau)= \frac{\sigma^{2}_{0}\,
\E^{- a\,\lvert \tau\rvert/\tilde{\tau}_{c}}}{(2\pi)^{d/2}\, r^{\nu}}  \int_{0}^{\infty}  k^{d/2}    {J_{\nu}(k  r)} \, \E^{ -b k^{2} \lvert \tau\rvert\,  /\tilde{\tau}_{c} - \beta k^2} \, \D k\,.
\]

Using the table of integrals <cit.> it follows that

C(r,τ)= σ^2_0  ^- a |τ|/τ̃_c/(2π)^d/2 ^- r^2 / 4( β+ b |τ|/τ̃_c ) /( 2b|τ|/τ̃_c+2β)^d/2  .
The kernel (<ref>) involves four free hyperparameters: $\sigma_{0}$, $\tilde{\tau}_{c}/a$, $\tilde{\tau}_{c}/b$ and $\beta$. The spatial and temporal marginal kernels are given respectively by

\begin{align}
\Cms(r) = & \frac{\sigma^{2}_{0}\, }{(4\pi\, \beta)^{d/2}}\,\E^{-r^{2}/4\beta}\,,
\\
\Cmt(\tau) = & \frac{\sigma^{2}_{0}\, \left(  1 + b\,\lvert \tau \rvert/\beta \tilde{\tau}_c\right)^{-d/2}}{(4\pi\, \beta)^{d/2}}\,\E^{-a\,\lvert \tau\rvert/\tilde{\tau}_{c}}\,  \,.
\end{align}

Hence, the spatial marginal covariance is given by the square exponential kernel while the temporal marginal is a modified exponential kernel.

* $A(k)=\E^{-\beta k}$, $B(k)=a + \xi\,k$, where $a,\xi,\beta>0$  are hyperparameters with units $[b]=[\beta]=[L]$, $[a]=[L]^{0}$:

The spectral  integral (<ref>) becomes ($\nu=d/2-1$):

\[
C({r},\tau)= \frac{\sigma^{2}_{0}\,
\E^{- a\,\lvert \tau\rvert/\tilde{\tau}_{c}}}{(2\pi)^{d/2}\, r^{\nu}}  \int_{0}^{\infty}  k^{\nu+1}    {J_{\nu}(k r)} \, \E^{ - k \left(\beta +\xi \,  \lvert \tau\rvert\,  /\tilde{\tau}_{c} \right)} \, \D k\,.
\]

Using the table of integrals <cit.> we obtain

C(r,τ)= σ^2_0  Γ(d+1/2) /π^(d+1)/2( βτ̃_c +ξ |τ|  )  ^- a |τ|/τ̃_c/τ̃_c  [ r^2 + ( β+ ξ|τ|/τ̃_c)^2 ]^(d+1)/2  .
The kernel (<ref>) involves four free hyperparameters: $\sigma_{0}$, $\tilde{\tau}_{c}/a$, $\tilde{\tau}_{c}/\xi$ and $\beta$. The spatial and temporal marginal kernels are given respectively by

\begin{align}
\Cms(r) = & \frac{\sigma^{2}_{0}\, \Gamma(\frac{d+1}{2})\,}{\pi^{(d+1)/2}}\frac{ \beta }{\left( r^{2} +  \beta^2 \right)^{(d+1)/2}} \,\,,
\\
\Cmt(\tau) = & \frac{\sigma^{2}_{0}\, \Gamma(\frac{d+1}{2})\,}{\pi^{(d+1)/2}}\frac{\,\E^{- a\,\lvert \tau\rvert/\tilde{\tau}_{c}}}{\left( \beta + \frac{\xi\lvert \tau\rvert}{\tilde{\tau}_{c}}\right)^{d}} \,.
\end{align}

The Ornstein-Uhlenbeck covariance kernels with square and linear $k$ dependence of the dispersion functions are illustrated in Fig. <ref>.

Normalized $C(r,\tau)$ and isolevel contour lines (red online) of the Ornstein-Uhlenbeck kernels in $d=2$ spatial dimensions. Left: obtained from the dispersion function (<ref>) with $k^2$ dependence using the  hyperparameters
$\tilde{\tau}_{c}=0.8$, $b=0.4$, $a=0.5$,  $\beta=8$.  Right: obtained from the dispersion functions (<ref>) with $k$ dependence using the  hyperparameters
$\tilde{\tau}_{c}=0.8$, $\xi=0.4$, $a=0.5$,  $\beta=8$.


subfig

# Fully-probabilistic Terrain Modelling with Stochastic Variational Gaussian
Process Maps

Ignacio Torroba1 Christopher Iliffe Sprague1 and John Folkesson1 1The authors
are with the Swedish Maritime Robotics Centre (SMaRC) and the Division of
Robotics, Perception and Learning at KTH Royal Institute of Technology, SE-100
44 Stockholm, Sweden {torroba, sprague<EMAIL_ADDRESS>
###### Abstract

Gaussian processes (GPs) are becoming a standard tool to build terrain
representations thanks to their capacity to model map uncertainty. This
effectively yields a reliability measure of the areas of the map, which can be
directly utilized by Bayes filtering algorithms in robot localization
problems. A key insight is that this uncertainty can incorporate the noise
intrinsic to the terrain surveying process through the GPs ability to train on
uncertain inputs (UIs). However, existing techniques to build GP maps with UIs
in a tractable manner are restricted in the form and degree of the input
noise. In this letter, we propose a flexible and efficient framework to build
large-scale GP maps with UIs based on Stochastic Variational GPs and Monte
Carlo sampling of the UIs distributions. We validate our mapping approach on a
large bathymetric survey collected with an AUV and analyze its performance
against the use of deterministic inputs (DI). Finally, we show how using UI
SVGP maps yields more accurate particle filter localization results than DI
SVGP on a real AUV mission over an entirely predicted area.

###### Index Terms:

Mapping, Marine Robotics, Localization, Gaussian Process

##  I Introduction

Seabed surveying with sonar-equipped, autonomous underwater vehicles (AUVs),
as opposed to surface vessels, provides higher-resolution bathymetric maps due
to a closer proximity of the AUV to the sea bottom. This higher resolution,
however, comes at the expense of a cumulative drift in the vehicle’s dead
reckoning (DR) estimate due to the lack of GPS signal underwater. Among
sonars, multibeam echosounders (MBES) have become standard in bathymetric
surveying thanks to the ease of interpreting the 3D data they provide, which
can be treated through existing point cloud manipulation techniques [1].
However, they have limited resolution and the beam-forming process is
inherently prone to noise. Due to this, bathymetric point cloud maps tend to
be memory-intensive, noisy and have gaps in the coverage.

Gaussian processes (GPs) have been demonstrated to be an invaluable technique
to tackle some of these problems when building maps, as they can regress the
original seabed structure from the noisy, sparse MBES beams collected during
the survey. The result is a lighter representation of the terrain when
compared to raw point clouds, which naturally smooth out MBES noise and can
predict unseen areas conditioned on the collected data. Furthermore, being a
Bayesian inference technique, GPs provide a measure of the variance of the
resulting map. This results in maps which provide a reliability measure of the
terrain modelled, from which localization algorithms based on Bayes filters
can benefit directly.

Another capability of GPs is that they can fuse the uncertainties emanating
from the vehicle’s DR and the sensor used in the data collection into their
learning process. These can be propagated into the sonar beams’ uncertainty to
create uncertain inputs (UIs), which the GP can learn. The work in [2] shows
how the use of UIs yields more robust predictions. The most common technique
currently used for training GP maps with UIs was introduced in [3] and has
since been applied in several works in the field of autonomous mapping and
informative path planning [4, 5]. This method is based on fusing the UIs
uncertainties directly into the GP kernel through a Gauss-Hermite quadrature
of the covariance function. However it is limited by the need to manually
adjust the number of sampling points $M$ for the quadrature to adequately
approximate the degree of uncertainty in the data. Moreover, this formulation
imposes the restriction on the UIs to be Gaussian.

In this work, we propose as an alternative the use of Monte Carlo (MC)
sampling of the UIs together with a Stochastic Variational Gaussian Process
(SVGP) [6] in order to propagate the input uncertainties into the map. We show
that this approach can handle large amounts of training data, and neither
imposes any constraint on the UIs distribution nor requires any prior
knowledge on the degree of uncertainty in the inputs. We test this method on a
very large bathymetric survey with an AUV and measure the advantages of the
use of UIs against deterministic inputs (DIs) for learning SVGP maps in the
tasks of bathymetry reconstruction and prediction. Finally, to show the
amenability of SVGPs with UIs to practical applications, we demonstrate AUV
localization on an entirely predicted area. We compare the results against
localization on DI SVGP maps and show how the UIs SVGP predictions result in a
more accurate AUV pose estimate. The framework developed for this work has
been released and can be found in this
repository111https://github.com/ignaciotb/UWExploration/.

##  II Related work

Gaussian processes are a non-parametric, Bayesian regression method capable of
modelling real-world problems in a fully probabilistic manner. Originally,
their main limitations was their cubic scalability with the number of training
inputs and their limited applicability only to Gaussian likelihoods [7].
Regarding scalability, the introduction of a sparse formulation based on a
low-rank approximation to the kernel [8] granted a reduction both in
computational complexity and memory requirements of the GP implementations
[9]. Regarding applicability, the adaption of variational inference techniques
to learning GP models in [10] allowed posing the training as an optimization
problem and approximating non-Gaussian GP likelihoods. Building upon this
approach, [6] presented a stochastic minibatch optimization framework for GP
training capable of handling datasets of hundreds of thousands of samples, the
Stochastic Variational GP.

Concerning the treatment of noisy training data, the GP formulation originally
assumed the input locations to be deterministic (DIs) in order to make the
training tractable. This, however, hindered its application to many real-world
problems that dealt with noise-corrupted inputs. In [2], Girard et al. derived
a closed-form formulation for GP kernels that accounted for Gaussian
uncertainty in the inputs. Although successful in fusing the inputs
uncertainty into the GP training, this approach was only applicable to a few
amenable Gaussian kernels. The later heteroscedastic approaches [11, 12] would
link the inputs’ variances to the outputs proportionally through the local,
approximated linear gradient of the GP posterior mean.

In the field of autonomous terrain modelling, GPs with UIs have been applied
to map regression from noisy sensory inputs in several contexts. Following the
conclusions from [2], in [3] O’Callaghan et al. employed a Gaussian-Hermite
quadrature approximation to the GP kernel to fuse the UIs’ uncertainty into
the GP training in a tractable manner. They showed the advantages of using UIs
in a grid mapping task with a mobile robot equipped with a laser scanner.
However, this approximation required the UIs to be Gaussian and manually
specifying the number of samples for the quadrature a priori, which depends on
the degree of uncertainty in the UIs. This would necessitate prior knowledge
of the training data and a similar degree of uncertainty across the whole
dataset for an optimal number of samples to be found. Besides, quadrature
approximations are known to scale poorly with the dimensions of the input
data. Despite these shortcomings, this technique is still being used in the
latest works in the field [5, 4].

SVGPs have been applied to map modelling both in the land [13] and the
underwater domains [14]. Their preferred use over other GP modalities was
justified by the need for a fast learning procedure capable of handling large
datasets in an online fashion. However, in both cases the inputs were assumed
to be noise-free and therefore DIs were used. To the best of our knowledge,
this work is the first to present results on SVGP mapping with UIs by
leveraging the minibatch optimization technique for SVGP training with the
unbiased MC approximation of the distribution of the uncertain inputs.

Finally, there are several instances of vehicle localization on GP maps based
on particle filters (PF) in the literature. Most notably, in [15] both the
prediction and measurements functions are modelled through GPs. In the
underwater domain, [16] presents an AUV PF localization framework in a GP map
and compares its performance in simulation against different bathymetry
reconstructions. Going beyond the aforementioned works, in this letter we
present the results of a real AUV PF localization test mission on an entirely
predicted area of seabed on the SVGP map.

##  III SVGP bathymetric maps with UIs

In this section, our approach to learn SVGP maps with UIs is presented.
Although we describe the method with the view on bathymetric mapping with AUVs
and sonar sensors, the framework is general and could be applied to any
equivalent instance of autonomous terrain construction under localization and
sensor uncertainty.

###  III-A Autonomous Bathymetric Surveying

Given an AUV equipped with a 3D mapping sensor, the aim is to autonomously
reconstruct an unknown area of seabed. The first-order Markov dynamics of the
6D AUV state at time $t$ is given by $r_{t}=g(r_{t-1},c_{t})$ where $c_{t}$ is
the control input to the motion model. The MBES measurement model
$p_{i}=h(r_{t},e_{i})$ provides a 3D sample of terrain $p_{i}=[x,y,z]$ from
the patch $e_{i}$ where the beam $i$ hits the seabed for a given vehicle pose
$r_{t}$. If we model the inherent noise in the real DR and measurement models
as white Gaussian with $\nu_{t}\sim\mathcal{N}(0,W)$ and
$\delta_{t}\sim\mathcal{N}(0,Q)$ respectively (with $W\in\mathbb{R}^{6\times
6}$ and $Q\in\mathbb{R}^{3\times 3}$) both $r_{t}$ and $p_{i}$ follow
probability distributions given by:

$\displaystyle r_{t}$ $\displaystyle\sim\mathcal{N}(g(r_{t-1},c_{t}),W)$ (1)
$\displaystyle P_{t}$ $\displaystyle\sim\mathcal{N}(h(r_{t},E_{t}),Q),$ (2)

where the measurement model has been generalized to gather the $n$ beams
collected in the same ping $t$, $P_{t}=\\{p_{i}\\}_{i=1}^{n}$ corresponding to
a set of patches of seabed $E_{t}=\\{e_{i}\\}_{i=1}^{n}$ in the map frame.

After a survey in which $N$ sensor observations or samples from the terrain
have been collected into a 3D point cloud, a dataset of the form
$D=\\{p^{xy}_{i},p^{z}_{i}\\}_{i=1}^{N}$ can be constructed. In the remainder
of the paper we will also use the notation $D=\\{X,Y\\}$ in place of the
previous definition to follow the standard notation in GP literature.

###  III-B Stochastic Variational Gaussian Processes Maps

The dataset $D$ can be used to regress a function of the form $f:X\to Y$ with
a SVGP, which will produce a continuous representation of the underlying
terrain from which the samples have been collected. Additionally, the
posterior of $f(X)$ on the observations will allow us to perform predictions
on unseen data conditioned on $D$. Applying Bayes theorem this posterior can
be computed as:

$p(f(X)|Y)=\frac{p(Y|f(X))p(f(X))}{p(Y)}$ (3)

We choose our prior to be a Gaussian process such that
$p(f(X))\sim\mathcal{N}(0,\mathcal{K}_{NN})$ modelled by its kernel
$\mathcal{K}_{NN}$, with a set of hyperparameters $\theta$. However, computing
that posterior is intractable due to the marginal likelihood $p(Y)$ and
therefore SVGP applies variational inference to approximate $p(f(X)|Y)\approx
q(f(X))$. Where $q(f(X))\sim\mathcal{N}(\tau,\beta)$ is a Gaussian variational
term.

Given a specific model, the optimal $\theta$ can be learned maximizing the
likelihood of observing our $N$ collected terrain samples for that given
model. In variational inference this is equivalent to maximizing the ELBO of
Eq. 3 with respect to the model parameters. The ELBO is derived minimizing the
KL divergence between our target posterior and the variational term:
$\text{argmin}_{\theta}KL[q(f(X)||p(f(X)|Y)]$. However, this optimization
involves the inversion of the very large kernel matrix $\mathcal{K}_{NN}$. In
order to overcome this, sparse GPs formulation makes use of a low-rank
approximation to the kernel
$\mathcal{K}_{NN}\approx\mathcal{K}_{NS}\mathcal{K}_{SS}^{-1}\mathcal{K}_{SN}$
using $S<<N$ inducing variables $u_{s}=f(Z_{s})$ in the model. $Z_{s}$ are the
locations of the inducing points, which belong in the same domain as $X$ in
this work. With the use of the inducing variables, the new target KL
divergence to optimize becomes $KL[q(f(Z)||p(f(Z)|u)]$, which can be
approximated using Jensen’s inequality:

$\mathcal{L}=\sum^{N}_{i}\mathbb{E}_{q(f(x_{i}))}\left[\log p(y_{i}\mid
f(x_{i}))\right]-\\\ \mathrm{KL}\left[q(u)\mid\mid p(u)\right]$ (4)

All the components in the ELBO above are Gaussian and can be derived.
Therefore, standard stochastic gradient descent (SGD) techniques can be
applied to its optimization. This, together with the use of the inducing
points formulation, allows us to handle very large training datasets through
the use of minibatch optimization. This is another key reason to advocate for
the use of SVGPs when constructing large bathymetric maps containing hundreds
of thousands of samples.

Finally, the ELBO in Eq. 4 can be derived with respect to the locations of the
inducing points $Z$ and so these can be added as hyperparameters within
$\theta$ and learned jointly. Despite adding computational workload to the
training, learning the optimal location of the inducing points will result in
a more accurate representation of the bathymetry [10]. Thus, it should always
be considered if the training time is not a constraint (such as in an offline
context).

###  III-C Vehicle Uncertainty Propagation to the SVGP Map

So far we have assumed that our training data is deterministic. But for a
fully probabilistic approach to autonomous mapping, all the system
uncertainties in the process need to be taken into account. There are two main
sources of noise in surveying: those emanating from the vehicle’s DR and those
related to physical sensors. Thus, we first introduce how the AUV pose
estimate and MBES uncertainties, modelled in Eq. 1, are propagated to the
beams and fused. After that, we show how standard SVGP can be readily used to
handle the resulting UIs through the use of MC integration.

If we now redefine the training inputs as samples from the terrain corrupted
with Gaussian noise, we obtain:

$p_{i}=\bar{p}_{i}+\epsilon_{pi}\;\;\text{with}\;\;\epsilon_{pi}\sim\mathcal{N}(0,\Sigma_{i})$
(5)

Through this $p_{i}\sim\mathcal{N}(\bar{p_{i}},\Sigma_{i})$ is a Gaussian
distribution and we can model through $\Sigma_{i}$ the uncertainty propagated
from the data collection process. We compute $\bar{p_{i}}$ and $\Sigma_{i}$
through the sigmapoint method originally presented in [17]. We choose to do so
following the results in [18], which argues that this method provides the most
accurate results for uncertainty propagation through a nonlinear measurement
model. We implement our uncertainty propagation framework with the tools to
manipulate poses and associated uncertainties presented in [18]. The reader is
referred to [18] for a full explanation of the method.

For every DR pose estimate in the AUV trajectory $\bar{r_{t}}$, and the $n$
seabed patches observed from that pose,
$e_{i,t}\sim\mathcal{N}(\bar{e}_{i,t},\Omega)$, a covariance matrix
$\Xi_{t}\in\mathbb{R}^{9\times 9}$ is compounded from $W_{t}$ (Eq. 1) and
$\Omega$. Then, the $l$ sigma point samples $\\{r_{l},e_{l}\\}$, which are
used to approximate the mean and variance of each measurement
$\\{\bar{p}_{i},\Sigma_{i}\\}$, are computed as follows:

$\displaystyle CC^{T}$ $\displaystyle=cholesky(\Xi_{t})$ (6)
$\displaystyle\lambda_{l}$ $\displaystyle=0\;\;l=0$ (7)
$\displaystyle\lambda_{l}$
$\displaystyle=\sqrt{L+\kappa}\>col_{l}C,\;\;l=1,...,L$ (8)
$\displaystyle\lambda_{l+L}$
$\displaystyle=-\sqrt{L+\kappa}\>col_{l}C,\;\;l=1,...,L$ (9)
$\displaystyle\begin{bmatrix}\xi_{l}\\\ \zeta_{l}\end{bmatrix}$
$\displaystyle=\lambda_{l}$ (10) $\displaystyle r_{l}$
$\displaystyle=exp(\xi^{\wedge}_{{}_{l}})\bar{r_{t}}$ (11) $\displaystyle
e_{l}$ $\displaystyle=\bar{e_{i}}+D\zeta_{l},$ (12)

where $col_{l}$ is the $l$ column of the matrix, $\xi_{l}\in\mathbb{R}^{6}$
and $\zeta_{l}\in\mathbb{R}^{3}$ are the perturbations applied to the AUV pose
and seabed patches respectively, $\wedge$ and $exp()$ are the operators to
turn them back and forth into $\mathbb{R}^{4\times 4}$ matrices belonging to
the $SE(3)$ group, $D$ is a dilation matrix and $L=9$ and $\kappa=0$. Each
sample is passed through the nonlinear MBES model in Eq 1,
$p_{l}=h(r_{l},e_{l})$ with $l=0,...,2L$ and the desired covariances for each
UI are computed as follows:

$\displaystyle\bar{p}_{i}$ $\displaystyle=\frac{1}{L+\kappa}\left(\kappa
p_{0}+\frac{1}{2}\sum_{l=1}^{2L}p_{l}\right)$ (13)
$\displaystyle\begin{split}\Sigma_{i}&=\frac{1}{L+\kappa}\left(\kappa(p_{0}-\bar{p}_{i})(p_{0}-\bar{p}_{i})^{T}\right.\\\
&+\left.\frac{1}{2}\sum_{l=1}^{2L}(p_{l}-\bar{p}_{i})(p_{l}-\bar{p}_{i})^{T}\right)\end{split}$
(14)

Considering now the training data as distributions instead of samples, our
dataset becomes $D=\\{\langle
p^{xy}_{i},\Sigma^{xy}_{i}\rangle,p^{z}_{i}\\}_{i=1}^{N}$, or equivalently
$D=\\{X,Y\\}$ with $X\sim\mathcal{N}(p^{xy}_{i},\Sigma^{xy}_{i})$. The way the
SVGP can handle these UIs is once again thanks to its capacity to operate on
large datasets, together with the approximation capabilities of Monte Carlo
sampling: Now for every learning step, a minibatch is constructed with $m<<N$
deterministic training samples drawn from the $m$ randomly selected UIs
distributions $X_{i}$ in a Monte Carlo fashion (one sample per UI, per
minibatch iteration). As the SVGP training progresses, the MC sampling
approaches the UIs original distributions while propagating their noise into
the learning process through stochastic sampling. Since the incoming data in
SVGP is ultimately of deterministic form, the SVGP does not need any
reformulation to handle UIs and therefore any standard implementation can be
used for this approach.

##  IV Particle Filter Localization on SVGP Maps

This section briefly presents the application of SVGP maps to the task of
vehicle localization with a particle filter. PF localization on GP maps
requires an accurate terrain reconstruction but also well-callibrated map
variances for a correct vehicle pose estimate, since the GP-PF utilizes both.
We argue that our method’s harnessing of the uncertainty in the map creation
process provides better estimates of both parameters, resulting in a more
accurate PF pose estimate than through the use of DI SVGP maps.

###  IV-A The SVGP-PF Algorithm

A particle filter approximates the posterior distribution of a vehicle pose
estimate through a set of $J$ finite sample states given by Eq. 15, following
Algorithm 1.

$R_{t}=\\{\langle r^{j}_{t},w^{j}_{t}\rangle|j=1,...,J\\}$ (15)

Algorithm 1 The SVGP Particle filter

1:$(R_{t},c_{t,t-1},P_{t},Q_{t}):$

2:$\hat{R_{t}}=R_{t}=\emptyset$

3:for $j=1$ to $J$ do

4: $r_{t}^{j}\sim\mathcal{N}(g(r_{t-1},c_{t}),W)$

5: $\hat{P}_{t}^{j}=meas(r_{t}^{j},P_{t})$

6:
$w_{t}^{j}=\mathcal{N}(P^{z}_{t};GP_{\mu}(\hat{P}_{t}^{j,xy}),GP_{\Sigma}(\hat{P}_{t}^{j,xy})+Q_{P^{z}_{t}})$

7: $add$ $\langle r^{j}_{t},w^{j}_{t}\rangle$ $to$ $\hat{R}_{t}$

8:end for

9:for $j=1$ to $J$ do

10: $draw$ $i$ $with$ $probability$ $\propto w_{k}^{i}$

11: $add$ $r_{t}^{i}$ $to$ $R_{t}$

12:end for

13:return $R_{t}$

In algorithm 1, the function $meas$ in line 4 is a transformation of the real
MBES ping at time $t$ to the pose of particle $j$ given by $r_{t}^{j}$. The
weight in line 5 models the probability of the difference between the depth
components in the ping $t$ and those predicted by the SVGP for particle $j$
being zero. These predicted measurements can be computed through the SVGP
predictive posterior from Eq. 3 on query points $X_{*}$, which has the form:

$p(f(X_{*})|Y)\sim\mathcal{N}(K_{*s}K^{-1}_{ss}\mu,\\\
(K_{*s}K^{-1}_{ss})\Sigma(K_{*s}K^{-1}_{ss})^{T}+K_{**}-K_{*s}K^{-1}_{ss}K^{T}_{*s}))$
(16)

Having $X_{*}=\hat{P}^{j,xy}$, line 5 can be formulated as
$w_{t}^{j}=p((P^{z}_{t}-GP(\hat{P}_{t}^{j,xy}))=0)$. This probability
distribution is normal and naturally utilizes both the $z$ component of the
noise model of the pings $Q_{P^{z}_{t}}$ and the uncertainty predicted by the
SVGP $GP_{\Sigma}(\hat{P}_{t}^{j,xy})$.

##  V Experiments

The experiments below compare the performance of SVGP maps trained with UIs
with the presented method against DIs in the tasks of map modelling and
prediction for different levels of noise in the system. Additionally, PF
localization experiments are carried out over predicted terrain to assess the
performance of UI SVGP against DI SVGP maps on a real application.

###  V-A The Bathymetric Dataset

For the experiments, a bathymetric survey was collected off the Swedish west
coast with a MBES Kongsberg 2040 mounted on an AUV Hugin 3000, following the
procedure explained in Section III. The raw bathymetric point cloud can be
seen in Fig. 2 a) and Fig. 1 shows the Hugin being deployed on the area. The
full mission lasted approximately 12.7 hours and covered a trajectory of 30.5
$km$ at an average speed of 1.5 $m/s$ . The resulting dataset $D$ contains
approximately 4.038.000 MBES beams over an area of approximately 1.4 $km^{2}$.

Figure 1: Hugin AUV being deployed on the survey area

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

Figure 2: Left to right: a) full GT bathymetry, b) corrupted DR point cloud,
c) UI SVGP posterior mean trained with b), and d) consistency error map
between b) and the posterior mean in c) sampled.

###  V-B The SVGP Training

The SVGP maps have been implemented within the Gpytorch framework [19]. A
Matern covariance function [20] with a smoothness parameter of $\tau=1/2$ was
selected as the SVGP kernel, as in [21]. The optimizer applied has been Adam
[22] with a learning rate of $1e^{-1}$ and a minibatch size of $m=4000$ points
per iteration. The same exponential moving average over the ELBO [23] has been
used as a terminating criterion for the optimization of SVGPs with DIs and
UIs. The number of inducing points has been manually set to $S=1000$ based on
prior knowledge of the dataset and the capabilities of the hardware. Existing
research on methods to include this number as a hyperparameter in the learning
process [24, 13] should be considered in an online setup. Finally, the SVGPs
training has been carried out in a laptop with an NVIDIA GeForce GTX 1060 GPU.

###  V-C SVGP Mapping Evaluation

In order to directly compare the SVGP maps trained with DIs and UIs against
the real bathymetry, we would need access to the ground truth map of the
seabed. Since this is generally not an option in underwater surveying, for our
experiments we have used the collected data as ground truth. Although not
ideal, this is a reasonable assumption based on the high quality of the DR
system on Hugin and the MBES model used. In order to generate lower quality
surveys, we have corrupted the vehicle DR with additive white noise and
propagated the vehicle and sensor uncertainties to the resulting bathymetric
point clouds for the GP training.

To assess the quality of the SVGPs trained under DR distortion when compared
to the GT bathymetry, we use the method presented in [25] and previously used
in [26]. This method computes the consistency error between the original
bathymetry and the point cloud generated from sampling the SVGP posterior on a
grid over the survey area with $100^{2}$ samples. We use the error between the
GT map and the raw point cloud generated from the disrupted survey as the
baseline error to compare the SVGP maps. This is because this error is
intrinsic to the data collection process and cannot be eliminated through
different map modelling techniques.

In order to evaluate both the reconstruction and prediction errors of the
SVGPs with DI and UIs, we have divided the survey into two sections. The
training dataset contains the outer loops of the survey, as seen in Fig. 2 b)
and consists of 2.650.000 MBES beams. The SVGPs have been trained on this data
and the reconstruction error has been measured directly comparing the SVGPs
posteriors to that section of the raw GT point cloud. The second survey
section contains the remaining MBES beams over the central area of the map.
This subset has not been used for training, only for validation of the SVGP
predictions over a large, unseen area. Thus, the prediction consistency error
has been computed by comparing the SVGP posterior point clouds to the full GT
bathymetry and subtracting the reconstruction error.

TABLE I: Average results comparing the performance of SVGP mapping with DIs and UIs. DR iid yaw noise (rad/s) | Reconstruction RMSE (m) | Prediction RMSE (m) | Learning steps | Tr($K_{ss}$)
---|---|---|---|---
DR PCL | SVGP Di | SVGP Ui | SVGP Di | SVGP Ui | SVGP Di | SVGP Ui | SVGP Di | SVGP Ui
0 | 0.750 | 1.226 | 1.179 | 1.660 | 1.678 | 8014.0 | 8015.0 | 58375.59 | 59230.035
$1e^{-3}$ | 2.354 | 2.461 | 2.437 | 0.766 | 0.566 | 8512.6 | 9045.8 | 60001.240 | 61666.679
$2e^{-3}$ | 3.366 | 3.309 | 3.308 | 0.697 | 0.554 | 8859.4 | 9189.2 | 66936.702 | 68230.707

###  V-D Experimental Setup

The framework used to replay, disrupt and visualize the AUV mission for the
uncertainty propagation and PF experiments has been implemented in ROS. Figure
Fig. 3 shows the visualization in RViz of the PF running on a SVGP map. The GT
Hugin pose is marked by the orange marker, with the particles represented in
grey.

Since we do not have access to the original Hugin localization uncertainty,
this has been approximated when generating the training datasets. The DR
uncertainty has been computed through an EKF with a 6 DOF point-mass motion
model, using as inputs the interpolated AUV velocities between time steps.
Regarding the PF experiments, the actual implementation of Eq. 16 of the SVGP
measurement model omits the computational loop in Algorithm 1 to improve the
computational complexity. Namely, instead of executing $J$ individual calls to
the GPU to compute the particles expected measurements, all the query points
are packed and sent together, optimizing GPU memory allocation. Finally, the
method used in line 9 of Algorithm 1 is the residual resampling, introduced in
[27].

Figure 3: Visualization of the AUV PF filter in a SVGP map

##  VI Results

This section presents and analyzes the results from the experiments described
above on the tasks of terrain modelling and AUV PF localization.

###  VI-A Comparison of SVGP Mapping Performance

Three instances of the training set have been generated differing on the noise
level applied to corrupt the DR. Only the heading of the vehicle has been
degraded, with the values used seen on Table I. The noise model of the
prediction step of the EKF ($W$ in Eq 1) has been adapted proportionally to
the yaw noise for the three cases. Due to the stochastic nature of both the
dataset generation and the SVGP training, five corrupted datasets have been
created per DR noise level in order to collect statistically representative
results.

Fig 2 (b-d) show a disrupted bathymetric point cloud from the DR collection,
the corresponding SVGP trained with UIs and the consistency error map between
its posterior and the full GT bathymetry. The average map consistency errors
across the tests are summarized in Table I. The RMSE values for different DR
disruption levels indicate very similar reconstruction errors for both DI and
UI SVGPs, with these errors being very close to the original DR point cloud
errors. Interestingly, this suggests that the use of UIs does not influence
greatly the mapping capacity of SVGPs. However, the prediction error, i.e.
error in the area where we removed the original training data, shows a
significant difference in performance between training methods. Indeed the
SVGP UIs predicted terrain is consistently more accurate for both levels of DR
disruption and almost equivalent for no disruption. As a concession, an
increase in the average training iterations when working with UIs can be seen
in the last column. This is an expected effect of the MC sampling of the UIs
and it amounts to rises of 6.26$\%$ and 3,72$\%$ respectively, which does not
hinder the application of UI SVGPs when compared to DI SVGPs. Finally, an
expected increase in the SVGP variance (monitored through the kernel trace) is
reported between SVGPs with DI and UIs. Nonetheless, when analysing the
distribution of the sampled posterior variance over the mapped area, it can be
seen that the variance values in the predicted area are in average lower for
UIs. See Fig. 4 for an example.

Figure 4: DI SVGPs (left) and UI SVGPs (right) posterior variances sampled
over a grid.

###  VI-B SVGP-PF on Predicted Terrain

An AUV localisation mission has been carried out over the second section of
the survey, corresponding to unmmaped terrain. The PF has been initialized
approximately 20 m away from the real AUV pose in order to test if it could
correct a DR drift equivalent to that accumulated through the first section of
the survey. Fig 5 shows an instance of the experiment. The Hugin has followed
a straight trajectory of approximately 500 m (in black) over the unseen area.
The DR trajectory is depicted in red and the PF result in green, with the
final PF variance in blue.

For the experiments, a PF has been fine-tuned to achieve the best performance
across the 5 UI maps generated on the previous section with yaw noise
$1e^{-3}rad/s$. The mission setup and filter parameters have been kept the
same and the measurement noise model has been again fine-tuned for best
performance on DI maps. The average results can be seen in Table II. The
average RMS errors in position and heading across the full trajectories of the
DR and the PFs are compared separately, as advised in [28]. It can be seen
that indeed the PF localisation in the UI maps yields a smaller error in the
pose estimate with no significant difference in the heading.

Figure 5: GT (black), DR (red) and PF (green) trajectories for one round of
AUV localisation over unseen terrain

The PF experiments were run with 60 particles at an average AUV speed of 2
knots per second (faster than the real AUV velocity in the mission) with
additive Gaussian noise in the yaw with covariance $1e-^{4}$, on a laptop with
an Intel processor i7-7700HQ CPU @ 2.80GHz and an NVIDIA GeForce GTX 1060 GPU.

TABLE II: PF localisation and DR errors against GT AUV localization | $RMSE_{xyz}$ | $RMSE_{\theta}$
---|---|---
DR | 12.323 | 0.070
DI SVGP-PF | 11.882 | 0.064
UI SVGP-PF | 10.711 | 0.064

##  VII Conclusions

We have presented a framework to learn SVGP maps with UIs, which harnesses the
flexibility of Monte Carlo integration together with the capability of SVGP to
handle very large datasets. When compared to the current state of the art
method being used in autonomous mapping under localization uncertainty [3, 29,
4], this approach is not constrained to Gaussian UIs and it does not require
to manually set any parameter for the UIs uncertainty approximation.
Furthermore it is more straight-forward to implement, as it does not require
modifying the GP formulation, and it scales better with the dimensionality of
the input data.

In order to test this method in a real bathymetric surveying scenario with an
AUV and a MBES sensor, we have implemented tools to propagate the AUV DR and
sensor uncertainties to the UIs in an efficient manner based on [18]. We have
analyzed the gains and disadvantages of learning bathymetric maps with SVGP
with DIs and UIs in both the tasks of terrain reconstruction and prediction.
For different levels of noise in the dataset, the resulting UI SVGP maps
consistently outperform the DI SVGP ones in the task of terrain prediction
while showing similar results in terrain reconstruction. The concession is
with regards to the training time, with a small increase in the number of
required learning steps with UI SVGPs, as expected. Furthermore, although the
resulting UI SVGP posterior variances are higher, an empirical analysis has
shown that their values on the predicted area are on average lower, which
correlates with the smaller prediction error. Finally we present favourable
results of the SVGP mapping technique with UIs against DIs in a real AUV
localisation scenario with a PF in an entirely predicted area. These reinforce
the conclusion that the UI SVGP has produced not only a more accurate
prediction of the terrain but also of the variance in the area, which the PF
has utilized to provide a more accurate pose estimate.

Overall, the results of the work presented, although limited to one dataset,
show the promising capabilities of the proposed method to apply UI SVGPs in
the tasks of mapping, terrain prediction and vehicle localisation over DI
SVGPs. Further research shall examine how this approach generalises to
datasets containing different variations of the seabed topography and their
full potential for PF localisation. To facilitate this, we release our
development setup in ROS for uncertainty propagation, SVGP mapping and AUV PF
localisation.

## Acknowledgement

The authors thank the Alice Wallenberg foundation for funding MUST, Mobile
Underwater System Tools, project that provided the Hugin AUV. This work was
supported by Stiftelsen fr Strategisk Forskning (SSF) through the Swedish
Maritime Robotics Centre (SMaRC) (IRC15-0046).

## References

  * [1] Ignacio Torroba, Christopher Iliffe Sprague, Nils Bore, and John Folkesson.  Pointnetkl: Deep inference for gicp covariance estimation in bathymetric slam.  IEEE Robotics and Automation Letters, 5(3):4078–4085, 2020.
  * [2] Agathe Girard.  Approximate methods for propagation of uncertainty with Gaussian process models.  University of Glasgow (United Kingdom), 2004.
  * [3] Simon T O’Callaghan and Fabio T Ramos.  Gaussian process occupancy maps.  The International Journal of Robotics Research, 31(1):42–62, 2012\.
  * [4] Marija Popović, Teresa Vidal-Calleja, Jen Jen Chung, Juan Nieto, and Roland Siegwart.  Informative path planning for active field mapping under localization uncertainty.  In 2020 IEEE International Conference on Robotics and Automation (ICRA), pages 10751–10757. IEEE, 2020.
  * [5] Maani Ghaffari Jadidi, Jaime Valls Miro, and Gamini Dissanayake.  Sampling-based incremental information gathering with applications to robotic exploration and environmental monitoring.  The International Journal of Robotics Research, 38(6):658–685, 2019\.
  * [6] James Hensman, Nicolo Fusi, and Neil D Lawrence.  Gaussian processes for big data.  arXiv preprint arXiv:1309.6835, 2013.
  * [7] Matthias Seeger.  Gaussian processes for machine learning.  International journal of neural systems, 14(02):69–106, 2004.
  * [8] Christopher Williams and Matthias Seeger.  Using the nyström method to speed up kernel machines.  In Proceedings of the 14th annual conference on neural information processing systems, number CONF, pages 682–688, 2001.
  * [9] Lehel Csató and Manfred Opper.  Sparse on-line gaussian processes.  Neural computation, 14(3):641–668, 2002.
  * [10] Michalis Titsias.  Variational learning of inducing variables in sparse gaussian processes.  In Artificial intelligence and statistics, pages 567–574. PMLR, 2009.
  * [11] Quoc V Le, Alex J Smola, and Stéphane Canu.  Heteroscedastic gaussian process regression.  In Proceedings of the 22nd international conference on Machine learning, pages 489–496, 2005.
  * [12] Andrew McHutchon and Carl Rasmussen.  Gaussian process training with input noise.  Advances in Neural Information Processing Systems, 24, 2011.
  * [13] Ransalu Senanayake, Simon O’Callaghan, and Fabio Ramos.  Learning highly dynamic environments with stochastic variational inference.  In 2017 IEEE International Conference on Robotics and Automation (ICRA), pages 2532–2539. IEEE, 2017.
  * [14] Weizhe Chen and Lantao Liu.  Long-term autonomous ocean monitoring with streaming samples.  In OCEANS 2019 MTS/IEEE SEATTLE, pages 1–8. IEEE, 2019.
  * [15] Jonathan Ko and Dieter Fox.  Gp-bayesfilters: Bayesian filtering using gaussian process prediction and observation models.  Autonomous Robots, 27(1):75–90, 2009.
  * [16] Ling Zhou, Xianghong Cheng, and Yixian Zhu.  Terrain aided navigation for autonomous underwater vehicles with coarse maps.  Measurement Science and Technology, 27(9):095002, 2016.
  * [17] Simon Julier and Jeffrey K Uhlmann.  A general method for approximating nonlinear transformations of probability distributions.  1996\.
  * [18] Timothy D Barfoot and Paul T Furgale.  Associating uncertainty with three-dimensional poses for use in estimation problems.  IEEE Transactions on Robotics, 30(3):679–693, 2014.
  * [19] Jacob R Gardner, Geoff Pleiss, David Bindel, Kilian Q Weinberger, and Andrew Gordon Wilson.  Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration.  arXiv preprint arXiv:1809.11165, 2018.
  * [20] Carl Edward Rasmussen.  Gaussian processes in machine learning.  In Summer school on machine learning, pages 63–71. Springer, 2003\.
  * [21] Christopher Iliffe Sprague and Petter Ögren.  Learning how to learn bathymetry.  In 2020 IEEE/OES Autonomous Underwater Vehicles Symposium (AUV)(50043), pages 1–2. IEEE.
  * [22] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [23] Maximilian Balandat, Brian Karrer, Daniel R Jiang, Samuel Daulton, Benjamin Letham, Andrew Gordon Wilson, and Eytan Bakshy.  Botorch: Programmable bayesian optimization in pytorch.  arxiv e-prints, pages arXiv–1910, 2019.
  * [24] Théo Galy-Fajou and Manfred Opper.  Adaptive inducing points selection for gaussian processes.  arXiv preprint arXiv:2107.10066, 2021.
  * [25] Chris Roman and Hanumant Singh.  Consistency based error evaluation for deep sea bathymetric mapping with robotic vehicles.  In Proceedings 2006 IEEE International Conference on Robotics and Automation, 2006. ICRA 2006., pages 3568–3574. Ieee, 2006.
  * [26] Nils Bore, Ignacio Torroba, and John Folkesson.  Sparse gaussian process slam, storage and filtering for auv multibeam bathymetry.  In 2018 IEEE/OES Autonomous Underwater Vehicle Workshop (AUV), pages 1–6. IEEE, 2018.
  * [27] Jun S Liu and Rong Chen.  Sequential monte carlo methods for dynamic systems.  Journal of the American statistical association, 93(443):1032–1044, 1998.
  * [28] Edwin Olson and Michael Kaess.  Evaluating the performance of map optimization algorithms.  In RSS Workshop on Good Experimental Methodology in Robotics, volume 15, 2009.
  * [29] Maani Ghaffari Jadidi, Jaime Valls Miro, and Gamini Dissanayake.  Mutual information-based exploration on continuous occupancy maps.  In 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 6086–6092. IEEE, 2015.


#  Confinement of many-body Bethe strings

Jiahao Yang  Tsung-Dao Lee Institute, Shanghai Jiao Tong University, Shanghai
201210, China  Tao Xie  Neutron Scattering Division, Oak Ridge National
Laboratory, Oak Ridge, TN 37831, USA  S. E. Nikitin<EMAIL_ADDRESS>Quantum Criticality and Dynamics Group, Paul Scherrer Institut, CH-5232
Villigen-PSI, Switzerland  Jianda Wu<EMAIL_ADDRESS>Tsung-Dao Lee
Institute, Shanghai Jiao Tong University, Shanghai 201210, China  School of
Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China
Shanghai Branch, Hefei National Laboratory, Shanghai 201315, China  A.
Podlesnyak  Neutron Scattering Division, Oak Ridge National Laboratory, Oak
Ridge, TN 37831, USA

###### Abstract

Based on Bethe-ansatz approach and inelastic neutron scattering experiments,
we reveal evolution of confinement of many-body Bethe strings in ordered
regions of quasi-one-dimensional antiferromagnet $\rm YbAlO_{3}$. In the
antiferromagnetic phase, the spin dynamics is dominated by the confined
length-1 Bethe strings, whose dominancy in the high-energy branch of the
excitation spectrum yields to the confined length-2 Bethe strings when the
material is tuned to the spin-density-wave phase. In the thermal-induced
disordered region, the confinement effect disappears, and the system restores
the conventional quantum integrable physics of the one-dimensional Heisenberg
model. Our results establish a unified picture based on Bethe string for the
spin dynamics in different magnetic phases of $\rm YbAlO_{3}$, and thus
provide profound insight into the many-body quantum magnetism.

Introduction. The one-dimensional (1D) spin-1/2 Heisenberg model, a
paradigmatic model for studying quantum many-body physics, exhibits rich
magnetic excitations such as magnon Bloch (1930); Holstein and Primakoff
(1940); Karabach  _et al._ (1997), spinon Jimbo and Miwa (1995); Caux  _et
al._ (2008); Castillo (2020); Caux and Hagemans (2006); Mourigal  _et al._
(2013), (anti)psinon Karbach  _et al._ (2000), and Bethe strings Bethe (1931);
Takahashi (1971); Gaudin (1971); Takahashi and Suzuki (1972); Takahashi
(1999); Yang  _et al._ (2019); Wang  _et al._ (2018, 2019); Bera  _et al._
(2020); Kohno (2009). Though the former three types of excitations have been
well studied via both theory and experiments, the Bethe strings, exotic many-
magnon bound states, are long-sought in real materials. Recently, with the aid
of Bethe-ansatz calculation Yang  _et al._ (2019), substantial progress has
been made by THz spectroscopy and inelastic neutron scattering (INS)
experiments on the quasi-1D Heisenberg-Ising antiferromagnet $\rm
SrCo_{2}V_{2}O_{8}$ (SCVO) Wang  _et al._ (2018); Bera  _et al._ (2020) and
$\rm BaCo_{2}V_{2}O_{8}$ (BCVO) Wang  _et al._ (2019). The progress reveals
the existence of the Bethe strings in SCVO and BCVO and identifies their vital
contributions to spin dynamics.

Recently, a rare-earth-based material, $\rm YbAlO_{3}$ (YAO), is evidenced to
be a quasi-1D Heisenberg antiferromagnet Wu  _et al._ (2019a); Podlesnyak  _et
al._ (2021). The magnetic ion $\rm Yb^{3+}$ carries on an effective spin
$S=1/2$ due to the combined effect between the spin-orbit coupling and the
crystal field effect. And among $\rm Yb^{3+}$ ions, the dominant spin exchange
interaction ($J$) is isotropic along the chain direction (the $c$ direction in
Fig. 1a) Wu  _et al._ (2019a); Nikitin  _et al._ (2021). The Tomonaga-
Luttinger liquid (TLL) theory has been successfully applied to understand
ground state and low-energy excitations of the material but failed to explain
the physics for the observed spectrum with rich structure beyond low-energy
region Wu  _et al._ (2019a); Nikitin  _et al._ (2021). Thus a unified physical
picture beyond TLL is desired to give a complete understanding for the
complete spectral response.

In this letter by the combined efforts of the INS experiment and the Bethe-
ansatz calculations, we report that the Bethe strings can provide a unified
physical picture to give a full understanding for all spin dynamical spectra
observed in different regions of YAO. Specifically, in ordered regions, the
confined length-1 strings dominate the dynamic spectrum in the
antiferromagnetic (AFM) phase, which gives way to the confined length-2
strings for the high energy branch of the spin dynamic when the material turns
into the spin-density-wave (SDW) phase. In the thermal-induced disordered
region Bethe strings are released from confinement. Our work not only uncovers
the existence and confinement of Bethe strings but also provides a unified
picture based on Bethe string to quantitatively describe the spin dynamics in
different phases of YAO, which is beyond the conventional static order-
parameter paradigm Fan  _et al._ (2020) and low-energy effective theories Wu
_et al._ (2019a); Nikitin  _et al._ (2021); Tomonaga (1950); Luttinger (1960);
Luther and Peschel (1975); Starykh and Balents (2014); Essler  _et al._
(1997).

Figure 1: Schematic illustration of the magnetic structure and phase diagram
of $\rm YbAlO_{3}$. a Below Néel temperature $T_{N}$, the magnetic moments of
Yb ions are antiferromagnetically (AFM) ordered along chain direction (the
$c$-axis), and the interchain coupling is ferromagnetic (FM) in the ab-plane.
b Field-temperature phase diagram with AFM phase, spin-density-wave (SDW)
phase, and presumably transverse-antiferromagnetic (TAF) phase, indicated by
blue, orange, and green, respectively. Two gray dotted lines show the
crossover regions between Tomonaga-Luttinger liquid, quantum critical, and
field-induced FM.

Experimental details. The inelastic neutron scattering (INS) experiments were
performed at the time-of-flight Cold Neutron Chopper Spectrometer (CNCS)
Ehlers  _et al._ (2011, 2016) at the Spallation Neutron Source, Oak Ridge
National Laboratory. A YAO crystal used in previous studies Wu  _et al._
(2019a); Nikitin  _et al._ (2021) with a mass of $\sim$ 0.5 g was mounted in
the $(0~{}K~{}L)$ scattering plane. The zero field data were collected in a
standard cryostat equipped with a dilution insert. The 0.4 T dataset was
measured in a vertical 5 T cryomagnet with a dilution refrigerator and field
was applied along the [1 0 0] direction (the easy axis). To obtain optimal
coverage in the energy and momentum space, the sample was rotated along the
vertical axis by 90∘. The data were collected with an incident neutron energy
$E_{\mathrm{i}}$ = 1.55 meV, with an energy resolution of 0.05 meV at the
elastic line. We used the software packages MantidPlot Arnold  _et al._ (2014)
and Horace Ewings  _et al._ (2016) for data reduction and analysis. The actual
temperature of the sample was calculated using the detailed balance principle,
$S(\mathbf{Q},E)=e^{-E/k_{\rm B}T}S(\mathbf{Q},-E)$, by comparing the spectral
intensity at positive and negative energy transfer.

Theoretical Models. The phase diagram of YAO (Fig. 1b) is obtained from
specific heat and magnetization measurements Nikitin  _et al._ (2021); Wu  _et
al._ (2019a); Fan  _et al._ (2020), which is qualitatively similar to that of
SCVO Bera  _et al._ (2020) and BCVO Klanjšek  _et al._ (2015). At low
temperatures, the interchain coupling helps to stabilize different long-range
orders (LROs) in the material. By tuning the external magnetic field, the
material can change from the AFM phase to the SDW phase after crossing through
the critical field $B_{c1}=0.32$ T. With the field further increasing to about
0.7 T the system begins to transform from the SDW phase to the presumably
transverse-antiferromagnetic (TAF) phase and eventually is fully polarized
when $B>B_{c2}=1.15$ T. In the disordered region of YAO, the interchain
coupling can be neglected following the chain mean-field treatment, and the
dominant magnetic property in YAO can then be described by the spin-1/2
Heisenberg model with external longitudinal field $H_{z}$ (along $a$ direction
in Fig. 1a),

$H_{0}=J\sum_{i}\mathbf{S}_{i}\cdot\mathbf{S}_{i+1}-H_{z}S_{i}^{z},$ (1)

where the AFM coupling $J=0.21$ meV, and $\mathbf{S}_{i}$ is the spin operator
at site $i$ with spin components $S^{\alpha}_{i}\ (\alpha=x,y,z)$. In the
ordered phases, different LROs are characterized by different ordering
wavevector $Q=(1-m)\pi$ where magnetization density $m$ is the ratio of
magnetization $M_{z}$ to its saturation value $M_{s}$. The LRO can effectively
induce a mean field ${\bf{h}}_{z}=h_{Q}\sum_{i}\cos(Qr_{i})\hat{z}$ which
couples to the spin chain Eq. (1) Wu  _et al._ (2019a); Nikitin  _et al._
(2021); Essler  _et al._ (1997); Starykh and Balents (2014),

$H=H_{0}-h_{Q}\sum_{i}\cos(Qr_{i})S^{z}_{i},$ (2)

where field strength $h_{Q}$ is generally dependent on the strength of LRO.
Apart from the AFM phase, the effective mean field ${\bf{h}}_{z}$ in general
is incommensurate, and becomes commensurate when $2\pi/Q$ is a rational number
(for example, ${\bf{h}}_{z}$ is commensurate in AFM phase with $Q=\pi$). In
this letter, following Bethe ansatz approach we analytically study the spin
dynamics based on the effective Hamiltonians of Eq. (1) and Eq. (2) for
disordered and ordered regions, respectively. The obtained results are then
compared in detail with the INS results, which reveals rich spin dynamics as
discussed below. Note, due to the strong Ising-anisotropy of $g$ tensor in YAO
the main contribution to the INS spectrum comes from the longitudinal spin
component Wu  _et al._ (2019a), which makes interpretation of the spectral
function much easier compared with BCVO and SCVO with more isotropic $g$
factor of the ground state doublet.

Figure 2: Pictorial spin configurations for different quasiparticles of 1D
spin-1/2 Heisenberg model and their band structures. a, Ground state (GS) with
SDW ordering and finite magnetization $M_{z}=N/2-M$, where $M$ is the number
of down spins with respect to the fully up-polarized state. b-d, Excitations
generated from the GS. Psinon-antipsinon pair $\psi\psi^{*}$ (b), two-string
$\chi^{2}$ (c), and three string $\chi^{3}$ (d). Viewing the down spin (red
arrow) and up spin (blue arrow) as the magnon and vacuum respectively, many-
body Bethe string $\chi^{j}$ ($j\geq 2$) contains $j$ bounded magnons and
lives in the corresponding $j$-string band in e. The 1-string band
accommodates unbound magnons, and GS has the lowest total energy in the
1-string band.  Figure 3:  Spin dynamics for YAO. INS spectra without magnetic
fields were collected at 1 K (a1) and 0.05 K (b1). The data in (a1) and (b1)
were integrated within $K=[0.3,0.5]$, $H=[-0.1,0.1]$ r.l.u.. No background
subtraction was applied to the data. The INS spectra at $B=0.4$ T at 0.5 K
(c1) and 0.1 K (d1). The data were integrated within all available
$\mathbf{Q}$ ranges along the orthogonal directions and the background
measured at 3 T was subtracted from both spectra. Correspondingly, thermal
DSFs of $H_{0}$ are determined in a2 and c2, and zero temperature DSFs of $H$
are determined in b2 and d2 where the ordering fields are $h_{\pi}=0.27J$ and
$h_{Q}=0.35J$ respectively.

Dynamical structure factor. Following standard linear response theory Negele
and Orland (1988); Chaikin and Lubensky (1995); Zhu (2005), the thermal
dynamical structure factor (DSF) for spin along longitudinal ($z$) direction
is given by ($k_{B}=1$)

$S^{zz}(T,q,\omega)=\frac{2\pi}{\mathcal{Z}}\sum_{\lambda,\mu}e^{-E_{\lambda}/T}|\langle\lambda|S_{q}^{z}|\mu\rangle|^{2}\delta(\hbar\omega-
E_{\mu}+E_{\lambda}),$ (3)

with partition function $\mathcal{Z}$, transfer momentum $q$ and transfer
energy $\hbar\omega$ between two eigenstates $|\lambda\rangle$ and
$|\mu\rangle$ with eigenenergies $E_{\lambda}$ and $E_{\mu}$. The double
summation goes over all the eigenstates of the effective Hamiltonians of Eq.
(1) or Eq. (2). At zero temperature the system is in the ground state, i.e.
$|\lambda\rangle=|GS\rangle$, and the thermal DSF reduces to a single
summation,

$S^{zz}(q,\omega)=2\pi\sum_{\mu}|\langle
GS|S_{q}^{z}|\mu\rangle|^{2}\delta(\hbar\omega-E_{\mu}+E_{GS}).$ (4)

In the following, we focus on the transfer energy range $\hbar\omega\leq 0.8$
meV in compliance with the experimental data.

Analysis and Discussions. For the effective Hamiltonian $H_{0}$ of YAO, its
excitations can be exactly obtained from the Bethe-ansatz method Franchini
(2017). In general, the excitation can be decomposed into Bethe strings with
different lengths. Typically, a Bethe string $\chi^{j}$ ($j\geq 2$) of length
$j$ contains $j$ bounded magnons, referred to as $j$-string. While if $j=1$,
the 1-string $\chi^{1}$ is just unbound magnon, as shown in Fig. 2. The psinon
($\psi$) and antipsinon ($\psi^{*}$) can be understood as the “particle” and
“hole” excitations from the ground state in the 1-string band, which are
always created in pairs ($\psi\psi^{*}s$). And the psinon-antipsinon (PAP)
pairs can adiabatically connect to fractionalized fermionic spinons at $m=0$
and bosonic magnons at $m\simeq 1$ Karbach  _et al._ (2000). All the
excitations can be regarded as the combination of these quasiparticles, i.e.
$n\chi^{j}n^{\prime}\psi\psi^{*}$ with integers $n,\ n^{\prime}\geq 0$. The
obtained dynamical spectra can produce characteristic spectrum continua due to
the many-body nature of Bethe states. Although the effective Hamiltonian $H$
is different from $H_{0}$ by an effective field $\textbf{h}_{z}$, it is beyond
the integrability and can not be solved exactly. Here we tackle this problem
by the truncated string state space method Sup . Following this method the
thermal DSF is calculated to describe the data collected at 1 K and 0.5 K,
while zero-temperature DSFs are obtained to compare with data at 0.05 K and
0.1 K.

In the disordered region with $B=0$ T (corresponding to $m=0$), the dominant
low-energy excitations are 1-strings. When sample temperature $T_{\rm s}=1\
\text{K}\sim 0.4J$, the thermal fluctuations are non-negligible and needs to
be considered in the calculation. From thermal DSF Eq. (3) it is
straightforward to observe that among thermal sampling states
$|\lambda\rangle$ the spectral weight is exponentially suppressed with
increasing eigenenergy $E_{\lambda}$ due to the Boltzmann factor. As such it
is reasonable to select a cutoff $E_{\lambda}\leq 4T_{\rm s}$ for thermal
sampling states $|\lambda\rangle$ (i.e. Boltzmann factor $\lesssim 1.8\%$).
With $T_{\rm s}=1$ K and $m=0$, the thermal DSF Eq. (3) of $H_{0}$ shows that
1-strings dominate the dynamical spectrum (Fig. 3a2) and exhibits a large
thermal broadening effect compared with zero temperature result (see
Supplemental Material Sup ), and the INS spectrum (Fig. 3a1) agrees well with
the theoretical spectrum (Fig. 3a2). In the AFM phase, the zero temperature
DSF Eq. (4) of $H$ with $h_{\pi}=0.27J$ shows that 1-strings are confined by
the effective staggered field and the corresponding INS spectrum (Fig. 3b1) is
confirmed by the theoretical result (Fig. 3b2). Note that the Bragg peak at
zero energy is present in our neutron scattering data at
$\mathbf{Q}=(0~{}0~{}1)$, but not visible in Fig. 3b2 because of the selected
integration range, $K=[0.3,0.5]$ r.l.u. It is worth noting that spinon
confinement observed in SCVO Wang  _et al._ (2016); Bera  _et al._ (2017),
BCVO Faure  _et al._ (2018) and $\rm CaCu_{2}O_{3}$ Lake  _et al._ (2010) is
also the confinement of 1-strings. However, the spinon picture only works at
$m=0$ and is invalid when $m>0$ Jimbo and Miwa (1995); Caux  _et al._ (2008);
Castillo (2020); Caux and Hagemans (2006), in contrast, the Bethe string
picture can describe excitations for both of $m=0$ and $m>0$. Thus, we shall
continuously discuss our results based on Bethe strings in the following.

In the disordered phase with $B=0.4$ T (corresponding to $m\simeq 12\%$),
2-strings $\chi^{2}$ are also favored in the spin dynamics in addition to
1-strings Sup ; Kohno (2009). When $T_{\rm s}=0.5\ \text{K}\sim 0.2J$, the
thermal fluctuation can not be neglected. Following the same strategy in the
case of $T_{\rm s}=1$ K, thermal DSF Eq. (3) is determined for $H_{0}$ at
$T_{\rm s}$ and $m\simeq 12\%$ with energy cutoff $E_{\lambda}\leq 4T_{\rm
s}$. In addition, with the presence of external field, the dipole interaction
in $ab$ plane would deviate from the “magic” line of vanishing dipole
interaction Wu  _et al._ (2019b), leading to non-negligible 3D fluctuations
that further frustrate the intrachain couplings within the framework of the
chain mean-field theory. From theoretical calculations, we find that at
$B=0.4$ T the obtained DSF with adjusted $J=0.18$ meV (Fig. 3c2) is consistent
with the INS experimental spectrum (Fig. 3c1). We note that the continuum of
2-strings is masked by that of 1-strings Sup  and there is no clear energy gap
between these two continua in Fig. 3c1,c2. In fact, this energy gap can be
tuned by external magnetic field $H_{z}$, Ising-anisotropy of the Heisenberg-
Ising model, and next-nearest-neighbor coupling $J_{2}$ in various DSFs. Kohno
(2009); Yang  _et al._ (2019); Keselman  _et al._ (2020).

In the SDW phase with $B=0.4$ T, the observed INS spectrum (Fig. 3d1) reveals
rich structures of string excitations from low to high energy. With
$h_{Q}=0.35J$, $J=0.18$ meV and $m\simeq 12\%$, the zero temperature DSF Eq.
(4) of $H$ shows that both 1- and 2-strings are confined by the effective
field ${\bf h}_{z}$, which plays the same role as that of the staggered field
for the confinement in the AFM phase. Accordingly, confined 2-strings exhibit
a characteristic M-shaped high-energy continuum that is gapped from the low-
energy counterpart of confined 1-strings (Fig. 3d2) Sup . The appearance of a
“gap” between these two continua reflects the energy cost for the confinement
of the 2-strings compared with the deconfined case (Fig. 3c1,c2). Both the
characteristic gap and the shape of continua can be quantitatively compared
between experimental (Fig. 3d1) and theoretical results (Fig. 3d2), which
explicitly confirms the confinement of the Bethe strings of spin chains in
YAO. In addition, at zero transfer energy, a series of satellite peaks appear
at $nQ\ (n=1,2,\cdots)$ in Fig. 3d2 Sup . This is because the effective field
${\bf h}_{z}$ can connect two Bethe states with momentum difference $\Delta
q=Q$ and stabilizes an SDW-ordered ground state. And this ground state has all
the Bethe states with momentum being integer times of $Q$, which makes
possible transitions at many different momentum-transfer. Therefore, the
satellite peaks only occur at $q=nQ$, which is consistent with the multi-
fermion scattering mechanism Nikitin  _et al._ (2021).

Figure 4:  Comparison of the INS scattering intensity with theoretical
calculations. Energy cuts at constant momentum transfers for different
magnetic fields and temperatures. Red points and blue lines are results from
INS measurements and Bethe-ansatz calculations, respectively. The green and
cyan arrows indicate the positions of confined 1-strings and 2-strings,
respectively. The INS data in light-gray regions ($\lesssim$ 0.1 meV) are
contaminated by elastic line contribution.

The energy cuts at constant momenta resolve the continua and peaks of Bethe
strings in Fig. 4. The broad peaks of the 1-string continuum are well captured
by theoretical calculations with consideration of thermal fluctuation (Fig.
4a1-3). From the AFM to the SDW phase, confined 1-strings (Fig. 4b1-3) are
suppressed by a magnetic field and give way to the confined 2-strings at high-
energy regions (Fig. 4d1-3). While in the disordered region at 0.4 T, the
confinement effect disappears, and the deconfined 2-strings are hidden inside
the 1-string continuum (Fig. 4c1-3). Note that the region with energy
$\gtrsim$ 0.1 meV, the theoretical results quantitatively agree with the
experimental data. Especially, the double peaks in experimental data are
evidently confirmed to be the confined 1- and 2-strings (indicated by green
and cyan arrows in Fig. 4d1-3)

Conclusion. To conclude, we establish a unified physical picture based on
Bethe string to understand spin dynamics in different phases of YAO. In the
AFM phase, excitations of confined 1-strings are predominant in the whole
dynamical spectrum, while confined 2-strings take control of the high-energy
branch of spin dynamics in the SDW phase. The excellent comparisons between
theoretical and experimental results explicitly justify the validity of the
Bethe-string-based picture and further confirm the confinement of many-body
Bethe strings in real material. Following the picture, it allows a uniform
characterization of spin dynamics at different phases of a material, which is
beyond the conventional understanding based on static order parameters and
low-energy effective theories. Our study may also inspire research on non-
integrable magnetic systems, and potentially provides a path toward a more
complete understanding of the many-body quantum magnetism.

Acknowledgements. We acknowledge R. Yu for helpful discussion and R. Jiang for
assistance in plotting Fig.1a. This work at Shanghai Jiao Tong University is
supported by National Natural Science Foundation of China No. 12274288 and the
Innovation Program for Quantum Science and Technology Grant No. 2021ZD0301900
and the Natural Science Foundation of Shanghai with grant No. 20ZR1428400
(J.Y. and J.W.). Work at Oak Ridge National Laboratory (ORNL) was supported by
the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences,
Materials Science and Engineering Division. We thank J. Keum for assistance
with X-ray Laue measurements. S. E. N. acknowledges financial support from
innovation program under Marie Skłodowska-Curie Grant No. 884104. This
research used resources at the Spallation Neutron Source, a DOE Office of
Science User Facility operated by the Oak Ridge National Laboratory. X-ray
Laue measurements were conducted at the Center for Nanophase Materials
Sciences (CNMS) (CNMS2019-R18) at ORNL, which is a DOE Office of Science User
Facility.

## References

  * Bloch (1930) F. Bloch, Zeitschrift für Physik 61, 206 (1930).
  * Holstein and Primakoff (1940) T. Holstein and H. Primakoff, Phys. Rev. 58, 1098 (1940).
  * Karabach  _et al._ (1997) M. Karabach, G. Müller, H. Gould, and J. Tobochnik, Comput. Phys. 11, 36 (1997).
  * Jimbo and Miwa (1995) M. Jimbo and T. Miwa,  _Algebraic Analysis of Solvable Lattice Models_ (American Mathematical Society, 1995).
  * Caux  _et al._ (2008) J.-S. Caux, J. Mossel, and I. P. Castillo, J. Stat. Mech.: Theory Exp. 2008, P08006 (2008).
  * Castillo (2020) I. P. Castillo, arXiv:2005.10729 [cond-mat] (2020).
  * Caux and Hagemans (2006) J.-S. Caux and R. Hagemans, J. Stat. Mech.: Theory Exp. 2006, P12013 (2006).
  * Mourigal  _et al._ (2013) M. Mourigal, M. Enderle, A. Klöpperpieper, J.-S. Caux, A. Stunault, and H. M. Rønnow, Nat. Phys. 9, 435 (2013).
  * Karbach  _et al._ (2000) M. Karbach, K. Hu, and G. Muller, arXiv:cond-mat/0008018 (2000).
  * Bethe (1931) H. Bethe, Zeitschrift für Physik 71, 205 (1931).
  * Takahashi (1971) M. Takahashi, Prog. Theor. Phys. 46, 401 (1971).
  * Gaudin (1971) M. Gaudin, Phys. Rev. Lett. 26, 1301 (1971).
  * Takahashi and Suzuki (1972) M. Takahashi and M. Suzuki, Prog. Theor. Phys. 48, 2187 (1972).
  * Takahashi (1999) M. Takahashi,  _Thermodynamics of One-Dimensional Solvable Models_ (Cambridge University Press, 1999).
  * Yang  _et al._ (2019) W. Yang, J. Wu, S. Xu, Z. Wang, and C. Wu, Phys. Rev. B 100, 184406 (2019).
  * Wang  _et al._ (2018) Z. Wang, J. Wu, W. Yang, A. K. Bera, D. Kamenskyi, A. T. M. N. Islam, S. Xu, J. M. Law, B. Lake, C. Wu, and A. Loidl, Nature 554, 219 (2018).
  * Wang  _et al._ (2019) Z. Wang, M. Schmidt, A. Loidl, J. Wu, H. Zou, W. Yang, C. Dong, Y. Kohama, K. Kindo, D. I. Gorbunov, S. Niesen, O. Breunig, J. Engelmayer, and T. Lorenz, Phys. Rev. Lett. 123, 067202 (2019).
  * Bera  _et al._ (2020) A. K. Bera, J. Wu, W. Yang, R. Bewley, M. Boehm, J. Xu, M. Bartkowiak, O. Prokhnenko, B. Klemke, A. T. M. N. Islam, J. M. Law, Z. Wang, and B. Lake, Nat. Phys. 16, 625 (2020).
  * Kohno (2009) M. Kohno, Phys. Rev. Lett. 102, 037203 (2009).
  * Wu  _et al._ (2019a) L. S. Wu, S. E. Nikitin, Z. Wang, W. Zhu, C. D. Batista, A. M. Tsvelik, A. M. Samarakoon, D. A. Tennant, M. Brando, L. Vasylechko, M. Frontzek, A. T. Savici, G. Sala, G. Ehlers, A. D. Christianson, M. D. Lumsden, and A. Podlesnyak, Nat. Commun. 10, 698 (2019a).
  * Podlesnyak  _et al._ (2021) A. Podlesnyak, S. E. Nikitin, and G. Ehlers, Journal of Physics: Condensed Matter 33, 403001 (2021), publisher: IOP Publishing.
  * Nikitin  _et al._ (2021) S. E. Nikitin, S. Nishimoto, Y. Fan, J. Wu, L. S. Wu, A. S. Sukhanov, M. Brando, N. S. Pavlovskii, J. Xu, L. Vasylechko, R. Yu, and A. Podlesnyak, Nat. Commun. 12, 3599 (2021).
  * Fan  _et al._ (2020) Y. Fan, J. Yang, W. Yu, J. Wu, and R. Yu, Phys. Rev. Research 2, 013345 (2020).
  * Tomonaga (1950) S.-i. Tomonaga, Prog. Theor. Phys. 5, 544 (1950).
  * Luttinger (1960) J. M. Luttinger, Phys. Rev. 119, 1153 (1960).
  * Luther and Peschel (1975) A. Luther and I. Peschel, Phys. Rev. B 12, 3908 (1975).
  * Starykh and Balents (2014) O. A. Starykh and L. Balents, Phys. Rev. B 89, 104407 (2014).
  * Essler  _et al._ (1997) F. H. L. Essler, A. M. Tsvelik, and G. Delfino, Phys. Rev. B 56, 11001 (1997).
  * Ehlers  _et al._ (2011) G. Ehlers, A. Podlesnyak, J. L. Niedziela, E. B. Iverson, and P. E. Sokol, Rev. Sci. Instrum. 82, 085108 (2011).
  * Ehlers  _et al._ (2016) G. Ehlers, A. Podlesnyak, and A. I. Kolesnikov, Rev. Sci. Instrum. 87, 093902 (2016).
  * Arnold  _et al._ (2014) O. Arnold, J. C. Bilheux, J. M. Borreguero, A. Buts, S. I. Campbell, L. Chapon, M. Doucet, N. Draper, R. F. Leal, M. A. Gigg, _et al._ , Nucl. Instrum. Methods Phys. Res. Sect. A 764, 156 (2014).
  * Ewings  _et al._ (2016) R. A. Ewings, A. Buts, M. D. Le, J. van Duijn, I. Bustinduy, and T. G. Perring, Nucl. Instrum. Methods Phys. Res. Sect. A 834, 3132 (2016).
  * Klanjšek  _et al._ (2015) M. Klanjšek, M. Horvatić, S. Krämer, S. Mukhopadhyay, H. Mayaffre, C. Berthier, E. Canévet, B. Grenier, P. Lejay, and E. Orignac, Phys. Rev. B 92, 060408 (2015).
  * Negele and Orland (1988) J. W. Negele and H. Orland,  _Quantum Many-Particle Systems_ (Perseus Books, 1988).
  * Chaikin and Lubensky (1995) P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ (Cambridge University Press, 1995).
  * Zhu (2005) Y. Zhu,  _Modern Techniques for Characterizing Magnetic Materials_ (Springer, Boston, MA, 2005).
  * Franchini (2017) F. Franchini,  _An Introduction to Integrable Techniques for One-Dimensional Quantum Systems_, Vol. 940 (Springer, Cham, 2017).
  * (38) See Supplemental Material, for details of truncated string state space method, and additional calculation results and comparisons.
  * Wang  _et al._ (2016) Z. Wang, J. Wu, S. Xu, W. Yang, C. Wu, A. K. Bera, A. T. M. N. Islam, B. Lake, D. Kamenskyi, P. Gogoi, H. Engelkamp, N. Wang, J. Deisenhofer, and A. Loidl, Phys. Rev. B 94, 125130 (2016).
  * Bera  _et al._ (2017) A. K. Bera, B. Lake, F. H. L. Essler, L. Vanderstraeten, C. Hubig, U. Schollwöck, A. T. M. N. Islam, A. Schneidewind, and D. L. Quintero-Castro, Phys. Rev. B 96, 054423 (2017).
  * Faure  _et al._ (2018) Q. Faure, S. Takayoshi, S. Petit, V. Simonet, S. Raymond, L.-P. Regnault, M. Boehm, J. S. White, M. Månsson, C. Rüegg, P. Lejay, B. Canals, T. Lorenz, S. C. Furuya, T. Giamarchi, and B. Grenier, Nature Physics 14, 716 (2018).
  * Lake  _et al._ (2010) B. Lake, A. M. Tsvelik, S. Notbohm, D. Alan Tennant, T. G. Perring, M. Reehuis, C. Sekar, G. Krabbes, and B. Büchner, Nat. Phys. 6, 50 (2010).
  * Wu  _et al._ (2019b) L. S. Wu, S. E. Nikitin, M. Brando, L. Vasylechko, G. Ehlers, M. Frontzek, A. T. Savici, G. Sala, A. D. Christianson, M. D. Lumsden, and A. Podlesnyak, Phys. Rev. B 99, 195117 (2019b).
  * Keselman  _et al._ (2020) A. Keselman, L. Balents, and O. A. Starykh, Phys. Rev. Lett. 125, 187201 (2020).

in 1,2 See pages \x of main_SM