# Solid-state 3He NMR of the superconducting rubidium endofulleride
Rb3(3He@C60).

Murari Soundararajan<EMAIL_ADDRESS>George R. Bacanu  Francesco
Giustiniano  Mark C Walkey  Gabriela Hoffman  Marina Carravetta  Department of
Chemistry, University of Southampton, Southampton SO17 1BJ, UK  Martin R. Lees
Physics Department, University of Warwick, Coventry CV4 7AL, UK  Richard J.
Whitby  Malcolm H. Levitt  Department of Chemistry, University of Southampton,
Southampton SO17 1BJ, UK

(23 August, 2023)

###### Abstract

A new variant of the superconducting fulleride Rb3C60 is presented, with 3He
atoms encapsulated in the C60 cages. The 3He nuclei act as sensitive NMR
probes embedded in the material. The superconducting and normal states are
characterised by 3He NMR. Evidence is found for co-existing vortex liquid and
vortex solid phases below the superconducting transition temperature. A strong
dependence of the spin-lattice relaxation time constant on spectral frequency
is observed in the superconducting state, as revealed by two-dimensional NMR
utilising an inverse Laplace transform. Surprisingly, this phenomenon
persists, in attenuated form, at temperatures well above the superconducting
transition.

##  I Introduction

Reactions of fullerene $\mathrm{C_{60}}$ with alkali metals give rise to a
wide range of salts called _fullerides_ , which consist of negatively-charged
$\mathrm{C_{60}^{n-}}$ anions ($n\in\\{1\ldots 6\\}$) separated by positive
alkali metal cations [1]. Superconductivity in the alkali metal fullerides was
first observed for $\mathrm{K_{3}C_{60}}$, which has a superconducting
transition temperature ($T_{C}$) of 18 K [2]. The transition temperature
$T_{C}$ increases as the size of the alkali metal cation (or unit cell volume)
increases, reaching 33 K in $\mathrm{Cs_{2}RbC_{60}}$ at ambient pressure [3].
The rubidium fulleride $\mathrm{Rb_{3}C_{60}}$, which is the topic of this
paper, has a superconducting transition temperature of around 30 K [4].

Despite their relatively high transition temperatures $T_{C}$ and upper
critical fields ($H_{C2}$), alkali metal fullerides have generally been
understood as s-wave superconductors following conventional Bardeen-Cooper-
Schrieffer (BCS) theory [5]. More recently, phase diagrams of Cs-containing
fullerides as well as other large volume fullerides over wide ranges of
temperature and pressure have been reported, displaying unconventional
metallic and superconducting states [6, 7]. These reports suggest an important
role played by strong electron correlations while retaining s-wave pairing,
placing fullerides in a category distinct from both conventional BCS
superconductors and high-temperature superconductors such as cuprates [8].

NMR is a valuable technique for studying superconductors both in their normal
and superconducting states. The shape and position of the NMR spectrum and
nuclear spin-lattice relaxation time ($T_{1}$) provide information about
electronic interactions and magnetic field distributions inside the material,
and also as direct probes of superconducting parameters such as the band gap
[9, 10, 11]. Rubidium fulleride $\mathrm{Rb_{3}C_{60}}$ has been extensively
studied by both $\mathrm{{}^{13}C}$ and $\mathrm{{}^{87}Rb}$ NMR [12].
However, the interpretation of the NMR data in terms of the electronic
structure and dynamics is made more complicated by the large chemical shift
anisotropy of the fullerene $\mathrm{{}^{13}C}$ nuclei, and by the electric
quadrupole interactions of the Rb nuclides.

Each molecule of fullerene $\mathrm{C_{60}}$ is a symmetrical carbon cage
enclosing a central cavity with an internal diameter of $\sim 3.4$ Å.
Supramolecular complexes known as _endofullerenes_ may be formed, in which
each fullerene cavity accommodates a single atom or molecule. These complexes
are denoted A@$\mathrm{C_{60}}$, where A is the endohedral species. Noble gas
endofullerenes such as He@$\mathrm{C_{60}}$ and Xe@$\mathrm{C_{60}}$ were
first synthesized in very low yield by high-temperature methods [13]. The
remarkable endofullerene N@$\mathrm{C_{60}}$ has been produced by ion
bombardment of $\mathrm{C_{60}}$ [14]. Recently, the multistep synthetic
procedure known as “molecular surgery” has made it possible to synthesize, in
high yield and high purity, a wide range of atomic and molecular
endofullerenes, incorporating species such as $\mathrm{H_{2}}$, $\mathrm{HD}$,
$\mathrm{H_{2}O}$, $\mathrm{HF}$, $\mathrm{CH_{4}}$, $\mathrm{He}$,
$\mathrm{Ne}$, $\mathrm{Ar}$ and $\mathrm{Kr}$ [15]. These systems have been
extensively studied by NMR, infrared spectroscopy, terahertz spectroscopy, and
inelastic neutron scattering [16, 17, 18, 19, 20, 21, 22, 23]. These studies
have revealed the spatial quantisation of the confined atoms or molecules, the
rich interactions of the quantised translational modes with the molecular
rotations and vibrations, and the spin isomerism of confined symmetrical
species such as $\mathrm{H_{2}}$ and $\mathrm{H_{2}O}$ [24, 25]. “Non-bonded”
internuclear J-couplings have been observed between the nucleus of endohedral
$\mathrm{{}^{3}He}$ and the $\mathrm{{}^{13}C}$ nuclei of the encapsulating
cage in $\mathrm{{}^{3}He@C_{60}}$ [26]. The potential energy function
describing the interactions of endohedral He atoms with the fullerene cages
was elucidated by terahertz spectroscopy and neutron scattering [27].

The availability of endofullerenes by molecular surgery suggests the
possibility of creating _endofullerides_ through the reaction of an
endofullerene with an alkali metal, generating materials in which
endofulleride anions such as $\mathrm{A@C_{60}^{3-}}$ are separated by alkali
metal cations. Endofulleride salts such as $\mathrm{Rb_{3}(H_{2}@C_{60})}$ and
$\mathrm{Rb_{3}(H_{2}O@C_{60})}$ have been synthesised and studied by solid-
state NMR, both in the metallic and superconducting states [28].

Endofullerides with A=$\mathrm{{}^{3}He}$ are particularly attractive for NMR
investigations, since the $\mathrm{{}^{3}He}$ nucleus is spin-1/2, has a large
gyromagnetic ratio, no quadrupole moment, and should lack chemical shift
anisotropy when located at the centre of a symmetrical $\mathrm{C_{60}}$ cage.
The zero natural abundance of $\mathrm{{}^{3}He}$ ensures the absence of
background signals from the sample container, probe parts, etc. which can be a
significant nuisance factor for other nuclides [29]. In this paper we describe
some preliminary solid-state $\mathrm{{}^{3}He}$ NMR observations of the He
endofulleride $\mathrm{Rb_{3}(^{3}He@C_{60})}$, over a range of temperatures
spanning the high-temperature metallic phase and the low-temperature
superconducting phase. We show that the endohedral $\mathrm{{}^{3}He}$ atoms
have a negligible influence on the superconducting properties of
$\mathrm{Rb_{3}C_{60}}$ and obtain some novel insights into superconductivity
in this material using relaxation and spectral information from
$\mathrm{{}^{3}He}$ NMR. The endohedral spin-1/2 $\mathrm{{}^{3}He}$ nuclei
serve as excellent probes of the electronic structure and dynamics, and the
internal magnetic fields of the material, with high NMR sensitivity and zero
background interference.

##  II Experimental Details

###  II.1 Samples

The endofullerene $\mathrm{{}^{3}He@C_{60}}$ used in the synthesis of the
endofulleride $\mathrm{Rb_{3}(^{3}He@C_{60})}$ was prepared as in [30].
$\mathrm{Rb_{3}(^{3}He@C_{60})}$ with $22\%$ of the $\mathrm{C_{60}}$
containing a $\mathrm{{}^{3}He}$ atom (a $22\%$ “filling factor”) was prepared
by a refinement of the method reported in [31] for $\mathrm{Rb_{3}C_{60}}$.
$\mathrm{Rb_{6}C_{60}}$ was prepared as a pure (by X-ray diffraction) free-
flowing powder by reacting sublimed $\mathrm{C_{60}}$ with excess Rb metal (10
equivalents, $400\ ^{\circ}\mathrm{C}$, 2 days) in a glass tube sealed under
vacuum, then removing the excess Rb by treatment in a thermal gradient (tube
furnace at $400\ ^{\circ}\mathrm{C}$, 3 days, condensing end of the tube 3 cm
out of the furnace). In a glove box the $\mathrm{Rb_{6}C_{60}}$ was weighed,
then the calculated stoichiometric amount of $\mathrm{He@C_{60}}$ ($44\%$
filled) needed to produce the ternary salt in its pure form added and well
mixed. After sealing in a glass tube the mixture was heated for 4 days at
$450\ ^{\circ}\mathrm{C}$ affording pure (by X-ray diffraction, see Figure
S.3) $\mathrm{Rb_{3}(^{3}He@C_{60})}$ with a filling factor of $22\%$. Samples
used for XRD were taken from the same synthesis batch as the NMR sample. No
amorphous signatures or spurious crystalline peaks are seen in the XRD,
indicating that impurities are limited to 5-10 % of the sample at most.
Reference $\mathrm{{}^{87}Rb}$ NMR spectra of the sample at room temperature
(see section 3 of the SI) display only the expected
$\mathrm{Rb_{3}(^{3}He@C_{60})}$ peaks, further indicating that no significant
quantity of impurities are present.

40 mg of $\mathrm{Rb_{3}(^{3}He@C_{60})}$ was filled into a 4 mm outer
diameter, 2 mm inner diameter borosilicate glass tube inside a glove box kept
under an inert nitrogen environment. The tube was filled with 0.25 bar of
helium gas ($\mathrm{{}^{4}He}$ to ensure good thermal contact with the
environment) and sealed with a flame. The same sample was used for NMR
measurements and for magnetic susceptibility measurements.

Figure 1: DC susceptibility ($\chi$) of $\mathrm{Rb_{3}(^{3}He@C_{60})}$ at an
externally applied magnetic field $\mu_{0}H\,=2.5\mathrm{\,mT}$. Field cooled
(upper trace) and zero field cooled (lower trace) measurements are plotted on
the same graph. a) $\chi$ versus temperature over the entire temperature sweep
range of 5 K to 35 K. b) $\chi$ versus temperature near 30 K from which a
$T_{C}$ of ($30.7\pm 0.1$) K can be estimated.

###  II.2 Instrumentation

####  II.2.1 Helium-3 NMR

NMR experiments were performed at 14.1 T on a Bruker AS 600 WB magnet equipped
with an AVANCE NEO console and fitted with an Oxford Instruments Spectrostat
flow cryostat. A cryogenic $\mathrm{{}^{3}He}$/X probe designed by NMR Service
GmbH and tuned to the $\mathrm{{}^{3}He}$ Larmor frequency of 457.401426 MHz
was used. A detailed description of the construction and NMR performance of
the probe is given in section 1 of the SI. The experiments used a
$\mathrm{{}^{3}He}$ radiofrequency amplitude corresponding to a nutation
frequency of $\sim$32 kHz at ambient temperature. A nutation frequency of
$\sim$26 kHz was used below 35 K to prevent probe arcing.

####  II.2.2 Magnetic susceptibility

A Quantum Design (QD) Magnetic Property Measurement System (MPMS) SQuID
magnetometer was used to collect DC magnetic susceptibility versus temperature
data between 5 and 300 K in zero-field-cooled warming and field-cooled cooling
modes in an applied magnetic field $\mu_{0}H$ of 2.5 mT.

##  III Results and Discussion

###  III.1 Susceptibility measurements

In order to probe the superconducting transition of
$\mathrm{Rb_{3}(^{3}He@C_{60})}$, the magnetic susceptibility of the sample
was measured as a function of temperature at an external applied magnetic
field $\mu_{0}H$ of 2.5 mT, well below the lower critical field
$\mu_{0}H_{C1}=11.4\mathrm{\,mT}$ of $\mathrm{Rb_{3}C_{60}}$ [32].

Field cooled (FC) and zero field cooled (ZFC) measurements of the DC
susceptibility of the $\mathrm{Rb_{3}(^{3}He@C_{60})}$ sample are presented in
Figure 1 in the temperature range of 5 to 35 K. The data indicates a
superconducting transition at $T_{C}=(30.7\pm 0.1)$ K. This is consistent with
the $T_{C}$ estimates of $\mathrm{Rb_{3}C_{60}}$ reported in the literature,
which lie between 25 and 31 K [32, 4, 33]. Hence these data strongly indicate
that the presence of $\mathrm{{}^{3}He}$ inside the fullerene cages does not
significantly influence the onset of superconductivity.

The shielding fraction extracted from the ZFC curve is about 27% and the
Meissner fraction extracted from the FC curve is about 7%. Low superconducting
fractions such as these are commonly observed in fulleride superconductors.
Politis et al.  [32] analyse the particle size distribution in their samples
of $\mathrm{Rb_{3}C_{60}}$ and conclude that a significant proportion of
particles are of comparable or smaller size than the magnetic field
penetration depth at 0 K. They state that the shielding fraction, even at low
applied magnetic fields, is an underestimate of the superconducting fraction.
The large discrepancy between Meissner and shielding fractions is a known
effect in type-II superconductors due to flux pinning in the mixed state just
below $T_{C}$ [34] and has been analysed in detail for high-temperature
superconductors like cuprates [35]. From the spectroscopic and diffraction
data, which gives no evidence of significant impurities, and the documented
difficulties in the interpretation of observed Meissner fractions for finely
divided powders, we contend that a large fraction of the sample is indeed in
the superconducting phase of $\mathrm{Rb_{3}(^{3}He@C_{60})}$ below 30 K,
despite the low observed Meissner fraction.

###  III.2 NMR spectra

$\mathrm{{}^{3}He}$ NMR measurements on the endofulleride
$\mathrm{Rb_{3}(^{3}He@C_{60})}$ were performed at temperatures between 15 K
and 290 K. The series of spectra is shown in Figure 2, which also includes the
$\mathrm{{}^{3}He}$ NMR spectrum of solid $\mathrm{{}^{3}He@C_{60}}$, obtained
in a separate measurement at 292 K. This signal, set to
$\delta=-6.3\mathrm{\,ppm}$ relative to $\mathrm{{}^{3}He}$ gas dissolved in
1-methylnaphthalene [29], was used as a secondary external reference for the
$\mathrm{Rb_{3}(^{3}He@C_{60})}$ spectra.

At room temperature, the $\mathrm{{}^{3}He}$ resonance of
$\mathrm{Rb_{3}(^{3}He@C_{60})}$ is shifted from that of
$\mathrm{{}^{3}He@C_{60}}$ by +25.9 ppm (i.e. in the direction of less
shielding, or higher absolute frequency). This shift is a superposition of a
chemical shift due to the local electronic environment of the endohedral
$\mathrm{{}^{3}He}$ atoms, and a Knight shift due to hyperfine couplings with
the conduction electrons in the high-temperature metallic phase. For
comparison, the $\mathrm{{}^{13}C}$ nuclei of $\mathrm{Rb_{3}C_{60}}$
experience a shift of $+60\mathrm{\,ppm}$ [12].

Figure 2: $\mathrm{{}^{3}He}$ NMR spectra of $\mathrm{Rb_{3}(^{3}He@C_{60})}$
from 15 K to 290 K presented as a waterfall plot along with a reference
spectrum of $\mathrm{{}^{3}He@C_{60}}$ at 292 K. For clarity, the vertical
scale of each spectrum has been adjusted to give a uniform vertical height
over the series. The superconducting transition is seen as a shift in the mean
peak position of the spectrum at 24 K, plotted in red. The split peak below
$T_{C}$ is attributed to the coexistence of vortex solid and vortex liquid
phases (see text).

NMR spectra below 24 K are broadened and shifted to low frequency (high
shielding), relative to those obtained at higher temperatures, as seen in
Figure 2. The shift and broadening is attributed to the diamagnetism of the
superconducting state and the inhomogeneous magnetic field associated with the
Abrisokov vortex lattice [36, 11]. The spectra below 21 K display a split
peak. The low frequency (more shielded) peak appears at $\sim$-5 ppm at 21 K,
and broadens and shifts to lower frequency (higher shielding) as the
temperature goes down, reaching -34 ppm at 15 K. A similar effect has been
reported in the 17O NMR of cuprate superconductors and has been interpreted in
terms of a vortex liquid state coexisting with a frozen vortex lattice [37].
The observed increase in the relative intensity of the lower frequency (high
shielding, right-hand) peak as the temperature is reduced below $T_{C}$ is
consistent with an increased fraction of the frozen vortex phase. Such vortex
phase coexistence in $\mathrm{Rb_{3}C_{60}}$ has been observed through
$\mathrm{{}^{87}Rb}$ and $\mathrm{{}^{13}C}$ transverse relaxation ($T_{2}$)
measurements [38] but not as a splitting in the NMR spectrum, likely due to
the poor resolution and the large $\mathrm{{}^{87}Rb}$ and $\mathrm{{}^{13}C}$
linewidths.

The $\mathrm{{}^{3}He}$ shift, evaluated as the shift of the centre of mass of
the NMR spectrum relative to the $\mathrm{{}^{3}He@C_{60}}$ spectrum, is shown
for temperatures between 15 K and 290 K in Figure 3(a). As the temperature is
decreased starting from room temperature, the $\mathrm{{}^{3}He}$ shift
increases monotonically until the superconducting transition is reached at 24
K. Further cooling causes a dramatic decrease in the $\mathrm{{}^{3}He}$
shift.

The temperature-dependence of the $\mathrm{{}^{3}He}$ shift above $T_{C}$ is
unexpected. As explained in Section 1 of the SI, this shift cannot be due to
instrumental effects. The shift is a superposition of the chemical shift and
the Knight shift. The Knight shift is expected to be temperature-independent
for metallic conductors within the Pauli approximation [39], while the
$\mathrm{{}^{3}He}$ chemical shift is also expected to be roughly temperature-
independent. Hence, the total $\mathrm{{}^{3}He}$ shift is expected to be
almost temperature-independent in the normal metallic state. However, this is
not the observation. A similar temperature-dependence of the shift in the
metallic state has been observed for $\mathrm{{}^{13}C}$ and for
$\mathrm{{}^{87}Rb}$ in several alkali metal fullerides [40, 12].

The superconducting transition temperature of 24 K is lower than the 30 K
indicated by the susceptibility measurements (Figure 1). This is presumably
due to the much higher magnetic field used in the NMR measurements (14.1 T).
The depression of $T_{C}$ by a strong applied field is a common feature for
all superconductors [41].

Figure 3: $\mathrm{{}^{3}He}$ NMR shifts and intensities of
$\mathrm{Rb_{3}(^{3}He@C_{60})}$ as a function of temperature. (a) Estimated
shifts as a function of temperature. The shifts are calculated as the relative
position of the centre-of-mass of the spectrum in
$\mathrm{Rb_{3}(^{3}He@C_{60})}$ relative to that of
$\mathrm{{}^{3}He@C_{60}}$. (b) Integrated intensities of the NMR spectra as a
function of temperature. The expected Curie law dependence of the
$\mathrm{{}^{3}He}$ polarisation at 14.1 T is shown in red.

The integrated signal intensity is presented as a function of inverse
temperature in Figure 3(b), alongside the expected functional form of the
$\mathrm{{}^{3}He}$ Boltzmann polarisation at 14.1 T, which is proportional to
$\tanh{(\hbar\omega_{0}/k_{B}T)}$. Over the relevant range of temperature,
this corresponds closely to the Curie law (polarisation proportional to
$T^{-1}$). The experimental signal intensities follow the Curie law closely
above $T_{C}\simeq 24\mathrm{\,K}$ but deviate sharply below $T_{C}$ becoming
almost temperature-independent at very low temperatures. The reason for this
apparent deviation from the Curie law is not known at this point. Since the
spectra below $T_{C}$ were acquired with 4 scans and a fixed recycling delay
of 300 s, one possibility is that the increase in $T_{1}$ with decreasing $T$
prevents the nuclear spin system from reaching complete thermal equilibrium
with the lattice at the lowest temperatures, within the timeframe of the
experiment. Another possibility is incomplete excitation of the NMR spectrum
as the linewidth increases considerably at 21 K and below. A related effect is
that as the temperature goes down, an increasing fraction of the sample
experiences magnetic fields that are so inhomogeneous that the
$\mathrm{{}^{3}He}$ nuclei become essentially NMR-silent.

###  III.3 Spin-lattice relaxation

Figure 4(a) shows the results of $T_{1}$ measurements performed using the
inversion-recovery method at temperatures ranging from 15 K to 290 K. The
recovery curves are accurately monoexponential above 21 K. However,
significant deviations from single-exponential magnetisation recovery
trajectories appear at temperatures of 21 K and below. Furthermore, as
discussed later, the $T_{1}$ values exhibit strong spectral inhomogeneity at
low temperatures.

Figure 4: Temperature dependence of the $\mathrm{{}^{3}He}$ $T_{1}$. (a) Black
circles: Experimentally determined relaxation rate constants ($T_{1}^{-1}$)
and their error bars. The data points below 24 K are fit to an Arrhenius
function (red line), yielding an energy gap $2\Delta$ of $\left(109\pm
6\right)$ K. The data points above 24 K are fit to a straight line (dashed
blue line) corresponding to a best fit Korringa product of
$(T_{1}T)^{-1}=0.0024\,\mathrm{(s\,K)^{-1}}$. (b) Korringa product
$(T_{1}T)^{-1}$ as a function of temperature.

The $T_{1}$ values plotted in Figure 4(a) are simply equal to the
monoexponential time constants in the high-temperature region. In the low-
temperature region, where the magnetisation recovery is clearly multi-
exponential, the data was fit to a biexponential recovery function. In this
regime, the plotted values of $T_{1}$ are an average of the two biexponential
time constants, weighted by their coefficients in the biexponential fit.

In the normal state, under the Pauli approximation and in the absence of
strong electron correlations, the $T_{1}$ at a temperature $T$ is related to
the Knight shift $K$ by the Korringa relation [39]:

$K^{2}T_{1}T=S,$ (1)

where the Korringa constant $S$ is given by

$S=\frac{\mu_{B}^{2}}{\pi\hbar k_{B}\gamma_{n}^{2}}$

and $\mu_{B}$ is the Bohr magneton, $\hbar$ is the reduced Planck constant,
$k_{B}$ is the Boltzmann constant and $\gamma_{n}$ is the gyromagnetic ratio
of the nucleus. For $\mathrm{{}^{3}He}$, the value of $S$ is $\sim 4.54\times
10^{-7}\mathrm{\,s\,K}$.

If the Knight shift is independent of temperature, as expected for a normal
metal, the Korringa product $(T_{1}T)^{-1}$ is temperature-independent. This
hypothesis is tested in Figure 4(b), which shows the Korringa product as a
function of temperature for the $\mathrm{{}^{3}He}$ resonance of
$\mathrm{Rb_{3}(^{3}He@C_{60})}$. The parameter $(T_{1}T)^{-1}$ is indeed
almost temperature-independent in the normal state, increasing only slightly
with increasing temperature, from a value of 0.0022 $\mathrm{(s\,K)^{-1}}$ at
27 K to 0.0025 $\mathrm{(s\,K)^{-1}}$ at 290 K. Similar deviations from the
Korringa law are often observed in fullerides [40].

In the superconducting state, for a superconducting energy gap $2\Delta$, the
relaxation rate can be modelled using an Arrhenius law [42] as

$T_{1}^{-1}\propto e^{-\Delta/T}.$

Fitting the relaxation data below $T_{C}$ to a function of this form, we
obtain an energy gap $2\Delta$ of $(109\pm 6)$ K. This yields a
$2\Delta/T_{C}$ of $\mathrm{4.5\pm 0.8}$, which agrees with the value of 4.1
estimated using $\mathrm{{}^{13}C}$ NMR by Tycko et. al. [40].

The peak splitting and multi-exponentiality in the recovery curves below 21 K
was explored further by examining the frequency dependence of the $T_{1}$
across the NMR spectrum at low temperature. The inverse Laplace transform
(ILT) is a numerical analysis technique which decomposes an input signal into
individual exponential decay components. It is often used in NMR to analyse
samples with overlapping spectral features that are unresolvable in the
frequency domain on the basis of their $T_{1}$ and $T_{2}$ relaxation
behaviour [43, 44, 45]. An extension to this approach was proposed in [46] in
which the ILT is performed pointwise across the NMR spectrum to obtain a so-
called relaxation-assisted separation (RAS) NMR spectrum with signal
resolution in the frequency and relaxation domains. RAS NMR spectra of our
sample at 15 K, 60 K and 290 K, with the noise reduced by averaging over a
moving window of 51 points in the frequency domain, are presented as contour
maps below. Additional RAS NMR spectra above $T_{C}$, at 30 K, 35 K and 40 K,
are presented in the supplementary information (Figures S.6-S.8).

Figure 5: $\mathrm{{}^{3}He}$ RAS NMR spectrum [46] at 15 K. The image is
presented as a contour map. The roughly diagonal ridge indicates a strong, and
continuous, correlation between the value of $T_{1}$ and the
$\mathrm{{}^{3}He}$ resonance frequency. Note the logarithmic scale used for
$T_{1}$. The 15 K NMR spectrum, extracted from the last slice of the inversion
recovery data used as the input dataset, is displayed above the contour map.
Figure 6: $\mathrm{{}^{3}He}$ RAS NMR spectrum [46] at 290 K. The 290 K NMR
spectrum, extracted from the last slice of the inversion recovery data used as
the input dataset, is displayed above the contour map. The relaxation time
constant $T_{1}$ does not display a dependence on spectral frequency at this
temperature. Figure 7: $\mathrm{{}^{3}He}$ RAS NMR spectrum [46] at 60 K. The
60 K NMR spectrum, extracted from the last slice of the inversion recovery
data used as the input dataset, is displayed above the contour map. The skewed
shape of the peak indicates that the relaxation time constant $T_{1}$ displays
a significant dependence on spectral frequency, even at a temperature more
than 30 K above $T_{C}$.

Figure 5 shows that at 15 K the value of $T_{1}$ varies continously across the
NMR spectrum from 30 s at 53 ppm to 110 s at -29 ppm. As the magnetic field
distribution inside type-II superconductors at externally applied fields in
the $H_{C1}<H_{0}<H_{C2}$ range is highly inhomogeneous due to the vortex
lattice [36], this frequency-dependent relaxation directly maps to a spatially
dependent relaxation. Such inhomogeneous $T_{1}$ behaviour is commonly
observed in cuprates and other d-wave superconductors [47, 11], where it acts
as a probe of vortex excitations. In s-wave superconductors, this is explained
in literature as fast relaxation within the vortex cores that spreads to the
bulk through spin diffusion, and has been observed in $\mathrm{{}^{31}P}$ NMR
on skutterudites [48].

The spectral frequency-dependence of $T_{1}$ is not seen at 290 K (Figure 6),
but begins to develop at 60 K, well above $T_{C}$ (Figure 7). Indeed, the
$T_{1}$ map at 60 K is qualitatively similar to that observed in the
superconducting state at 15 K, albeit with a much smaller range of $T_{1}$
values.

Pennington et al. [49] reported a superficially similar variation in normal-
state $\mathrm{{}^{13}C}$ $T_{1}$ across the spectrum, which they attributed
to a dipolar contribution to the hyperfine interactions for
$\mathrm{{}^{13}C}$, causing a correlated anisotropy in the spectral shift and
the spin-lattice relaxation. However, in our case, the $T_{1}$ values are
_inversely_ correlated with frequency, contrary to the results reported in
ref. [49]. Furthermore, the dipole-dipole mechanism proposed in [49] for the
case of $\mathrm{{}^{13}C}$ appears to be unavailable for the
$\mathrm{{}^{3}He}$ nuclei, which are located at the centres of the fulleride
cages.

The mechanism of the frequency-dependent $T_{1}$ is not known at the current
time. Nevertheless, the detection of an unusual effect that is very strong
below $T_{C}$ but which persists in an attenuated form over a wide temperature
range above $T_{C}$ is potentially significant. The phenomenon is reminiscent
of superconducting fluctuations at temperatures much higher than $T_{C}$,
recently detected by transport measurements on $\mathrm{K_{3}C_{60}}$ [50].
Determining whether these phenomena have a common origin is a matter for
future research.

##  IV Conclusion

A novel “endofulleride” material $\mathrm{Rb_{3}(^{3}He@C_{60})}$ was
synthesised, and the structure and bulk magnetic properties of the parent
fulleride were shown to be essentially unperturbed by the presence of the
endohedral $\mathrm{{}^{3}He}$ atom. The $\mathrm{{}^{3}He}$ nucleus was used
as an NMR probe in both the normal and the superconducting states of the
material, and minor deviations from ideal Knight shift and Korringa behaviour
were observed in the normal state, consistent with prior literature on
$\mathrm{Rb_{3}C_{60}}$. Superconducting state measurements display features
in the $\mathrm{{}^{3}He}$ NMR spectrum which are consistent with known
features of vortex dynamics, such as the apparent co-existence of a vortex
solid and vortex liquid phase. Similar behaviour has been observed by
$\mathrm{{}^{13}C}$ and $\mathrm{{}^{87}Rb}$ relaxation and 2D exchange
experiments in $\mathrm{Rb_{3}C_{60}}$ [38]. However, to our knowledge, this
is the first time such phenomena have been directly observed in the NMR
spectra of fulleride superconductors. The $T_{1}$ relaxation times in the
superconducting state vary strongly across the NMR spectrum, a phenomenon that
has been previously observed in d-wave and some s-wave superconductors but not
in fullerides. Surprisingly, this was seen to persist even above $T_{C}$.
While recent reports suggest the possibility of superconducting fluctuations
above $T_{C}$ in $\mathrm{K_{3}C_{60}}$ [50], more detailed study is needed to
conclusively determine the origin of the spectral distribution of $T_{1}$
values in $\mathrm{Rb_{3}(^{3}He@C_{60})}$.

Endohedral $\mathrm{{}^{3}He}$ has proven to be a valuable and sensitive
addition to the existing NMR tools for studying the properties of fullerides.
While $\mathrm{Rb_{3}C_{60}}$ is at the limit of conventionally understood
superconductivity in fullerides, materials like $\mathrm{Cs_{3}C_{60}}$ and
$\mathrm{Rb_{x}Cs_{3-x}C_{60}}$ are known to exhibit more unusual properties
including strong electron correlations. Introducing endohedral
$\mathrm{{}^{3}He}$ into such materials raises interesting possibilities for
their NMR characterisation.

## Acknowledgments

The authors are grateful to Dr. Richard Bounds, Dr. Karel Kouřil, and Dr.
Roland Thoma for their prior work in the group on similar endofulleride
systems, and to Dr. Mark E. Light for assistance with the X-ray diffraction
measurements. The authors would like to thank Prof. Kosmas Prassides for
valuable discussions and insight into alkali metal fulleride systems.

## Declarations

### Funding

This work was supported by the Engineering and Physical Sciences Research
Council (UK), grant numbers EP/T004320/1, EP/P009980/1, EP/M001962/1,
EP/K00509X/1 and EP/P030491/1. The powder diffraction facilities (Rigaku
Smartlab) were supported through grant number EP/K009877/1.

### Authors’ contributions

MS and GRB performed the NMR experiments and data analysis. MCW and GH
synthesised the fullerene precursor and FG synthesised the fulleride material.
FG performed the X-ray diffraction and MRL performed the magnetometry
measurements. MHL, RJW and MC were involved in conceptualising and acquiring
funding for the project and supervised the research. MS wrote the first draft.
MS, GRB, MRL, RJW and MHL reviewed and edited the manuscript.

## References

  * [1] Holczer, K., Klein, O., Huang, S.-M., Kaner, R.B., Fu, K.-J., Whetten, R.L., Diederich, F.: Alkali-Fulleride Superconductors: Synthesis, Composition, and Diamagnetic Shielding. Science 252(5009), 1154–1157 (1991). https://doi.org/10.1126/science.252.5009.1154
  * [2] Hebard, A.F., Rosseinsky, M.J., Haddon, R.C., Murphy, D.W., Glarum, S.H., Palstra, T.T.M., Ramirez, A.P., Kortan, A.R.: Superconductivity at 18 K in potassium-doped C60. Nature 350(6319), 600–601 (1991). https://doi.org/10.1038/350600a0
  * [3] Gunnarsson, O.: Superconductivity in fullerides. Reviews of Modern Physics 69(2), 575–606 (1997). https://doi.org/10.1103/RevModPhys.69.575
  * [4] Rosseinsky, M.J., Ramirez, A.P., Glarum, S.H., Murphy, D.W., Haddon, R.C., Hebard, A.F., Palstra, T.T.M., Kortan, A.R., Zahurak, S.M., Makhija, A.V.: Superconductivity at 28 K in RbxC60. Physical Review Letters 66(21), 2830–2832 (1991). https://doi.org/10.1103/PhysRevLett.66.2830
  * [5] Pennington, C.H., Stenger, V.A.: Nuclear magnetic resonance of C60 and fulleride superconductors. Reviews of Modern Physics 68(3), 855–910 (1996). https://doi.org/10.1103/RevModPhys.68.855
  * [6] Zadik, R.H., Takabayashi, Y., Klupp, G., Colman, R.H., Ganin, A.Y., Potočnik, A., Jeglič, P., Arčon, D., Matus, P., Kamarás, K., Kasahara, Y., Iwasa, Y., Fitch, A.N., Ohishi, Y., Garbarino, G., Kato, K., Rosseinsky, M.J., Prassides, K.: Optimized unconventional superconductivity in a molecular Jahn-Teller metal. Science Advances 1(3), 1500059 (2015). https://doi.org/10.1126/sciadv.1500059
  * [7] Takabayashi, Y., Ganin, A.Y., Jeglič, P., Arčon, D., Takano, T., Iwasa, Y., Ohishi, Y., Takata, M., Takeshita, N., Prassides, K., Rosseinsky, M.J.: The Disorder-Free Non-BCS Superconductor Cs3C60 Emerges from an Antiferromagnetic Insulator Parent State. Science 323(5921), 1585–1590 (2009). https://doi.org/10.1126/science.1169163
  * [8] Takabayashi, Y., Prassides, K.: Unconventional high-_T_ c superconductivity in fullerides. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 374(2076), 20150320 (2016). https://doi.org/10.1098/rsta.2015.0320
  * [9] MacLaughlin, D.E.: Magnetic Resonance in the Superconducting State. Solid State Physics, vol. 31, pp. 1–69. Academic Press (1976). https://doi.org/10.1016/S0081-1947(08)60541-X
  * [10] F. Smith, D., P. Slichter, C.: The Study of Mechanisms of Superconductivity by NMR Relaxation. In: Beig, R., Beiglböck, W., Domcke, W., Englert, B.-G., Frisch, U., Hänggi, P., Hasinger, G., Hepp, K., Hillebrandt, W., Imboden, D., Jaffe, R.L., Lipowsky, R., Löhneysen, H.V., Ojima, I., Sornette, D., Theisen, S., Weise, W., Wess, J., Zittartz, J., Dolinšek, J., Vilfan, M., Žumer, S. (eds.) Novel NMR and EPR Techniques vol. 684, pp. 243–295. Springer Berlin Heidelberg, Berlin, Heidelberg (2006). https://doi.org/10.1007/3-540-32627-8_9
  * [11] Mounce, A.M., Oh, S., Halperin, W.P.: Nuclear magnetic resonance studies of vortices in high temperature superconductors. Frontiers of Physics 6(4), 450–462 (2011). https://doi.org/%****␣main.bbl␣Line␣275␣****10.1007/s11467-011-0237-5
  * [12] Zimmer, G., Helmle, M., Mehring, M., Rachdi, F., Reichenbach, J., Firlej, L., Bernier, P.: Analysis of 87Rb and 13C Hyperfine Interaction in Rb3C60. Europhysics Letters (EPL) 24(1), 59–64 (1993). https://doi.org/10.1209/0295-5075/24/1/010
  * [13] Saunders, M., Jimenez-Vazquez, H.A., Cross, R.J., Mroczkowski, S., Gross, M.L., Giblin, D.E., Poreda, R.J.: Incorporation of helium, neon, argon, krypton, and xenon into fullerenes using high pressure. Journal of the American Chemical Society 116(5), 2193–2194 (1994). https://doi.org/10.1021/ja00084a089
  * [14] Almeida Murphy, T., Pawlik, Th., Weidinger, A., Höhne, M., Alcala, R., Spaeth, J.-M.: Observation of Atomlike Nitrogen in Nitrogen-Implanted Solid C60. Physical Review Letters 77(6), 1075–1078 (1996). https://doi.org/10.1103/PhysRevLett.77.1075
  * [15] Bloodworth, S., Whitby, R.J.: Synthesis of endohedral fullerenes by molecular surgery. Communications Chemistry 5(1), 121 (2022). https://doi.org/10.1038/s42004-022-00738-9
  * [16] Aouane, M., Armstrong, J., Walkey, M., Hoffman, G., Bacanu, G.R., Whitby, R.J., Levitt, M.H., Rols, S.: A combined inelastic neutron scattering and simulation study of the 3He@C60 endofullerene. Physical Chemistry Chemical Physics (2023). https://doi.org/10.1039/D3CP02253F
  * [17] Bacanu, G.R., Hoffman, G., Amponsah, M., Concistrè, M., Whitby, R.J., Levitt, M.H.: Fine structure in the solution state 13C-NMR spectrum of C60 and its endofullerene derivatives. Physical Chemistry Chemical Physics 22(21), 11850–11860 (2020). https://doi.org/10.1039/D0CP01282C
  * [18] Bloodworth, S., Sitinova, G., Alom, S., Vidal, S., Bacanu, G.R., Elliott, S.J., Light, M.E., Herniman, J.M., Langley, G.J., Levitt, M.H., Whitby, R.J.: First Synthesis and Characterization of CH4@C60. Angewandte Chemie International Edition 58(15), 5038–5043 (2019). https://doi.org/10.1002/anie.201900983
  * [19] Hoffman, G., Bacanu, G.R., Marsden, E.S., Walkey, M.C., Sabba, M., Bloodworth, S., Tizzard, G.J., Levitt, M.H., Whitby, R.J.: Synthesis and 83Kr NMR spectroscopy of Kr@C60. Chemical Communications 58(80), 11284–11287 (2022). https://doi.org/10.1039/D2CC03398D
  * [20] Jafari, T., Shugai, A., Nagel, U., Bacanu, G.R., Aouane, M., Jiménez-Ruiz, M., Rols, S., Bloodworth, S., Walkey, M., Hoffman, G., Whitby, R.J., Levitt, M.H., Rõõm, T.: Ne, Ar, and Kr oscillators in the molecular cavity of fullerene C60. The Journal of Chemical Physics 158(23), 234305 (2023). https://doi.org/10.1063/5.0152628
  * [21] Jafari, T., Bacanu, G.R., Shugai, A., Nagel, U., Walkey, M., Hoffman, G., Levitt, M.H., Whitby, R.J., Rõõm, T.: Terahertz spectroscopy of the helium endofullerene He@C60. Physical Chemistry Chemical Physics 24(17), 9943–9952 (2022). https://doi.org/10.1039/D2CP00515H
  * [22] Krachmalnicoff, A., Bounds, R., Mamone, S., Alom, S., Concistrè, M., Meier, B., Kouřil, K., Light, M.E., Johnson, M.R., Rols, S., Horsewill, A.J., Shugai, A., Nagel, U., Rõõm, T., Carravetta, M., Levitt, M.H., Whitby, R.J.: The dipolar endofullerene HF@C60. Nature Chemistry 8(10), 953–957 (2016). https://doi.org/10.1038/nchem.2563
  * [23] Shugai, A., Nagel, U., Murata, Y., Li, Y., Mamone, S., Krachmalnicoff, A., Alom, S., Whitby, R.J., Levitt, M.H., Rõõm, T.: Infrared spectroscopy of an endohedral water in fullerene. The Journal of Chemical Physics 154(12), 124311 (2021). https://doi.org/10.1063/5.0047350
  * [24] Turro, N.J., Martí, A.A., Chen, J.Y.-C., Jockusch, S., Lawler, R.G., Ruzzi, M., Sartori, E., Chuang, S.-C., Komatsu, K., Murata, Y.: Demonstration of a Chemical Transformation Inside a Fullerene. The Reversible Conversion of the Allotropes of H2@C60. Journal of the American Chemical Society 130(32), 10506–10507 (2008). https://doi.org/10.1021/ja804311c
  * [25] Meier, B., Mamone, S., Concistrè, M., Alonso-Valdesueiro, J., Krachmalnicoff, A., Whitby, R.J., Levitt, M.H.: Electrical detection of ortho–para conversion in fullerene-encapsulated water. Nature Communications 6(1), 8112 (2015). https://doi.org/10.1038/ncomms9112
  * [26] Bacanu, G.R., Rantaharju, J., Hoffman, G., Walkey, M.C., Bloodworth, S., Concistrè, M., Whitby, R.J., Levitt, M.H.: An Internuclear _J_ -Coupling of 3He Induced by Molecular Confinement. Journal of the American Chemical Society 142(40), 16926–16929 (2020). https://doi.org/10.1021/jacs.0c08586
  * [27] Bacanu, G.R., Jafari, T., Aouane, M., Rantaharju, J., Walkey, M., Hoffman, G., Shugai, A., Nagel, U., Jiménez-Ruiz, M., Horsewill, A.J., Rols, S., Rõõm, T., Whitby, R.J., Levitt, M.H.: Experimental determination of the interaction potential between a helium atom and the interior surface of a C60 fullerene molecule. The Journal of Chemical Physics 155(14), 144302 (2021). https://doi.org/10.1063/5.0066817
  * [28] Bounds, R.: Nuclear magnetic resonance of novel type-II superconductors. PhD thesis, University of Southampton (2016). https://eprints.soton.ac.uk/401112/
  * [29] Saunders, M., Jimenez-Vazquez, H.A., Cross, R.J., Mroczkowsklt, S., Anett, F.A.L.: Probing the interior of fullerenes by 3He NMR spectroscopy of endohedral 3He@C60 and 3He@C70. Nature 367(6460), 256–258 (1994). https://doi.org/10.1038/367256a0
  * [30] Hoffman, G., Walkey, M.C., Gräsvik, J., Bacanu, G.R., Alom, S., Bloodworth, S., Light, M.E., Levitt, M.H., Whitby, R.J.: A Solid-State Intramolecular Wittig Reaction Enables Efficient Synthesis of Endofullerenes Including Ne@C60, 3He@C60, and HD@C60. Angewandte Chemie International Edition 60(16), 8960–8966 (2021). https://doi.org/10.1002/anie.202100817
  * [31] McCauley, J.P., Zhu, Q., Coustel, N., Zhou, O., Vaughan, G., Idziak, S.H.J., Fischer, J.E., Tozer, S.W., Groski, D.M.: Synthesis, structure, and superconducting properties of single-phase Rb3C60. A new, convenient method for the preparation of M3C60 superconductors. Journal of the American Chemical Society 113(22), 8537–8538 (1991). https://doi.org/10.1021/ja00022a060
  * [32] Politis, C., Sokolov, A.I., Buntar, V.: Penetration depth and coherence length in superconducting fullerene Rb3C60. Modern Physics Letters B 06(06), 351–357 (1992). https://doi.org/10.1142/S0217984992000430
  * [33] Louis, J.M., Chouteau, G., Ksari, Y., Martinez, G., Hou, J.G., Crespi, V.H., Xiang, X.D., Vareka, W.A., Briceño, G., Zetl, A., Cohen, M.L.: Superconducting Properties of K3C60 and Rb3C60 Single Crystals in High Fields. Molecular Crystals and Liquid Crystals Science and Technology. Section A. Molecular Crystals and Liquid Crystals 245(1), 333–337 (1994). https://doi.org/10.1080/10587259408051710
  * [34] Dahlke, P., Denning, M.S., Henry, P.F., Rosseinsky, M.J.: Superconductivity in Expanded fcc C 60 3- Fullerides. Journal of the American Chemical Society 122(49), 12352–12361 (2000). https://doi.org/10.1021/ja002861d
  * [35] Tomioka, Y., Naito, M., Kishio, K., Kitazawa, K.: The Meissner and shielding effects of high-temperature oxide superconductors. Physica C: Superconductivity 223(3-4), 347–360 (1994). https://doi.org/10.1016/0921-4534(94)91279-3
  * [36] Brandt, E.H.: Magnetic-field variance in layered superconductors. Physical Review Letters 66(24), 3213–3216 (1991). https://doi.org/10.1103/PhysRevLett.66.3213
  * [37] Reyes, A.P., Tang, X.P., Bachman, H.N., Halperin, W.P., Martindale, J.A., Hammel, P.C.: Vortex melting in polycrystalline YBa2Cu3O7 from 17O NMR. Physical Review B 55(22), 14737–14740 (1997). https://doi.org/10.1103/PhysRevB.55.R14737
  * [38] Zimmer, G., Mehring, M., Rachdi, F., Fischer, J.E.: Vortex dynamics in Rb3C60 observed by 87Rb and 13C NMR. Physical Review B 54(6), 3768–3771 (1996). https://doi.org/10.1103/PhysRevB.54.R3768
  * [39] van der Klink, J.J., Brom, H.B.: NMR in metals, metal particles and metal cluster compounds. Progress in Nuclear Magnetic Resonance Spectroscopy 36(2), 89–201 (2000). https://doi.org/10.1016/S0079-6565(99)00020-5
  * [40] Tycko, R., Dabbagh, G., Rosseinsky, M.J., Murphy, D.W., Ramirez, A.P., Fleming, R.M.: Electronic properties of normal and superconducting alkali fullerides probed by 13C nuclear magnetic resonance. Physical Review Letters 68(12), 1912–1915 (1992). https://doi.org/10.1103/PhysRevLett.68.1912
  * [41] Helfand, E., Werthamer, N.R.: Temperature and Purity Dependence of the Superconducting Critical Field, Hc2. II. Physical Review 147(1), 288–294 (1966). https://doi.org/10.1103/PhysRev.147.288
  * [42] Masuda, Y., Redfield, A.G.: Nuclear Spin-Lattice Relaxation in Superconducting Aluminum. Physical Review 125(1), 159–163 (1962). https://doi.org/10.1103/PhysRev.125.159
  * [43] Lee, J.H., Labadie, C., Springer, C.S., Harbison, G.S.: Two-dimensional inverse Laplace transform NMR: altered relaxation times allow detection of exchange correlation. Journal of the American Chemical Society 115(17), 7761–7764 (1993). https://doi.org/10.1021/ja00070a022
  * [44] Moraes, T.: Transformada Inversa de Laplace Para Análise de Sinais de Ressonância Magnética Nuclear de Baixo Campo. Química Nova 44(8), 1020–1027 (2021). https://doi.org/10.21577/0100-4042.20170751
  * [45] Telkki, V.-V., Urbańczyk, M., Zhivonitko, V.: Ultrafast methods for relaxation and diffusion. Progress in Nuclear Magnetic Resonance Spectroscopy 126–127, 101–120 (2021). https://doi.org/10.1016/j.pnmrs.2021.07.001
  * [46] Lupulescu, A., Kotecha, M., Frydman, L.: Relaxation-Assisted Separation of Chemical Sites in NMR Spectroscopy of Static Solids. Journal of the American Chemical Society 125(11), 3376–3383 (2003). https://doi.org/10.1021/ja021173m
  * [47] Wortis, R., Berlinsky, A.J., Kallin, C.: Spin-lattice relaxation in the mixed state of YBa2Cu3O7-$\delta$ and Doppler-shifted d-wave quasiparticles. Physical Review B 61(18), 12342–12351 (2000). https://doi.org/%****␣main.bbl␣Line␣975␣****10.1103/PhysRevB.61.12342
  * [48] Nakai, Y., Hayashi, Y., Ishida, K., Sugawara, H., Kikuchi, D., Sato, H.: Site-selective NMR study of the vortex state in an s-wave superconductor. Physica B: Condensed Matter 403(5-9), 1109–1111 (2008). https://doi.org/10.1016/j.physb.2007.10.261
  * [49] Pennington, C.H., Stenger, V.A., Recchia, C.H., Hahm, C., Gorny, K., Nandor, V., Buffinger, D.R., Lee, S.M., Ziebarth, R.P.: 13C NMR hyperfine couplings, T1 anisotropy, and Korringa relations in Rb2CsC60 : Search for effects of strong correlation. Physical Review B 53(6), 2967–2970 (1996). https://doi.org/10.1103/PhysRevB.53.R2967
  * [50] Jotzu, G., Meier, G., Cantaluppi, A., Cavalleri, A., Pontiroli, D., Riccò, M., Ardavan, A., Nam, M.-S.: Superconducting Fluctuations Observed Far above Tc in the Isotropic Superconductor K3C60. Physical Review X 13(2), 021008 (2023). https://doi.org/10.1103/PhysRevX.13.021008


# Anticipatory Navigation in Crowds by Probabilistic Prediction of Pedestrian
Future Movements

Weiming Zhi1, Tin Lai1, Lionel Ott1, Fabio Ramos1,2
<EMAIL_ADDRESS>School of Computer Science, The University
of Sydney, Australia2 NVIDIA, USA

###### Abstract

Critical for the coexistence of humans and robots in dynamic environments is
the capability for agents to understand each other’s actions, and anticipate
their movements. This paper presents Stochastic Process Anticipatory
Navigation (SPAN), a framework that enables nonholonomic robots to navigate in
environments with crowds, while anticipating and accounting for the motion
patterns of pedestrians. To this end, we learn a predictive model to predict
continuous-time stochastic processes to model future movement of pedestrians.
Anticipated pedestrian positions are used to conduct chance constrained
collision-checking, and are incorporated into a time-to-collision control
problem. An occupancy map is also integrated to allow for probabilistic
collision-checking with static obstacles. We demonstrate the capability of
SPAN in crowded simulation environments, as well as with a real-world
pedestrian dataset.

##  I Introduction

Robots coexisting with humans often need to drive towards a goal while
actively avoiding collision with pedestrians and obstacles present in the
environment. To safely manoeuvre in crowded environments, humans have been
found to anticipate the movement patterns in the environment and adjust
accordingly [1]. Likewise, we aim to endow robots with the ability to adopt a
probabilistic view of their surroundings, and interact with likely future
trajectories of pedestrians. Static obstacles in the environment are encoded
in an occupancy map, and avoided during local navigation. We present the
Stochastic Process Anticipatory Navigation (SPAN) framework, which generates
local motion trajectories that takes into account predictions of the
anticipated movements of dynamics objects.

SPAN leverages a probabilistic predictive model to predict movements of
pedestrians. Motion prediction for dynamic objects, in particular for human
pedestrians, is an open problem and has typically been studied in isolation
without direct consideration of its use in a broader planning framework. Many
recent methods use learning-based approaches [2, 3, 4, 5] to generalise motion
patterns from collected data. To adequately capture the uncertainty of these
predictions, representations of future motion patterns of pedestrians need to
be probabilistic. Furthermore, collision-checking with future pedestrian
position needs to be done at a relatively high time resolution. To address
these requirements, we represent uncertain future motions as continuous-time
stochastic processes. Our representation captures prediction uncertainty and
can be evaluated at an arbitrary time-resolution. We train a neural network to
predict parameters of the stochastic process, conditional on the observed
motion of pedestrians in the environment. Navigating through an environment
requires online generation of collision-free and kinematically-feasible
trajectories. We combine learned probabilistic anticipations of future motion
and an occupancy map into a non-holonomic control problem formulation. To
account for environment interactions a time-to-collision term [6] is
integrated into the control problem. The formulated problem is non-smooth, and
a derivative-free constrained optimiser is used to efficiently solve the
problem to obtain control actions. A receding horizon approach is taken and
the problem is re-optimised at fixed frequency to continuously adapt to
pedestrian prediction updates.

SPAN is novel in jointly integrating data-driven probabilistic predictions of
pedestrian motions and static obstacles in a control problem formulation. The
technical contributions of this paper include: (1) A continuous stochastic
process representation of pedestrian futures, that is compatible with neural
networks, and allows for fast chance-constrained collision-checking, at
flexible resolutions; (2) formulation of a control problem that utilises
predicted stochastic processes as anticipated future pedestrian positions, to
navigate through crowds. We demonstrate that the optimisation problem can be
solved efficiently to allow the robot to navigate through crowded
environments. SPAN also opens opportunities for advances of neural network-
based motion prediction methods to be utilised for anticipatory navigation.
The stochastic process pedestrian motion is compatible with many other neural
network models as the output component of the network.

Figure 1: A non-holonomic robot drives towards a goal while interacting with
pedestrians and static obstacles, and actively adapts to probabilistic
predictions of surrounding movements.

##  II Related Work

SPAN actively predicts the motion patterns of pedestrians in the surrounding
area and avoids collision while navigating towards a goal. Simple solutions to
motion prediction include constant velocity or acceleration models [7]. Recent
methods of motion prediction utilise machine learning to generalise motion
patterns from observed data [2, 5, 3, 4]. These methods have a particular
focus on capturing the uncertainty in the predicted motion, as well as on
conditioning on a variety of environmental factors, such as the position of
neighbouring agents. Generally, past literature examines the motion prediction
in isolation, without addressing issues arising from integrating the
predictions into robot navigation, where anticipation of pedestrian futures is
valuable.

We seek to generate collision-free motion, accounting for predicted motion
patterns. Planning collision-free paths is well-studied, with methods at
differing levels of locality. At one end of the spectrum, probabilistically-
complete planning methods [8, 9] explore the entire search-space to find a
globally optimal solution. At the other end, local methods [10, 11, 12] aim to
perturb collision-prone motion away from obstacles, trading-off global
optimality guarantees for faster run-time. Our work focuses on the local
aspect, where we aim to quickly obtain a sequence of control commands to
navigate through dynamic environments. Various non-linear controllers [13, 14]
also generate control sequences to avoid collision, but often assume obstacles
are not time-varying.

Local collision-avoidance methods for navigation for dynamic environments have
been explored. Methods such as velocity obstacles (VO) [15] and optimal
reciprocal collision avoidance (ORCA) [16] find a set of feasible velocities
for the robot, considering the current velocities of other agents. Partially
observable Markov decision process (POMDP) solvers have been used in [17, 18]
to control linear acceleration for autonomous vehicles, but require steering
controls from a global planner. Methods based on OCRA that indirectly solve
for control sequences, such as [19], are known to be overly conservative. To
address this issue, an approach to directly optimise in the space of controls
was proposed in [20], where a subgradient-based method utilises a time-to-
collision value [6], under constant velocity assumptions. Similarly, SPAN also
optimises against time-to-collision, but using derivative-free constrained
optimisation solvers. We extend the usage of time-to-collision to
probabilistic predictions of pedestrian positions, provided by learning
models, as well as integrating static obstacles captured by occupancy maps.

##  III Problem Overview

Learned Model Stochastic process pedestrian futures Occupancy map Control
problem Derivative-free solver Next controls Robot states Loop Predict Solve
for Update states

Figure 2: An overview of SPAN.

A general overview of SPAN is provided in fig. 2. (1) First, we train a model
offline to predict probabilistic future pedestrian positions, conditioning on
previous positions; (2) The predictive model and occupancy map are incorporate
in an anticipatory control problem; (3) The control problem is solved online,
in a receding horizon manner, with derivative-free solvers to obtain controls
for the next states. This process is looped, with the learned model is
continuously queried, and the problem continuously formulated and solved to
obtain the next the next states.

###  III-A Problem Formulation

This paper addresses the problem of controlling a robot, to navigate towards a
goal $\bm{g}$, in an environment with moving pedestrians and static obstacles.
The state of the robot at a given time $t$ , is given by
$\mathbf{x}(t)=\begin{smallmatrix}[\,x(t)&y(t)&\theta(t)\,]\end{smallmatrix}^{\top}$,
denoting the robot’s two dimensional spatial position and an additional
angular orientation. We focus our investigation on wheeled robots, with the
robot following velocity-controlled non-holonomic unicycle dynamics, where the
system dynamics is given by:

$\displaystyle\dot{\mathbf{x}}=f(\mathbf{x},\mathbf{u})=[v\cos(\theta),v\sin(\theta),\omega]^{\top},$
$\displaystyle\mathbf{u}=[v,\omega]^{\top},$ (1)

where $\dot{\mathbf{x}}$ denotes the state derivatives. The robot is
controlled via its linear velocity $v$ and angular velocity $\omega$. We
denote the controls as $\mathbf{u}$ and the non-linear dynamics as
$f(\mathbf{x},\mathbf{u})$.

Information about _static obstacles_ in the environment is provided by an
occupancy map, denoted as $f_{m}(\cdot)$, which represents the probability of
being occupied $p(\mathrm{Occupied}|\bm{\hat{x}})\in[0,1]$ at a coordinate
$\bm{\hat{x}}\in\mathbb{R}^{2}$. The set of positions of $N$ pedestrian
obstacles, at a given time $t$ from the current time, is denoted by the set
$\mathcal{O}=\\{\bm{o}^{1}(t),\ldots,\bm{o}^{N}(t)\\}$. We are assumed to have
movement data, allowing for training of a probabilistic model to predict
future positions of pedestrians.

##  IV Probabilistic and Continuous Prediction of Pedestrian Futures

In this section, we introduce a continuous _stochastic process_ representation
of future pedestrian motion, and outline how the representation integrates
into a neural network learning model.

###  IV-A Pedestrian Futures as Stochastic Processes

Representations of future pedestrians positions need to capture uncertainty.
Additionally, the time-resolution of the pedestrian position forecasts needs
to synchronise with the frequency of collision-checks. These factors motivate
the use of a continuous-time stochastic process (SP) for future pedestrian
positions. A SP can be thought of as a distribution over functions.
Furthermore, the continuous nature allows querying of future pedestrian
position at an arbitrary time, rather than at fixed resolution, without
additional on-line interpolation.

We start by considering a deterministic pedestrian motion trajectory, before
extending to a _distribution over motion trajectories_. A trajectory is
modelled by a continuous-time function, given as the weighted sum of $m$ basis
functions. We write the $n^{\text{th}}$ pedestrian $\bm{o}^{n}\in\mathcal{O}$
coordinate at $t$ as:

$\displaystyle\mathbf{o}^{n}(t)={\mathbf{W}^{n}}^{\top}\bm{\Phi}(t),$
$\displaystyle\bm{\Phi}(t)=[\bm{\phi}(t,\mathbf{t}^{\prime}),\bm{\phi}(t,\mathbf{t}^{\prime})],$
(2)

where ${\mathbf{W}^{n}}^{\top}\in\mathbb{R}^{m\times 2}$ is a matrix of
weights that defines the trajectory. $\bm{\Phi}(t)$ contains two vectors of
basis function evaluations, $\bm{\phi}(t,\mathbf{t}^{\prime})$, one for each
coordinate axis.
$\bm{\phi}(t,\mathbf{t}^{\prime})=[\phi(t,t^{\prime}_{1}),\ldots,\phi(t,t^{\prime}_{m})]^{\top}$
is a vector of basis function evaluations at $m$ evenly distanced time points,
$\mathbf{t}^{\prime}\in\mathbb{R}^{m}$. We use the squared exponential basis
function given by $\phi(t,t^{\prime})=\exp(-\gamma\|t-t^{\prime}\|^{2})$,
where $\gamma$ is a length-scale hyper-parameter.

To extend our representation from a single trajectory to a distribution over
future motion trajectories, we assume that the weight matrix is not
deterministic, but random with a matrix normal ($MN$) distribution:

$\displaystyle\mathbf{W}^{n}\sim
MN(\mathbf{M}^{n},\mathbf{U}^{n},\mathbf{V}^{n}).$ (3)

The matrix normal distribution is a generalisation of the normal distribution
to matrix-valued random variables [21], parameters mean matrix
$\mathbf{M}^{n}\in\mathbb{R}^{m\times 2}$ and scale matrices
$\mathbf{U}^{n}\in\mathbb{R}^{m\times m}$ and
$\mathbf{V}\in\mathbb{R}^{2\times 2}$. $\mathbf{M}^{n}$ is analogous to the
mean of a normal distribution, and $\mathbf{U}^{n},\mathbf{V}^{n}$ capture the
covariance. A log-likelihood of the matrix normal distribution is given later
in eq. 6. Predicting a SP requires predicting
$\mathbf{M}^{n},\mathbf{U}^{n},\mathbf{V}^{n}$.

Figure 3: Sample trajectories (red) from a predicted distribution over future
pedestrian motion, conditional on previously observed pedestrian coordinates
(blue).

###  IV-B Learning Pedestrian Motion Stochastic Processes

This subsection outlines how to train a model to predict SP representations of
pedestrian motion. In the following subsection, we denote the weight matrix as
$\mathbf{W}$ and parameters of the matrix normal distribution as
$\mathbf{M},\mathbf{U},\mathbf{V}$, without referring to an obstacle index of
a specific pedestrian. We aim to condition on a short sequence of recently
observed pedestrian positions, up to the present
$\\{\mathbf{\hat{o}}_{1},\ldots,\mathbf{\hat{o}}_{p}\\}$, and predict
$\mathbf{M},\mathbf{U},\mathbf{V}$, which define a SP of the motion
thereafter. A predictive model is trained offline with collected data, then
queried online.

####  IV-B1 From timestamped data to continuous functions

Pedestrian motion data is typically in the form of a discrete sequence of $n$
timestamped coordinates $\\{t_{i},\mathbf{\hat{o}}_{i}\\}_{i=1}^{n}$, where
$t_{i}$ gives time and $\mathbf{\hat{o}}_{i}\in\mathbb{R}^{2}$ gives
coordinates. Timestamped coordinates can be converted into a continuous
function by finding a deterministic $\mathbf{W}$. For a pre-defined basis
function $\bf{\Phi}(\cdot)$, we solve the least squares problem:

${\min_{\mathbf{W}}\Big{\\{}\sum_{i=1}^{n}(\mathbf{\hat{o}}_{i}-\underbrace{\mathbf{W}^{\top}\bm{\Phi}(t_{i})}_{\mathbf{o}(t_{i})})^{2}+\lambda\|\mathrm{vec}(\mathbf{W})\|\Big{\\}}}$
(4)

where $\lambda$ is a regularisation hyper-parameter, and $\mathrm{vec}(\cdot)$
is the vectorise operator. By converting timestamped data to continuous
functions, we can obtain a dataset, $\mathcal{D}$, with $N_{D}$ pairs,
$\\{(\\{\mathbf{\hat{o}}_{1},\ldots,\mathbf{\hat{o}}_{p}\\}_{d},\mathbf{W}_{d})\\}_{d=1}^{N_{D}}$,
where $\\{\mathbf{\hat{o}}_{1},\ldots,\mathbf{\hat{o}}_{p}\\}_{d}$ are
previously observed coordinates, and $\mathbf{W}_{d}$ defines the continuous
trajectory from $\mathbf{\hat{o}}_{p}$ onwards. $\mathbf{W}_{d}$ is obtained
from timestamped coordinate data via eq. 4. This dataset is used to train the
predictive model.

####  IV-B2 Training a Neural Network

We wish to learn a mapping, $F$, from a sequence of observed pedestrian
coordinates, $\\{\mathbf{\hat{o}}_{1},\ldots,\mathbf{\hat{o}}_{p}\\}$, to the
parameters over distribution $\mathbf{W}$, that defines future motion beyond
the sequence. We can use a neural network as,

$\overbrace{(\mathbf{M},\mathbf{U},\mathbf{V})}^{\mathclap{\text{defines SP
over future motion}}}=\underbrace{F}_{\mathclap{\text{neural network function
approximator}}}(\overbrace{\\{\mathbf{\hat{o}}_{1},\ldots,\mathbf{\hat{o}}_{p}\\}}^{\mathclap{\text{past
coordinates}}})$ (5)

and train via the negative log-likelihood [21] of the matrix normal
distribution, over $N_{D}$ training samples of dataset $\mathcal{D}$. The loss
function of the network is given as:

$\displaystyle\mathcal{L}=\\!-\\!\sum_{d=1}^{N_{D}}\\!\log\frac{\exp\\{-\frac{1}{2}\mathrm{tr}[{\mathbf{V}}^{-1}\\!(\mathbf{W}_{d}-\mathbf{M})^{\top}\\!\mathbf{U}(\mathbf{W}_{d}-\mathbf{M})]\\}}{(2\pi)^{m}|\mathbf{V}||\mathbf{U}|^{\frac{m}{2}}},$
(6)

where $\mathrm{tr}(\cdot)$ is the trace operator. In our experiments a simple
3 layer fully-connected neural network with width 100 units, and output layer
giving vectorised $\mathbf{M},\mathbf{U},\mathbf{V}$ is sufficiently flexible
to learn $F$. In this work, we focus on outlining a representation of motion
able to be learned and used for collision-avoidance. More complex network
architectures can be used to learn $F$, including recurrent network-based
models [2].

####  IV-B3 Querying the Predictive Model

After training, the neural network outputs $\mathbf{M},\mathbf{U},\mathbf{V}$,
conditional on observed sequences. $\mathbf{M},\mathbf{U},\mathbf{V}$ define
random weight matrix $\mathbf{W}$, and the SP model via eq. 2. Individual
trajectories can be realised from the SP by sampling $\mathbf{W}$. Figure 3
shows sampled trajectories from a predicted SP.

##  V Anticipatory Navigation with Probabilistic Predictions

###  V-A Time-to-Collision Cost for Control

This subsection introduces a control formulation which makes use of
probabilistic anticipations of future pedestrian positions. Human interactions
have been found to be governed by an inverse power-law relation with the time-
to-collision [22], and can be integrated in robot control [20]. To this end,
we formulate an optimal control problem, with the inverse time-to-collision in
the control cost:

$\displaystyle\min_{\mathbf{u}}\quad$
$\displaystyle\Big{\\{}{\|\bm{\hat{x}}(T)-\bm{g}\|_{2}}+\frac{\kappa}{\tau(\mathbf{x}_{0},\mathbf{u},\bm{o},f_{m})}\Big{\\}}$
(7a) s.t. $\displaystyle\dot{\mathbf{x}}(t)=f(\mathbf{x}(t),\mathbf{u}),\quad
t\in[0,T],$ (7b) $\displaystyle\mathbf{x}(0)=\mathbf{x}_{0},$ (7c)
$\displaystyle\mathbf{u}_{lower}\leq\mathbf{u}\leq\mathbf{u}_{upper},$ (7d)

where $T$ is the time horizon; $\bm{\hat{x}}$, $\bm{g}$, $\mathbf{x}$,
$\mathbf{u}$ are the robot coordinates, goal coordinates, states, and target
controls; $\tau$ is the time-to-collision value; $\kappa$ is a weight hyper-
parameter controlling the collision-averseness; $\bm{o}$ denotes the position
of pedestrians, $f_{m}$ denotes an occupancy map;
$\mathbf{u}_{upper},\mathbf{u}_{lower}$ denote the control limits, and $f$
refers to the system dynamics given in eq. 1. Like [20], we define time-to-
collision as the first time a collision occurs under constant controls over
the time horizon, enabling efficient optimisation.

###  V-B Chance-constrained Collision with Probabilistic Pedestrian Movements

This subsection outlines a method to check for collisions between the robot
and SP representations of pedestrian movement, to evaluate a time-to-collision
value. Collision occurs when the distance between the robot and a pedestrian
is below their collision radii. The set of coordinates where collision occurs,
at time $t$, with the $n^{th}$ pedestrian is given by:

$\mathcal{C}^{n}_{t}=\left\\{\bm{x}\in\mathbb{R}^{2}\;\middle|\;\|\bm{x}-\bm{o}^{n}(t)\|_{2}<r_{\bm{\hat{x}}}+r_{{o}_{n}}\right\\},$
(8)

where $r_{\bm{\hat{x}}}$, $r_{\bm{o}_{n}}$ denote the radii of the robot and
the pedestrian, and $\bm{o}^{n}(t)$ the pedestrian position. The position of
pedestrians at a given time is uncertain, and given by the stochastic process
defined in eq. 2. Therefore, we specify collision at $t$ in a chance
constraint manner, where a collision between the robot and a pedestrian is
assumed to occur if the probability of the robot position, $\bm{\hat{x}}(t)$,
being in set $C^{n}_{t}$ is above threshold $\epsilon$, i.e.
$p(\bm{\hat{x}}(t)\in\mathcal{C}^{n}_{t})>\epsilon$. Time-to-collision is
taken as the first instance that a collision occurs, i.e., when the
probability of collision with pedestrians exceeds $\epsilon$. Thus, for the
chance-constraint formulation, the set of robot motions colliding at $t$ is
given by,

$\displaystyle\mathcal{C}_{t}=\\{\bm{\hat{x}}(t)|\max_{n}\\{p(\bm{\hat{x}}(t)\in\mathcal{C}^{n}_{t})\\}>\epsilon\\}.$
(9)

As the weight parameters, $\mathcal{W}$, of the stochastic process is a matrix
normal distribution, the position distribution at a specified time is a
multivariate Gaussian. By considering the weight matrix distribution of the
$n^{th}$ pedestrian, $\mathbf{W}^{n}$ and eq. 2, we have:

$\bm{o}^{n}(t)\sim\mathcal{N}({\mathbf{M}^{n}}^{\top}\bm{\Phi}(t),\mathbf{V}^{n}\bm{\Phi}^{\top}(t)\mathbf{U}^{n}\bm{\Phi}(t)),$
(10)

where $\mathbf{M}^{n},\mathbf{U}^{n},\mathbf{V}^{n}$ are the parameters of the
distribution of weight matrix $\mathbf{W}^{n}$, and $\bm{\Phi}(t)$ contains
basis function evaluations. The position of the robot $\bm{\hat{x}}(t)$ is
deterministic, thus $\bm{\hat{x}}(t)-\bm{o}^{n}(t)$ is also multi-variate
Gaussian, with a shifted mean. The probability of collision with the $n^{th}$
pedestrian is the integral of a Gaussian over a circle of radius
$r_{\bm{\hat{x}}}+r_{\bm{o}_{n}}$. Let
$\bm{d}^{n}=\bm{\hat{x}}(t)-\bm{o}^{n}(t)$, we have the expressions:

$\displaystyle
p(\bm{\hat{x}}(t)\in\mathcal{C}^{n}_{t})=\int_{\|\bm{d}^{n}\|_{2}<r_{\bm{\hat{x}}}+r_{\bm{o}_{n}}}p(\bm{d}^{n})\mathrm{d}\bm{d}^{n},$
(11) $\displaystyle
p(\bm{d}^{n})=\mathcal{N}(\bm{\hat{x}}(t)-{\mathbf{M}^{n}}^{\top}\bm{\Phi}(t),\mathbf{V}^{n}\bm{\Phi}^{\top}(t)\mathbf{U}^{n}\bm{\Phi}(t)).$
(12)

This expression is intractable to evaluate analytically, and an upper bound to
integrals of this kind is given in [23] by:

$p(\bm{\hat{x}}(t)\in\mathcal{C}^{n}_{t})\leq\frac{1}{2}\bigg{\\{}1+\mathrm{erf}\bigg{(}\frac{r_{\bm{\hat{x}}}+r_{\bm{o}_{n}}-\bm{a}^{\top}\bm{d}^{n}(t)}{\sqrt{2\bm{a}^{\top}\Sigma^{n}(t)\bm{a}}}\bigg{)}\bigg{\\}},$
(13)

where $\bm{a}=\frac{\bm{d}^{n}(t)}{\|\bm{d}^{n}(t)\|_{2}}$,
$\Sigma^{n}(t)=\mathbf{V}^{n}\bm{\Phi}^{\top}(t)\mathbf{U}^{n}\bm{\Phi}(t)\in\mathbb{R}^{2\times
2}$ is the covariance of the Gaussian $\bm{o}^{n}(t)$, and
$\mathrm{erf}(\cdot)$ is the error function. Using the defined eq. 13, we can
efficiently evaluate the collision probability upper-bound between the robot
and each pedestrian, and take the first time probability upper-bound exceeds
$\epsilon$ as the time-to-collision with pedestrians. By eqs. 9 and 13 a
collision occurs at $t$ if:

$\max_{n=1,\ldots,N}\Bigg{\\{}\frac{1}{2}\bigg{[}1+\mathrm{erf}\bigg{(}\frac{r_{\bm{\hat{x}}}+r_{\bm{o}_{n}}-\bm{a}^{\top}\bm{d}^{n}(t)}{\sqrt{2\bm{a}^{\top}\Sigma^{n}(t)\bm{a}}}\bigg{)}\bigg{]}\Bigg{\\}}>\epsilon.$
(14)

That is, collisions arise if the probability of collision with any pedestrian
is above threshold $\epsilon$. The robot-pedestrian time-to-collision,
$\tau_{\bm{o}}$, is the earliest time collision with a pedestrian occurs,
namely $\tau_{\bm{o}}=\min_{t\in[0,T]}(t)$, subject to eq. 14.

###  V-C Chance-constrained Collision with Occupancy Maps

Static obstacles in the environment are often captured via occupancy maps.
Time-to-collision values accounting for occupancy maps can be integrated into
SPAN. An occupancy map can be thought of as a function, $f_{m}(\bm{\hat{x}})$,
mapping from position, $\bm{\hat{x}}\in\mathbb{R}^{2}$, to the probability the
position is occupied, $p(\mathrm{Occupied}|\bm{\hat{x}})$. We ensure the
probability of being occupied at any coordinate on the circle of radius,
$r_{\bm{\hat{x}}}$, around our robot position is below the chance constraint
threshold, $\epsilon$. By considering the occupancy map as an implicit
surface, a collision occurs at time $t$ if,

$\displaystyle\max_{\phi\in[-\pi,\pi)}\bigg{\\{}f_{m}\bigg{(}\bm{\hat{x}}(t)+r_{\bm{\hat{x}}}\begin{bmatrix}\sin(\phi)\\\
\cos(\phi)\end{bmatrix}\bigg{)}\bigg{\\}}>\epsilon.$ (15)

Maximising the independent variable, $\phi$, sweeps a circle around the robot
to check for collisions. As the domain of $\phi$ is relatively limited and
querying from occupancy maps is highly efficient, taking under $10^{-6}$s each
query, a brute force optimisation can be done. The time-to-collision due to
static obstacles is then the first time-step where a collision occurs, namely
$\tau_{m}=\min_{t\in[0,T]}(t)$, subject to eq. 15. We take the overall time-
to-collision value for the optimisation in eq. 7a as the minimum of the time-
to-collision for both static obstacles and pedestrians, specifically
$\tau=\mathrm{minimum}(\tau_{\bm{o}},\tau_{m})$.

###  V-D Solving the Control Problem

We take a receding horizon control approach, continuously solving the control
problem defined in eq. 7a \- eq. 7d, and updating the controls for a single
time-step. Algorithm 1 outlines our framework, bringing together the
predictive model, the formulation and solving of control problem. A predictive
model of obstacle movements is trained offline with motion data. The
predictive model and a static map of the environment are provided to evaluate
the time-to-collision, $\tau$, in the control problem. We take a “single-
shooting” approach and evaluate the integral in the control cost via an Euler
scheme with step-size $\Delta t=0.1s$. At each time step we check for
collision by eqs. 14 and 15, and $\tau$ is the first time-step a collision
occurs. The control problem is non-smooth, so we use the derivative-free
solver COBYLA [24] to solve the optimisation problem, and obtain controls
$\mathbf{u}$. As the control-space is relatively small, the optimisation can
be computed quickly. We can re-solve with several different initial solutions
in the same iteration to refine and evade local minima. Controls are updated
at $10$ Hz.

Input:
$\mathbf{x}_{0},\dot{\mathbf{x}}_{0},\mathbf{g},\mathcal{D},f_{m},\Delta t$

Offline: Train neural network, $F$, to predict with dataset $\mathcal{D}$
(section IV-B)

Initialise: $t\leftarrow 0$; $\mathbf{x}(0)\leftarrow\mathbf{x}_{0}$;
$\dot{\mathbf{x}}(0)\leftarrow\dot{\mathbf{x}}_{0}$

while  _goal not reach_ do

Observe $N$ pedestrians’ coordinates, within a past time window,
$\\{\\{\mathbf{\hat{o}}_{1},\ldots,\mathbf{\hat{o}}_{p}\\}\\}_{i=0}^{N}$;

$\mathcal{O}\leftarrow\emptyset$ ;

// collects predictions

Predict
$\\{\mathbf{M}^{n},\mathbf{U}^{n},\mathbf{V}^{n}\\}_{n=0}^{N}=F(\\{\mathbf{\hat{o}}_{1},\ldots,\mathbf{\hat{o}}_{p}\\}_{i})$;

_ _for  _$n=1$ to $N$// For all peds. _

do

// Construct pedestrian futures as SP

$\bm{o}^{n}(t)={\mathbf{W}^{n}}^{\top}\bm{\Phi}(t)$,
$\mathbf{W}^{n}\sim\mathrm{MN}(\mathbf{M}^{n},\mathbf{U}^{n},\mathbf{V}^{n})$;

$\mathcal{O}\leftarrow\mathcal{O}\cup\\{\bm{o}^{n}(t)\\}$;

end for

// optimise for control with pred.

$\mathbf{u}\leftarrow$ Solve eqs. 7a, 7b, 7c and 7d, with
$\mathcal{O},f_{m},\mathbf{x}(t)$ to get controls;

$\mathbf{x}(t+\Delta t)\leftarrow\mathrm{NextState(\mathbf{x}(t),\mathbf{u})}$
// Execute controls

$t\leftarrow t+\Delta t$;

if  _$\|\bm{\hat{x}}(t)-\mathbf{g}\|_{2} <\epsilon_{goal}$_ then

Terminates;

// goal reached

end if

end while

Algorithm 1 Stochastic Process Anticipatory Navigation

##  VI Experimental Evaluation

We empirically evaluate the ability of SPAN to navigate in environments with
pedestrians and static obstacles. Results and desirable emergent behaviour are
discussed below.

###  VI-A Experimental Setup

Both simulated pedestrian behaviour and real-world pedestrian data are
considered in evaluations. We use a power-law crowd simulator [22, 25], and
real-world indoor pedestrian data from [26]. We construct two simulated
environments, _Sim 1_ and _Sim 2_ , both with 24 pedestrians. _Sim 1_ models a
crowded environment with noticeable congestion; _Sim 2_ models a relatively
open environment with less frequent pedestrian interaction. The simulated
pedestrians have a maximum speed of $1m/s$. Pedestrian trajectories collected
at different times in dataset are overlaid, increasing the crowd density. The
maps used in our experiments are shown in fig. 4. We navigate the robot from
starting point to a goal. Gazebo simulator [27] was used for visualisation.
The following metrics are used:

  1. 1.

Time-to-goal (TTG): The time, in seconds, taken for the robot to reach the
goal;

  2. 2.

Duration of collision (DOC): The total time, in seconds, the robot is in
collision.

We evaluate SPAN against the following methods:

  * •

Closed-loop Rapidly-exploring Random Trees∗ (CL-RRT∗). An asymptotic optimal
version of [28] with a closed-loop prediction and a nonlinear pure-pursuit
controller [29]. [28] demonstrates an ability to handle dynamic traffic, and
is a DARPA challenge entry. At each iteration, a time budget of $0.1s$ is
given to the planner to find a free trajectory towards the goal, with
pedestrians represented as obstacles. If a feasible solution is not found, the
robot will remain at its current pose.

  * •

A reactive controller. The reactive controller continuously solves to find
optimal collision-free controls for the next time-step, without anticipating
future interaction.

In our experiments, we set the time-to-collision weight, $\kappa$=100; the
trajectory length-scale hyper-parameter, $\gamma=0.01$; the control limits on
linear/ angular velocity as $-1\leq\bm{v}\leq 1$ and $0\leq\omega\leq 1$. Each
iteration we solve for a “look-ahead” time horizon of $T=4.0s$. The collision
threshold is set at $\epsilon=0.25$, and the radii for the robot and the
pedestrians are given as $r_{\bm{\hat{x}}},r_{\bm{o}}=0.4m$. In the simulated
setups, the predictive model is trained with additionally generated data from
the crowd simulator; in the dataset setup, $80\%$ of the data was used to
train the predictive model, with the remainder for evaluation. We condition on
observed pedestrian movements in the last 0.5s, to predict the next 4.0s.
Occupancy maps provided are implemented with the method in [30]. Controls were
solved with different initialisations 40 times per iteration. The predictive
model is implemented in Tensorflow [31], control formulation and COBYLA [24]
are in FORTRAN, interfacing with Python2. Experiments run on a laptop with an
Intel i7 CPU, and 32GB RAM.

###  VI-B Experimental Results

Figure 4: Motion of pedestrians (black) and robot (red) overlaid on an
occupancy map of the environments for Sim1 (left), Sim2 (center), and Dataset
(right). Note that crossing trajectories does not necessarily indicate a
collision, as moving agents could be at the same position at different times.
The starting position of the pedestrians and robot are indicated by circles.
The robot navigates smoothly and sensibly, remaining collision-free.

Figure 5: _Sim 1_ setup. (Left): The reactive controller often drives head-on
into the crowd, leading the robot to become stuck and resulting in minor
collisions. A considerable amount of time is needed to resolve the congestion.
(Right): The robot controlled by SPAN evades congestion when possible, and
follows behind pedestrians that are moving in the same direction.

The resulting motion trajectories by the controlled robot in the experiments
are illustrated in fig. 4, where the motion trajectories of the robot is
marked in red, while those of pedestrians are marked in black. At a glance, we
see that the controlled robot navigates smoothly and sensibly. The evaluated
metrics are tabulated in table I. The average computational times of an entire
iteration of SPAN, including prediction and obtaining controls, are 27ms,
31ms, 33ms in _Dataset_ , _Sim 1_ , _Sim 2_ setups respectively, with no
iteration taking more than 50ms in any setup.

SPAN is able to control the robot reach the goal, without collisions, in all
of the experiment environments. In particular, relative to the compared
methods, SPAN demonstrates markedly safer navigation through the _Sim 1_ and
_Dataset_ environment setups. Both of these environments are relatively
cluttered, with regions of congestion arising. We observe that both CL-RRT∗
and the reactive controller has a greater tendency to freeze amongst large
crowds, while SPAN gives more fluid navigation amidst crowds. Figure 5 shows
how the reactive controller (left figure) and SPAN (right figure) navigate the
robot through crowds in the dense _Sim 1_ environment. The reactive controller
often controls the robot to drive towards incoming pedestrians head-on. This
behaviour often results in driving into bottleneck positions and congestion
arising, with the robot then unable to escape from a trapped position. On the
other hand, the robot controlled by SPAN avoids getting stuck in crowds, as
the future position of pedestrians are anticipated. Incoming pedestrians are
actively evaded, while the robot follows other pedestrians through the crowd.

In the _Dataset_ setup, both the compared CL-RRT∗ and the reactive controller
produce very collision-prone navigation trajectories. This is due to the
completely dominate pedestrian behaviour, which follows collected trajectory
data, whereas the simulated pedestrian crowds also attempts to avoid collision
with the robot. Both CL-RRT∗ and the reactive controller give a relatively
straight trajectory towards the goal and are in collision for more than $3s$,
colliding with three different pedestrians. Meanwhile SPAN produces a safe and
meandering trajectory. The trajectories of the controlled robot by CL-RRT∗ and
SPAN is illustrated in fig. 7, along with plots of how the distance to goal
changes in time. We argue, in the _Dataset_ setup, although the time-to-goal
of SPAN is marginally longer, the navigated path is much safer with no
collisions, while the compared methods have more than 3s of collision.
Additionally, all three of the evaluated methods perform similarly in the
comparatively open _Sim 2_ environment, producing collision-free and efficient
trajectories. This is due to the ample space in the environment, with relative
little congestion occurring.

Figure 6: Anticipatory behaviour of our controlled robot can be observed.
(Left, center) The robot follows (red arrow) pedestrians anticipated to move
ahead, while preemptively avoiding incoming pedestrians. (Right) Robot as red
circle and pedestrians as black. Sampled trajectories (in purple, next 4.0s)
visualise predicted SP, and the robot trajectory is in red. The robot evades
an incoming group of pedestrians preemptively, and moves toward pedestrians
predicted to move away.

cl-rrt∗span

Figure 7: (Left, Center) shows time against the distance between the robot and
goal, for the CL-RRT∗ and SPAN respectively. The times that are in collision
are given in red, while the collision-free times are in blue. (Right) shows
the motion of the robot as controlled by CL-RRT∗ in blue, and that as
controlled by ours in green. Time-steps where the robot is in collision are in
red. We see that SPAN is able to provide a safe path for the controlled robot,
while CL-RRT∗ provides an aggressive path, resulting in over $3s$ of
collision.

###  VI-C Emergent Behaviour

The formulated control problem in eqs. 7a, 7b, 7c and 7d does not explicitly
optimise for the following of other agents. However, we observe following
behaviour, behind other pedestrians, also moving in the goal direction,
through crowds. Figure 6 (left, center) show the controlled robot following
moving pedestrians through crowds, while avoiding those moving towards the
robot. This is also detailed in fig. 6 (right), where predicted future
positions are visualised by sampling and overlaying trajectories (in purple)
from the SP model. The entire navigated motion trajectory taken by the robot
is outlined in red. We observe the robot taking a swerve away from the group
of incoming pedestrians preemptively, and turning to follow the other group
predicted to move away. At that snapshot, the pedestrians in the lower half of
the map blocks the path through to the goal, but are all predicted to move
away, giving sufficient space for our robot to follow and pass through. The
emergence of the crowd-following behaviour facilitates smooth robot navigation
through crowds.

Method | Metrics | Dataset | Sim 1 (crowded) | Sim 2 (open)
---|---|---|---|---
SPAN (ours) | TTG | 21.7 | 25.9 | 17.4
| DOC | 0.0 | 0.0 | 0.0
CL-RRT∗ | TTG | 19.0 | 36.8 | 17.8
| DOC | 3.7 | 1.1 | 0.0
Reactive | TTG | 19.1 | 43.5 | 17.1
Controller | DOC | 3.4 | 0.4 | 0.0

TABLE I: The evaluations of SPAN against compared methods. TTG gives the time
(in seconds) taken by the robot to reach the goal, and DOC gives the total
duration (in seconds) the robot is in collision. SPAN finds a collision-free
solution in all three evaluated set-ups.

##  VII Conclusions and Future Work

We introduce Stochastic Process Anticipatory Navigation (SPAN), a framework
for anticipatory navigation of non-holonomic robots through environments
containing crowds and static obstacles. We learn to anticipate the positions
of pedestrians, modelling future positions as continuous stochastic processes.
Then predictions are used to formulate a time-to-collision control problem. We
evaluate SPAN with simulated crowds and real-world pedestrian data. We observe
smooth collision-free navigation through challenging environments. A desired
emergent behaviour of following behind pedestrians, which move in the same
direction, arises. The stochastic process representation of pedestrian motion
is compatible as an output for many neural network models, bringing together
advances in motion prediction literature with robot navigation. Particularly,
conditioning on more environmental factors to produce higher quality
predictions within SPAN, is a clear direction for future work.

## References

  * [1] T. Higuchi, “Visuomotor control of human adaptive locomotion: Understanding the anticipatory nature,” Frontiers in Psychology, 2013.
  * [2] A. Alahi, K. Goel, V. Ramanathan, A. Robicquet, L. Fei-Fei, and S. Savarese, “Social lstm: Human trajectory prediction in crowded spaces,” IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 961–971, 2016\.
  * [3] W. Zhi, L. Ott, and F. Ramos, “Kernel trajectory maps for multi-modal probabilistic motion prediction,” Conference on Robot Learning, 2019\.
  * [4] W. Zhi, R. Senanayake, L. Ott, and F. Ramos, “Spatiotemporal learning of directional uncertainty in urban environments with kernel recurrent mixture density networks,” IEEE Robotics and Automation Letters, 2019\.
  * [5] A. Gupta, J. Johnson, L. Fei-Fei, S. Savarese, and A. Alahi, “Social gan: Socially acceptable trajectories with generative adversarial networks,” in IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2018\.
  * [6] J. C. Hayward, “Near-miss determination through use of a scale of danger,” Highway Research Record, 1972.
  * [7] A. Rudenko, L. Palmieri, M. Herman, K. M. Kitani, D. M. Gavrila, and K. O. Arras, “Human motion trajectory prediction: a survey,” The International Journal of Robotics Research, 2020.
  * [8] S. M. Lavalle, “Rapidly-exploring random trees: A new tool for path planning,” tech. rep., Department of Computer Science, Iowa State University, 1998.
  * [9] T. Lai, F. Ramos, and G. Francis, “Balancing global exploration and local-connectivity exploitation with rapidly-exploring random disjointed-trees,” in 2019 International Conference on Robotics and Automation (ICRA), 2019.
  * [10] O. Khatib, “Real-time obstacle avoidance for manipulators and mobile robots,” IEEE International Conference on Robotics and Automation, 1985.
  * [11] C.-A. Cheng, M. Mukadam, J. Issac, S. Birchfield, D. Fox, B. Boots, and N. Ratliff, “Rmpflow: A computational graph for automatic motion policy generation,” in International Workshop on Algorithmic Foundations of Robotics (WAFR), 2018.
  * [12] N. Ratliff, M. Zucker, J. A. Bagnell, and S. Srinivasa, “Chomp: Gradient optimization techniques for efficient motion planning,” in IEEE International Conference on Robotics and Automation, 2009.
  * [13] M. Neunert, C. de Crousaz, F. Furrer, M. Kamel, F. Farshidian, R. Siegwart, and J. Buchli, “Fast nonlinear model predictive control for unified trajectory optimization and tracking,” 2016 IEEE International Conference on Robotics and Automation (ICRA), 2016.
  * [14] G. Williams, P. Drews, B. Goldfain, J. M. Rehg, and E. Theodorou, “Aggressive driving with model predictive path integral control,” in IEEE International Conference on Robotics and Automation (ICRA), 2016.
  * [15] P. Fiorini and Z. Shiller, “Motion planning in dynamic environments using velocity obstacles,” The International Journal of Robotics Research, 1998\.
  * [16] J. V. D. Berg, S. J. Guy, M. Lin, and D. Manocha, “Reciprocal n-body collision avoidance,” in ISRR, 2009.
  * [17] Y. Luo, P. Cai, A. Bera, D. Hsu, W. S. Lee, and D. Manocha, “Porca: Modeling and planning for autonomous driving among many pedestrians,” IEEE Robotics and Automation Letters, 2018.
  * [18] P. Cai, Y. Luo, D. Hsu, and W. S. Lee, “Hyp-despot: A hybrid parallel algorithm for online planning under uncertainty,” International Journal of Robotics Research, 2020.
  * [19] D. Bareiss and J. van den Berg, “Generalized reciprocal collision avoidance,” The International Journal of Robotics Research, 2015.
  * [20] B. Davis, I. Karamouzas, and S. J. Guy, “Nh-ttc: A gradient-based framework for generalized anticipatory collision avoidance,” Robotics: Science and Systems, 2020.
  * [21] P. Dutilleul, “The mle algorithm for the matrix normal distribution,” Journal of Statistical Computation and Simulation, 1999.
  * [22] I. Karamouzas, B. Skinner, and S. J. Guy, “Universal power law governing pedestrian interactions,” Phys. Rev. Lett., 2014.
  * [23] H. Zhu and J. Alonso-Mora, “Chance-constrained collision avoidance for mavs in dynamic environments,” IEEE Robotics and Automation Letters, 2019.
  * [24] M. Powell, “A direct search optimization method that models the objective and constraint functions by linear interpolation,” Advances in Optimization and Numerical Analysis, 1994.
  * [25] I. Karamouzas, N. Sohre, R. Narain, and S. J. Guy, “Implicit crowds: Optimization integrator for robust crowd simulation,” ACM Transactions on Graphics, vol. 36, July 2017.
  * [26] A. Rudenko, T. P. Kucner, C. S. Swaminathan, R. T. Chadalavada, K. O. Arras, and A. J. Lilienthal, “Thör: Human-robot navigation data collection and accurate motion trajectories dataset,” IEEE Robotics and Automation Letters, 2020.
  * [27] N. Koenig and A. Howard, “Design and use paradigms for gazebo, an open-source multi-robot simulator,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2004.
  * [28] Y. Kuwata, J. Teo, S. Karaman, G. Fiore, E. Frazzoli, and J. How, “Motion planning in complex environments using closed-loop prediction,” in AIAA Guidance, Navigation and Control Conference and Exhibit, p. 7166, 2008.
  * [29] O. Amidi and C. E. Thorpe, “Integrated mobile robot control,” in Mobile Robots V, vol. 1388, pp. 504–523, International Society for Optics and Photonics, 1991.
  * [30] F. Ramos and L. Ott, “Hilbert maps: Scalable continuous occupancy mapping with stochastic gradient descent,” The International Journal of Robotics Research, 2016.
  * [31] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, M. Kudlur, J. Levenberg, R. Monga, S. Moore, D. G. Murray, B. Steiner, P. Tucker, V. Vasudevan, P. Warden, M. Wicke, Y. Yu, and X. Zheng, “Tensorflow: A system for large-scale machine learning,” OSDI’16, 2016.