# Gaussian Process based Stochastic Model Predictive Control for Cooperative
Adaptive Cruise Control

Sahand Mosharafian†, Mahdi Razzaghpour∗, Yaser P. Fallah∗, Javad Mohammadpour
Velni†

†School of Electrical & Computer Engineering, Univ. of Georgia, Athens, GA
∗Dept. of Electrical & Computer Engineering, Univ. of Central Florida,
Orlando, FL

###### Abstract

Cooperative driving relies on communication among vehicles to create
situational awareness. One application of cooperative driving is Cooperative
Adaptive Cruise Control (CACC) that aims at enhancing highway transportation
safety and capacity. Model-based communication (MBC) is a new paradigm with a
flexible content structure for broadcasting joint vehicle-driver predictive
behavioral models. The vehicle’s complex dynamics and diverse driving
behaviors add complexity to the modeling process. Gaussian process (GP) is a
fully data-driven and non-parametric Bayesian modeling approach which can be
used as a modeling component of MBC. The knowledge about the uncertainty is
propagated through predictions by generating local GPs for vehicles and
broadcasting their hyper-parameters as a model to the neighboring vehicles. In
this research study, GP is used to model each vehicle’s speed trajectory,
which allows vehicles to access the future behavior of their preceding vehicle
during communication loss and/or low-rate communication. Besides, to overcome
the safety issues in a vehicle platoon, two operating modes for each vehicle
are considered; free following and emergency braking. This paper presents a
discrete hybrid stochastic model predictive control, which incorporates system
modes as well as uncertainties captured by GP models. The proposed control
design approach finds the optimal vehicle speed trajectory with the goal of
achieving a safe and efficient platoon of vehicles with small inter-vehicle
gap while reducing the reliance of the vehicles on a frequent communication.
Simulation studies demonstrate the efficacy of the proposed controller
considering the aforementioned communication paradigm with low-rate
intermittent communication.

###### Index Terms:

Cooperative adaptive cruise control, Model predictive control, Hybrid
stochastic automata, Non-parametric Bayesian inference, Gaussian process,
Model-based communication

##  I Introduction

Adaptive cruise control (ACC) is a radar-based system, which is designed to
enhance driving comfort and safety by adjusting a vehicle’s speed to match the
speed of the preceding vehicle. However, ACC only has a small impact on the
highway capacity [1]. The objective of cooperation in a highway scenario is to
ensure that all vehicles in a lane move at the same speed while maintaining a
desired formation geometry, which is specified by a desired inter-vehicle gap
policy. The cooperative driving with constant spacing policy, which is called
platooning, makes vehicles maintain a constant distance from their immediate
predecessor while for the cooperative adaptive cruise control (CACC) constant
time headway gap is used, in which the desired following distance should be
proportional to the speed of the vehicle; the higher the speed, the larger the
distance. CACC and platooning have the potential to increase the highway
capacity when they reach a high market penetration [2]. It is shown that
platooning is more sensitive to communication losses than the CACC is, mainly
due to its very close coupling between vehicles [3].

Figure 1: A simple representation of the system model. The distance between
vehicle $i$ and vehicle $i-1$ is denoted by $d_{i}$, and $x_{i}$ is the
location of the rear bumper of the $i^{th}$ vehicle. Vehicle $i$ receives data
from its preceding vehicle through wireless communication.

Model-based communication (MBC) is a recently-explored communication
scalability solution, which has shown a promising potential to reduce the
channel congestion [4]. The fundamental intention behind the MBC scheme is to
utilize a more flexible content structure for the broadcast packets based on
the joint vehicle-driver predictive behavioral models in comparison with the
Basic Safety Message (BSM) content structure defined by J2735 standard. MBC
can potentially shrink the payload size by extracting an abstract
representation of the vehicle’s state. In addition, it reduces the
transmission rate by enabling the recipient vehicles to predict their
neighbors mobility more accurately for a longer future time horizon. As a
result of reduced load, the MBC experiences lower rate of packet collision
compared to its baseline counterpart in different traffic densities [5, 6].
MBC can utilize different methods of modeling vehicle movement behavior. Non-
parametric Bayesian inference techniques, particularly Gaussian Processes
(GPs) are among the promising methods for analytically tractable modeling of
joint vehicle-driver behavior, which at the same time is not limited to some
certain criteria. The driver behavioral models are functions of different
factors such as the driver’s personal driving style, road traffic, weather
condition, etc. Therefore, movement models may become very complex. Gaussian
process regression is a powerful non-parametric tool used to infer values of
an unknown function given previously collected measurements [7]. In this work,
we intend to use the Gaussian process regression to derive the model of the
remote vehicle and its driver as a unique object. In addition to exhibiting
very good generalization properties, a major advantage of GPs is that they
come equipped with a measure of model uncertainty, making them particularly
beneficial for safety critical applications. The GP-based MBC module trains
the GP based on the last received information. This procedure results in
generating a new situational awareness messages which carry the last updated
abstract model of the vehicle’s state.

The main goal of the V2V communications is to enable every vehicle in a
Vehicular Ad-hoc NETwork (VANET) to frequently inform the surrounding vehicles
about its most recent dynamic states. However, its performance degrades with
channel load increase (i.e., high number of surrounding vehicles). Due to
safety reasons, CACC methods should be robust against packet loss and
communication failure [8]. CACC would degrade to ACC if no action is taken to
compensate for packet losses [9]. According to [10], three general solutions
are available to cope with the safety problem during communication failure:

  * •

Improving the communication protocol (e.g., see [11]);

  * •

Increasing the vehicle headway during failure (e.g., see [12]);

  * •

Adapting CACC algorithms to cope with the safety problem (e.g., using a model
predictive control approach).

The method proposed in this work aims at improving the robustness of the
controller to the communication losses, and successfully operating in low-rate
communication scenarios. The method proposed in this paper fuses the first and
third solutions by taking advantage of two important factors, namely, model-
based communication and emergency braking.

In this paper, the CACC design problem consists of two operating modes for
each vehicle: free following and emergency braking. In free following mode,
each vehicle simply follows its preceding vehicle while the vehicle in the
emergency braking mode performs hard braking to avoid possible collision and
preserve safety. These operating modes add binary variables to the system, and
MBC contains stochastic information about preceding vehicle speed profile. To
formulate the underlying problem, a hybrid stochastic model is needed. Hence,
discrete hybrid stochastic automata, introduced initially in [13] suits the
aforementioned stochastic system with both continuous and discrete variables.
Finally, each vehicle employs Model Predictive Control (MPC) to optimize their
speed trajectory considering system constraints. By employing MBC, during
communication loss and low-rate communication, CACC will not degrade to ACC
since the preceding vehicle trajectory can be produced through GP-based MBC,
which assures the efficacy and the safety of the system even with a low
communication rate.

The contributions of the paper are as follows. The paper fuses model-free and
model-based communication paradigms to assure the CACC safety and performance
during communication losses. At each successful communication event, vehicles
directly share their future speed profile, and they also share/update a model
which allows vehicles to predict their preceding vehicle’s behavior so that
each autonomous vehicle is able to plan its speed profile properly. Besides,
integrating two operating modes (namely free following and emergency braking)
improves safety while avoiding unnecessary braking through penalizing
excessive braking.

The rest of the paper is organized as follows. The system model and Gaussian
process are described in Section II. System model considering constraints and
operating modes are presented in mixed logical dynamical from in Section III.
The underlying MPC problem is then explained in Section IV. The performance of
the proposed controller is evaluated through simulation studies in section V,
and concluding remarks are finally made in Section VI.

##  II Preliminaries and System Model

This section first provides a state-space representation of the vehicle model
used for the purpose of CACC design. Then, an introduction to Gaussian
processes is given.

###  II-A System Model

In this paper, we aim at designing controllers for CACC that build a string
stable vehicle platoon with longitudinal movement while reducing the reliance
of the system on frequent communication. For safety reasons, emergency braking
is taken into account to avoid possible collisions that might occur due to the
sudden changes in the preceding vehicle’s speed or random packet loss. The
CACC enables reducing the distance between CAVs in a platoon thereby
increasing the road capacity while preserving safety.

A CACC system with $N_{v}$ vehicles is considered here, and the index
$i\in\\{0,1,\ldots,N_{v}-1\\}$ is used to represent the $i^{th}$ vehicle with
$i=0$ being the leader vehicle. The distance between $i^{th}$ vehicle and its
preceding vehicle at time $t$ is denoted by $d_{i}(t)$ (see Fig. 1) and
defined as

$\displaystyle d_{i}(t)=x_{i-1}(t)-x_{i}(t)-l_{i},$ (1)

where $x_{i}$ and $l_{i}$ are the location of the $i^{th}$ vehicle rear
bumper, and the length of the $i^{th}$ vehicle, respectively. Hence, the
desired spacing policy for vehicle $i$ shown by $d^{*}_{i}(t)$ can be defined
as

$\displaystyle d^{*}_{i}(t)=\tau\,v_{i}(t)+d^{s}_{i},$ (2)

where $v_{i}(t)$ is the vehicle speed at time $t$, $\tau$ is the time gap, and
$d^{s}_{i}$ is the stand still distance. Using constant time headway gap as in
(2) improves string stability and safety[14]. The difference between the
distance and its desired value is defined as

$\displaystyle\Delta d_{i}(t)=d_{i}(t)-d^{*}_{i}(t),$ (3)

and speed difference between the $i^{th}$ vehicle and its preceding vehicle is
considered as $\Delta v_{i}(t)=v_{i-1}(t)-v_{i}(t)$ and therefore
$\Delta\dot{d}_{i}(t)=\Delta v_{i}(t)-\tau\,a_{i}(t)$ and
$\Delta\dot{v}_{i}(t)=a_{i-1}(t)-a_{i}(t)$, where $a_{i}(t)$ is the
acceleration of the $i^{th}$ vehicle at time $t$. Let us consider a linear
model for the $i^{th}$ vehicle as

$\displaystyle\dot{a_{i}}(t)=-\mathnormal{f}_{i}a_{i}(t)+\mathnormal{f}_{i}u_{i}(t),$
(4)

where the constant $\mathnormal{f}_{i}$ represents drive-line dynamics and its
unit is $s^{-1}$. State-space representation of the given system for the
$i^{th}$ vehicle is as follows

$\dot{\mathbf{x}}_{i}(t)=A_{i}\,\mathbf{x}_{i}(t)+B_{i}\,u_{i}(t)+D\,a_{i-1}(t)\\\
\\\ =\begin{bmatrix}0&1&-\tau\\\ 0&0&-1\\\
0&0&-\mathnormal{f}_{i}\end{bmatrix}\mathbf{x}_{i}(t)+\begin{bmatrix}0\\\ 0\\\
\mathnormal{f}_{i}\end{bmatrix}u_{i}(t)+\begin{bmatrix}0\\\ 1\\\
0\end{bmatrix}a_{i-1}(t)$ (5)

where $\mathbf{x}_{i}(t)=[\Delta d_{i}(t)\,\,\Delta
v_{i}(t)\,\,a_{i}(t)]^{T}$. For the leader ($i=0$), $a_{i-1}(t)$ in (5) can be
assumed to be zero.

Using forward-time approximation for the first-order derivative, (5) can be
written in discrete-time form. The discrete-time state space model for each
follower is as follows

$\mathbf{x}_{i}(k+1)=(I+t_{s}\,A_{i})\,\mathbf{x}_{i}(k)+t_{s}\,B_{i}\,u_{i}(k)+t_{s}\,D\,a_{i-1}(k)$
(6)

where $I$ is the identity matrix and $t_{s}$ is the sampling time.

###  II-B Gaussian Process

A Gaussian process is a collection of random variables, any finite number of
which have a joint Gaussian distribution. Gaussian process regression (GPR) is
a non-parametric Bayesian approach to provide uncertainty measurements on the
predictions. A Gaussian process is completely specified by its mean function
and covariance function, which is called the kernel function in the GP context
and defines the trend of target function based on the similarity pattern among
the observed values. We define mean function $m(t)$ and the kernel function
$k(t,t^{\prime})$ of a real process $\mathcal{V}(t)$ as

$m(\mathbf{t})=\mathbb{E}[\mathcal{V}(\mathbf{t})]$ (7)
$k\left(\mathbf{t},\mathbf{t}^{\prime}\right)=\mathbb{E}\left[(\mathcal{V}(\mathbf{t})-m(\mathbf{t}))\left(\mathcal{V}\left(\mathbf{t}^{\prime}\right)-m\left(\mathbf{t}^{\prime}\right)\right)\right]$
(8)

and write the Gaussian process as

$\mathcal{V}(\mathbf{t})\sim\mathcal{G}\mathcal{P}\left(m(\mathbf{t}),k\left(\mathbf{t},\mathbf{t}^{\prime}\right)\right).$
(9)

Different vehicle dynamics can be considered as separate time series which
should be regressed using an appropriate supervised learning method. The
regression problem here is equivalent to inferring the characteristics of the
unknown target functions which have generated these time series. In our case,
the random variables represent the speed of the vehicle at time $t$. We are
interested in incorporating the knowledge that the training data provides
about the function and its future values. The joint distribution of the past
values, $\mathcal{V}$, and the future values $\mathcal{V}^{\ast}$ according to
the prior is

$\left[\begin{array}[]{l}\mathbf{\mathcal{V}}\\\
\mathbf{\mathcal{V}}^{*}\end{array}\right]\sim\mathcal{N}\left(\mathbf{0},\left[\begin{array}[]{ll}K(t,t)&K\left(t,t^{*}\right)\\\
K\left(t^{*},t\right)&K\left(t^{*},t^{*}\right)\end{array}\right]\right).$
(10)

If there are $n$ training points and $n^{\ast}$ test points, then
$K(X,X^{\ast})$ denotes the $n\times n^{\ast}$ covariance matrix evaluated at
all pairs of training and test points, and similarly for the other entries
$K(X,X)$, $K(X^{\ast},X^{\ast})$ and $K(X^{\ast},X)$. To get predictions at
unseen time of interest $t^{\ast}$, the predictive distribution can be
calculated by weighting all possible predictions by their calculated posterior
distribution or in probabilistic terms

$\displaystyle\mathbf{\mathcal{V}}^{*}\mid
t^{*},t,\mathbf{\mathcal{V}}\sim\mathcal{N}$
$\displaystyle\left(K\left(t^{*},t\right)K(t,t)^{-1}\mathbf{\mathcal{V}},\right.$
(11) $\displaystyle K\left(t^{*},t^{*}\right)$
$\displaystyle\left.-K\left(t^{*},t\right)K(t,t)^{-1}K\left(t,t^{*}\right)\right).$

Function values $\mathcal{V}^{\ast}$ (corresponding to test inputs $t^{\ast}$)
can be sampled from the joint posterior distribution by evaluating the mean
and covariance matrix from (11).

We choose the most commonly used kernel in machine learning which is the
Gaussian form radial basis function (RBF) kernel. It is also commonly referred
to as the exponentiated quadratic or squared exponential kernel. The kernel’s
parameters are estimated using the maximum likelihood principle. The advantage
of GP is its capability to substantially increase the forecasting accuracy
over longer prediction horizons without reception of new raw information or
model updates during a certain period. The computational complexity of a GP
regression strongly depends on the number of data points $N$. Training window
size has been set to 5 latest equally spaced (last 0.5 second)
observed/received speed samples in time.

##  III Vehicle Model and Constraints in Mixed Logical Dynamical Form

According to [13], a stochastic system with both binary and
continuous/discrete-time variables and inputs can be modeled using a discrete
hybrid stochastic automata (DHSA) which consists of four components: a
switched affine model, an event generator, a mode selector, and a finite state
machine. Representing a system in DHSA is explained in details in [13]. CACC
problem in this paper consists of two operating modes: free following and
emergency braking. The system is modeled using DHSA in which the
uncontrollable events are defined in terms of uncertainty in speed prediction
using GP regression. Then, DHSA is presented in mixed logical dynamical (MLD)
form [15], and finally, mixed-integer quadratic programming is employed to
find optimal control input(s). In the remainder of this section, system
inequalities in the form of an MLD are derived.

The constraints on the system include bounds on the acceleration, input, road
speed limit, and distance between vehicles (note that a negative distance
implies collision and therefore should not occur). The following inequalities
should always hold true

$\displaystyle a_{i}^{min}\leq a_{i}(k)\leq a_{i}^{max},$ (12a) $\displaystyle
u_{i}^{min}\leq u_{i}(k)\leq u_{i}^{max},$ (12b) $\displaystyle v_{i}(k)\leq
v_{max},$ (12c) $\displaystyle d_{i}(k)>0.$ (12d)

Besides, for passenger comfort, system input changes are bounded as follows

$\begin{gathered}t_{s}\,u_{i}^{min}\leq u(k+1)-u(k)\leq
t_{s}\,u_{i}^{max}.\end{gathered}$ (13)

To determine the mode of the system, the auxiliary binary variable
$\xi_{i}^{e}(k)$ is considered such that

$\displaystyle\xi_{i}^{e}(k)=1\iff\text{Emergency braking constraint is
activated}.$

First, emergency braking constraint needs to be defined. Here, for safety
purpose and collision avoidance, a constraint is considered such that if
$\Delta d_{i}(k)$ goes below a fixed level ($\,\underline{d}\,$), $i^{th}$
vehicle should brake with the minimum possible acceleration
($u_{i}(k)=u_{i}^{min}$), thereby

$\displaystyle\Delta d_{i}(k)+\underline{d}\leq 0\iff\xi_{i}^{e}(k)=1,$ (14)

which can be described in MLD form using the following inequalities

$\begin{gathered}\Delta
d_{i}(k)+\underline{d}\leq\overline{b_{i}^{e}}[1-\xi_{i}^{e}(k)],\\\ \Delta
d_{i}(k)+\underline{d}\geq\varepsilon+\xi_{i}^{e}(k)[\underline{b_{i}^{e}}-\varepsilon],\end{gathered}$
(15)

where $\varepsilon$ is the machine precision, and

$\begin{gathered}\underline{b_{i}^{e}}\leq\Delta
d_{i}(k)+\underline{d}\leq\overline{b_{i}^{e}}.\end{gathered}$

Next, to enforce hard braking, an upper bound constraint on $u_{i}(k)$ is
added such that

$\displaystyle
u_{i}(k)\leq\xi_{i}^{e}(k)\,u_{i}^{min}+[1-\xi_{i}^{e}(k)]\,u_{i}^{max}.$ (16)

Thus, as far as the system operates in the free following mode, the upper
bound on input is $u_{i}^{max}$. However, when emergency braking mode is
activated, the upper bound becomes $u_{i}^{min}$. Since the emergency braking
mode forces the input $u_{i}(k)=u_{i}^{min}$, the speed of the vehicle may
become negative. To handle this issue, the constraint $v_{i}(k)\geq 0$ is
added to the system. However, this may result in infeasibility because of
relation of speed at time $k$ and $k+1$; $v_{i}(k)\geq 0$ while
$v_{i}(k+1)=v_{i}(k)+t_{s}\,a_{i}(k)$ may become negative. To handle this edge
case, two new auxiliary binary variables $\xi_{i}^{v}(k)$ and $\xi_{i}^{E}(k)$
are defined, and (16) rewritten as

$\begin{gathered}u_{i}(k)\leq\,[\xi_{i}^{E}(k)]\,u_{i}^{min}+[1-\xi_{i}^{E}(k)]\,u_{i}^{max},\end{gathered}$
(17)

where

$\displaystyle v(k)<\underline{v}_{i}\iff\xi_{i}^{v}(k)=0,$ (18a)
$\displaystyle\xi_{i}^{E}(k)=\xi_{i}^{e}(k)\,\xi_{i}^{v}(k).$ (18b)

This implies that if the emergency braking event is activated but the vehicle
speed is below the threshold, there is no need to enforce the vehicle input to
get to its minimum. With this assumption and by simply choosing
$\underline{v}_{i}=1\,m/s$, this issue is addressed. Consequently, the vehicle
performs emergency braking if and only if $\xi_{i}^{e}(k)\,\xi_{i}^{v}(k)=1$.
The statement in (18a) in the MLD form turns into

$\begin{gathered}v_{i}(k)-\underline{v}_{i}\leq\overline{b_{i}^{v}}\,\xi_{i}^{v}(k),\\\
v_{i}(k)-\underline{v}_{i}\geq\varepsilon+[1-\xi_{i}^{v}(k)](\underline{b_{i}^{v}}-\varepsilon),\end{gathered}$

where

$\begin{gathered}\underline{b_{i}^{v}}\leq
v_{i}(k)-\underline{v}_{i}\leq\overline{b_{i}^{v}},\end{gathered}$

while (18b) is equivalent to the following inequalities

$\begin{gathered}\xi_{i}^{E}(k)\geq\xi_{i}^{e}(k)+\xi_{i}^{v}(k)-1,\\\
\xi_{i}^{E}(k)\leq\xi_{i}^{e}(k),\\\
\xi_{i}^{E}(k)\leq\xi_{i}^{v}(k).\end{gathered}$ (19)

Whenever a vehicle performs emergency braking, the constraints for passenger
comfort should be ignored. Hence, (13) is updated as follows

$\begin{gathered}(1-\xi_{i}^{E}(k))\,t_{s}\,u_{i}^{min}-\xi_{i}^{E}(k)\,\bar{u_{i}}\leq
u(k+1)-u(k),\\\
u(k+1)-u(k)\leq\xi_{i}^{E}(k)\,\bar{u_{i}}+(1-\xi_{i}^{E}(k))\,t_{s}\,u_{i}^{max},\end{gathered}$
(20)

where $\bar{u_{i}}=u_{i}^{max}-u_{i}^{min}$.

As described earlier, each vehicle needs to access the preceding vehicle’s
future speed trajectory (or the future acceleration profile), which is
available in either of two ways; each vehicle will receive it through
communication, or the vehicle will generate the profile using the last GP
model received from the preceding vehicle. In the case that the speed profile
is generated using GP, the variance of the GP output should be also taken into
account. Hence, the uncertain part of the future velocity (the kernel
introduced in (8)) can be discretized to some levels
$\\{n_{i-1}^{1}(k),n_{i-1}^{2}(k),...,n_{i-1}^{m_{i-1}}(k)\\}$ with known
probabilities $\\{p_{i-1}^{1},p_{i-1}^{2},...,p_{i-1}^{m_{i-1}}\\}$.
Therefore, the preceding vehicle velocity uncertainty can be formulated as

$\begin{gathered}\eta_{i-1}(k)=\begin{bmatrix}n_{i-1}^{1}&n_{i-1}^{2}&\ldots&n_{i-1}^{m_{i-1}}\end{bmatrix}\begin{bmatrix}w_{i-1}^{1}(k)\\\
w_{i-1}^{2}(k)\\\ \vdots\\\ w_{i-1}^{m_{i-1}}(k)\end{bmatrix},\end{gathered}$

where $w_{i-1}^{j}(k)\,,j\in\\{1,2,\ldots,m_{i-1}\\}$ are auxiliary binary
variables that represent uncontrollable events and
$P[w_{i-1}^{j}(k)=1]=p_{i-1}^{j}$. Based on the above parametrization of the
uncertainty, the following equality should always hold true

$\sum_{j=1}^{m_{i-1}}w_{i-1}^{j}(k)=1.$

Using the auxiliary variables added to the system, the discrete-time state
space model in (6) is reformulated as

$\begin{gathered}\mathbf{x}_{i}(k+1)=\begin{bmatrix}1&t_{s}&-\tau\,t_{s}\\\
0&1&-t_{s}\\\
0&0&1-t_{s}\mathnormal{f}_{i}\end{bmatrix}\mathbf{x}_{i}(k)+\begin{bmatrix}0\\\
0\\\ t_{s}\mathnormal{f}_{i}\end{bmatrix}{u}_{i}(k)\\\
~{}~{}~{}~{}~{}+\begin{bmatrix}0\\\ 1\\\
0\end{bmatrix}\eta_{i-1}(k)+\begin{bmatrix}0\\\ t_{s}\\\
0\end{bmatrix}{a}_{i-1}(k).\end{gathered}$ (21)

If the acceleration profile is received through communication, constants
$n_{i-1}^{j}(k)$ will be simply set to zero (the problem turns into discrete
hybrid automata). Whenever GP model is used (between communication events or
during communication loss), $n_{i-1}^{j}(k)$ values are chosen based on the GP
output, and $a_{i-1}(k)$ can be easily calculated as
$\left[\hat{v}_{i-1}(k+1)-\hat{v}_{i-1}(k)\right]/t_{s}$ where
$\hat{v}_{i-1}(k)$ is the mean of the GP in (7).

##  IV Discrete Hybrid Stochastic MPC Problem for CACC Design

To achieve a string stable CACC system, $\mathbf{x}_{i}(k)$ should converge to
zero for every follower vehicle. This implies that the distance between
$i^{th}$ vehicle and its preceding vehicle converges to its desired value
while vehicles move with equal and constant speed. For the CACC design
problem, when the number of vehicles increases, centralized MPC is not
computationally efficient[16]. Besides, communication losses would make
centralized MPC unreliable and even dangerous for such safety-critical
systems. Instead, distributed MPC can be used to reach a string stable CACC.
The MPC problem for each vehicle is formulated as follows

$\min_{\textbf{u}_{i},\textbf{w}_{i-1},\textbf{z}_{i}}\\\
\sum_{k=0}^{N-1}\Big{[}(\mathbf{x}_{i}(k)-R_{i})^{T}\,Q_{i}\,(\mathbf{x}_{i}(k)-R_{i})\Big{]}-q_{i}\,\ln(\pi(\textbf{w}_{i-1}))\\\
\text{subject to:}~{}~{}~{}~{}~{}\text{MLD system equations,}\\\
\ln(\pi(\textbf{w}_{i-1}))\geq\ln(\tilde{p_{i}})~{}~{}~{}~{}~{}~{}~{}~{}$ (22)

where $\textbf{w}_{i-1}$ is a vector including all uncontrollable events,
$\textbf{u}_{i}$ and $\textbf{z}_{i}$ are the system inputs and the vector of
auxiliary variables from $k=0$ to $k=N-1$, respectively, and
$\pi(\mathbf{w}_{i-1})$ denotes the probability of a trajectory occurred by
$\mathbf{w}_{i-1}$ (see [13] for details).

In this paper, the CACC leverages the model-based communication in which
vehicles share their latest GP model with their follower vehicle. It is
assumed that each vehicle uses a low-rate communication to share its future
acceleration profile (model-free communication) with its follower vehicle.
Consequently, every vehicle has access to its preceding vehicle’s future speed
trajectory through either the communication (every $t_{c}$ seconds if packet
loss does not occur) or the GP model (until the next successful communication
event). Consequently, during communication loss or low-rate data exchange,
CACC will not degrade to ACC. If a vehicle uses the GP model, its optimization
problem is a discrete hybrid stochastic MPC while if the future speed profile
is received through communication, the optimization problem turns into a
discrete hybrid MPC. For simplicity, the controller that solely relies on the
model-free communication is called DHMPC, the controller that only relies on
the model-based communication is called DHSMPC, and the one that takes
advantage of both types of communication is named DH-DHSMPC.

##  V Simulations Results and Discussion

Simulations are conducted considering 10 vehicles and using multiple
scenarios. CVXPY package in Python is used for implementing the optimization
problem and Gurobi optimization package is used as the solver for the mixed
integer programs [17, 18, 19]. The desired speed trajectory for the leading
vehicle in the simulation is considered to be

$\begin{gathered}v^{*}_{0}(t)=\begin{cases}27&\,\,\,\,t<15\,s,\\\
0&\,\,\,\,15\,s\leq t<30\,s,\\\ 25&\,\,\,\,t\geq
30\,s.\end{cases}\end{gathered}$ (23)

Parameters used in these simulations can be found in Table I. Three different
case studies are investigated here. In the first two case studies, the
performance of the pure DHSMPC (using GP model, solving discrete hybrid
stochastic MPC) is compared with the DHMPC (sharing future acceleration
profile through communication, solving a discrete hybrid MPC). Finally, the
performance of the DH-DHSMPC, which leverages both model-based and model-free
communication, is compared with DHMPC when communication rate is low (1 Hz)
and an iid packet loss is considered.

TABLE I: Model and optimization parameters used in the simulations. parameter |  value |  parameter |  value
---|---|---|---
$N$ |  $7$ |  $t_{s}$ |  $0.1\,s$
$l_{i}$ |  $5\,m$ |  $d^{s}_{i}$ |  $2\,m$
$\underline{d}$ |  $1\,m$ |  $f_{i}$ |  $10\,s^{-1}$
$a_{i}^{min}$ |  $-4\,m/s^{2}$ |  $a_{i}^{max}$ |  $3\,m/s^{2}$
$u_{i}^{min}$ |  $-4\,m/s^{2}$ |  $u_{i}^{max}$ |  $3\,m/s^{2}$
$\hat{p}$ |  $0.01^{N}$ |  $q_{i}$ |  $10$
Figure 2: Comparing the performance of the DHSMPC (GP-based communication)
with DHMPC (model-free communication) when the time gap of $1\,s$ is
considered. From time $0\,s$ to $60\,s$, the platoon operates based on the GP
model while the environment resets at time $60\,s$ and after that vehicles
share their speed trajectory. It is observed that the performance of both
methods is almost the same, which demonstrates the efficacy of the GP-based
method. Figure 3: Comparing the performance of the DHSMPC (GP-based
communication) with DHMPC (model-free communication) when the time gap is
$0.7\,s$. From time $0\,s$ to $60\,s$, the platoon operates using the GP model
while the environment resets at time $60\,s$ and after that vehicles share
their speed trajectory. Although there are overshoots during acceleration in
the GP-based method and emergency braking occurs more frequently during
deceleration, the vehicles are still able to operate safely.

Figure 4: Performance of the DH-DHSMPC with 10 vehicles, $T=0.7\,s$,
$t_{c}=1\,s$, and communication loss with the probability of 0.25. Using the
GP model results in compensating both the low-rate communication and
communication failure.

Figure 5: Performance of the DHMPC with 10 vehicles, $T=0.7\,s$, $t_{c}=1\,s$,
and communication loss with the probability of 0.25. As observed, the vehicles
are not able to cope with the sudden deceleration, and an accident occurs.

For the DHMPC case, each vehicle employs the time-shifted acceleration
information of its preceding vehicle to efficiently use the last received data
between two successful communication events. Assuming that the last
communication occurred at time instant $k_{0}$, at $k_{1}<k_{0}+N-1$,
${a}_{i-1}(k)$ in (21) is replaced by

$\hat{a}_{i-1}(k)=\begin{cases}a_{i-1}(k)&\,k_{1}\leq k\leq k_{0}+N-1\\\
b(k)&\,\text{otherwise}\\\ \end{cases}$ (24)

where $b(k)$ can be estimated using a linear extrapolation of the last two
samples of $a_{i-1}(k)$. If there is a communication loss and $i^{th}$ vehicle
does not receive information from its preceding vehicle, then the vehicle
would resort to operating in adaptive cruise control (ACC) mode until it
receives information from its preceding vehicle and returns back to CACC mode.

In Figs. 2-5, the first subplot shows every vehicle’s distance from its
predecessor ($d_{i}(t)$) while the second subplot shows each vehicle’s
velocity ($v_{i}(t)$). The third subplot depicts each vehicle’s acceleration
information, and the last one shows the emergency braking status
($\xi^{E}_{i}(t)$). As shown in Fig. 2, in the first case study, when vehicles
try to keep a time gap of $1\,s,$ the performance of the both DHSMPC and DHMPC
is almost the same, and the vehicles do not need to perform emergency braking,
which shows the efficacy of the model-based communication. To evaluate the
performance of the model-based communication (DHSMPC) with smaller vehicle
headway, the time gap is reduced to $0.7\,s$ in the second case study. As
shown in Fig. 3, the vehicles with the proposed communication paradigm are
able to safely follow their preceding vehicle. However, compared to the model-
free communication case, vehicles need to perform emergency braking frequently
during deceleration, and there are noticeable overshoots during acceleration
in the distance/velocity profiles.

In the last simulation study, the performance of DHMPC is compared with DH-
DHSMPC. In this scenario, the time gap is $0.7\,s$, and vehicles send their
predictive acceleration trajectory every $1\,s$ ($t_{c}=1$) while the
communication success probability is $0.75$. Fig. 5 shows the performance of
the DH-DHSMPC method. Since this method takes advantage of the model-based
communication, vehicles are able to operate safely and efficiently with a low
communication rate and even in the presence of a communication failure. During
the deceleration event, a few vehicles enter the emergency braking mode for a
few seconds to assure safety. By comparing the results in Figs. 3 and 5, it is
perceived that the performance of the DH-DHSMPC with low-rate communication is
very similar to the performance of the DH-MPC with frequent and perfect
communication (second case study). On the other hand, as shown in Fig. 5, the
DHMPC is not able to preserve the safety of the platoon and during the sudden
deceleration, the second follower vehicle is not able to react properly and
hence crashes into its front vehicle.

##  VI Conclusion

In this paper, a discrete hybrid stochastic MPC design method was proposed for
CACC applications by leveraging model-based communication. The proposed method
aims at reducing the reliance of the vehicles on communication, i.e.,
operating safely during communication failure. It was assumed that vehicles
share their future acceleration profile as well as an updated model for their
speed profile (using Gaussian process) at each successful communication event.
For safety purposes, vehicles may operate in either free following mode or
emergency braking mode. The vehicle operating mode was chosen based on the
predictive speed profile of the preceding vehicle (either received through
communication or generated using the last GP model available). The performance
of the proposed controller was evaluated through simulation studies, which
validated the efficacy of the proposed method even with a low-rate
intermittent communication.

##  VII Acknowledgement

This research was supported by the National Science Foundation under grants
CNS-1931981 and CNS-1932037.

## References

  * [1] Joel Vander Werf, Steven E. Shladover, Mark A. Miller, and Natalia Kourjanskaia.  Effects of adaptive cruise control systems on highway traffic flow capacity.  Transportation Research Record, 1800(1):78–84, 2002.
  * [2] Steven E. Shladover, Dongyan Su, and Xiao-Yun Lu.  Impacts of cooperative adaptive cruise control on freeway traffic flow.  Transportation Research Record, 2324(1):63–70, 2012.
  * [3] Mahdi Razzaghpour, Shahriar Shahram, Rodolfo Valiente, and Yaser P. Fallah.  Impact of communication loss on mpc based cooperative adaptive cruise control and platooning, 2021.
  * [4] Yaser P. Fallah.  A model-based communication approach for distributed and connected vehicle safety systems.  In 2016 Annual IEEE Systems Conference (SysCon), pages 1–6, 2016\.
  * [5] Hossein Nourkhiz Mahjoub, Behrad Toghi, S M Osman Gani, and Yaser P. Fallah.  V2x system architecture utilizing hybrid gaussian process-based model structures.  In 2019 IEEE International Systems Conference (SysCon), pages 1–7, 2019.
  * [6] Ehsan Moradi-Pari, Hossein Nourkhiz Mahjoub, Hadi Kazemi, Yaser P. Fallah, and Amin Tahmasbi-Sarvestani.  Utilizing model-based communication and control for cooperative automated vehicle applications.  IEEE Transactions on Intelligent Vehicles, 2(1):38–51, 2017.
  * [7] H. N. Mahjoub, A. Raftari, R. Valiente, Y. P. Fallah, and S. K. Mahmud.  Representing realistic human driver behaviors using a finite size gaussian process kernel bank.  In 2019 IEEE Vehicular Networking Conference (VNC), pages 1–8, 2019\.
  * [8] Jeroen Ploeg.  Analysis and design of controllers for cooperative and automated driving.  2014\.
  * [9] Jeroen Ploeg, Elham Semsar-Kazerooni, Guido Lijster, Nathan van de Wouw, and Henk Nijmeijer.  Graceful degradation of cooperative adaptive cruise control.  IEEE Transactions on Intelligent Transportation Systems, 16(1):488–497, 2014.
  * [10] Chaoxian Wu, Yuan Lin, and Azim Eskandarian.  Cooperative adaptive cruise control with adaptive kalman filter subject to temporary communication loss.  IEEE Access, 7:93558–93568, 2019.
  * [11] Kelvin K Lee and Samuel T Chanson.  Packet loss probability for real-time wireless communications.  IEEE Transactions on Vehicular Technology, 51(6):1569–1575, 2002\.
  * [12] Youssef Abou Harfouch, Shuai Yuan, and Simone Baldi.  An adaptive switched control approach to heterogeneous platooning with intervehicle communication losses.  IEEE Transactions on Control of Network Systems, 5(3):1434–1444, 2017.
  * [13] Alberto Bemporad and Stefano Di Cairano.  Model-predictive control of discrete hybrid stochastic automata.  IEEE Transactions on Automatic Control, 56(6):1307–1321, 2010.
  * [14] Gerrit JL Naus, Rene PA Vugts, Jeroen Ploeg, Marinus JG van De Molengraft, and Maarten Steinbuch.  String-stable cacc design and experimental validation: A frequency-domain approach.  IEEE Transactions on vehicular technology, 59(9):4268–4279, 2010\.
  * [15] Alberto Bemporad and Manfred Morari.  Control of systems integrating logic, dynamics, and constraints.  Automatica, 35(3):407–427, 1999.
  * [16] Yang Zhou, Meng Wang, and Soyoung Ahn.  Distributed model predictive control approach for cooperative car-following with guaranteed local and string stability.  Transportation research part B: methodological, 128:69–86, 2019\.
  * [17] Steven Diamond and Stephen Boyd.  CVXPY: A Python-embedded modeling language for convex optimization.  Journal of Machine Learning Research, 17(83):1–5, 2016.
  * [18] Akshay Agrawal, Robin Verschueren, Steven Diamond, and Stephen Boyd.  A rewriting system for convex optimization problems.  Journal of Control and Decision, 5(1):42–60, 2018.
  * [19] LLC Gurobi Optimization.  Gurobi optimizer reference manual, 2021.


# Superfluid rings as quantum pendulums

Antonio Muñoz Mateo  Departamento de Física, Universidad de La Laguna, 38200
Tenerife, Spain  Grigory E. Astrakharchik  Departament de Física, Universitat
Politècnica de Catalunya, 08034 Barcelona, Spain  Departament de Física
Quàntica i Astrofísica, Universitat de Barcelona, Martí Franquès 1, 08028
Barcelona, Spain.  Institut de Ciències del Cosmos (ICCUB), Universitat de
Barcelona, Martí Franquès, 1, 08028 Barcelona, Spain.  Bruno Juliá-Díaz
Departament de Física Quàntica i Astrofísica, Universitat de Barcelona, Martí
Franquès 1, 08028 Barcelona, Spain.  Institut de Ciències del Cosmos (ICCUB),
Universitat de Barcelona, Martí Franquès, 1, 08028 Barcelona, Spain.

###### Abstract

A feasible experimental proposal to realize a non-dispersive quantum pendulum
is presented. The proposed setup consists of an ultracold atomic cloud,
featuring attractive interatomic interactions, loaded into a tilted ring
potential. The classical and quantum domains are switched on by tuned
interactions, and the classical dynamical stabilization of unstable states,
i.e. a la Kapitza, is shown to be driven by quantum phase imprinting. The
potential use of this system as a gravimeter is discussed.

Ring geometries are omnipresent in physics. Mathematically, they endow systems
with periodic boundary conditions; physically, they realize the minimal block
of cyclic transport, which would become perpetual if there were no
dissipation. Approaching the dissipationless limit, superconductors and
superfluids are capable of making the cyclic transport of charge or particles,
if not perpetual at least persistent, a particularly striking demonstration of
which is the persistent flow of superconducting gravimeters Van Camp  _et al._
(2017). In this regard, the first realizations of ring geometries in ultracold
gases opened new avenues for experiments with persistent currents of highly
controlled Bose-Einstein condensates Ryu  _et al._ (2007); Heathcote  _et al._
(2008); Ramanathan  _et al._ (2011); Moulder  _et al._ (2012). However, in
this case, the gravitational pull is an apparent hindrance to stationary
flows, since a small tilting of the ring axis with respect to the vertical
direction gives rise to an unwanted azimuthal potential for the trapped atoms
Heathcote  _et al._ (2008); Pandey  _et al._ (2019). Interestingly, although
suppression of the tilting is necessary for the study of unobstructed
currents, by letting the tilting occur the ring is transformed into a pendulum
(see Fig. 1); for instance, horizontally setting the axis of a typical ring of
radius $R=20\,\mu$m converts it into a pendulum of angular frequency
$\omega=\sqrt{g/R}=700$ rad/s, where $g$ is the gravity acceleration. While
for a repulsively interacting, extended Bose-Einstein Condensate (BEC) the
resemblance to the classical pendulum would only apply to the motion of the
center of mass, one expects to find its quantum analogue in attractively
interacting, localized BECs. The present work is devoted to exploring the
validity of this analogy.

Figure 1: (a) Schematic picture, side (left) and front (right) views, of a
tilted ring potential (purple isocontour) inside which a BEC (cream blob)
moves in the presence of the constant gravitational force $F_{\mbox{\tiny
G}}$. (b) Ground-state density profile in a ring trap, characterized by the
characteristic frequency $\omega=2\hbar/mR^{2}$ and potential depth $s=32$
(see text), for various strengths of the interparticle attraction,
parameterized by $\tilde{\gamma}=|\gamma|mNR/\hbar^{2}$. (c) Oscillations of
the center of mass (top panel) and the average momentum (bottom panel) of a
bright soliton of width $\xi=R/4$ (or, equivalently, $\tilde{\gamma}=8$) for
different initial positions. The shaded areas around the solid lines represent
the mean width in position $\sigma$ and momentum $\sigma_{k}$ spaces. (d)
Trajectories in the phase space. States at $t=0$ (triangles) and $t=T$
(circles) are indicated by open symbols.

The pendulum dynamics has been addressed under multiple perspectives in the
context of ultracold atomic systems. Its presence is implicitly evident in the
coherent tunneling of particles, as observed in the Josephson effect Smerzi
_et al._ (1997); Fialko  _et al._ (2015); Pigneur and Schmiedmayer (2018).
Additionally, the proposal for the dynamic generation of nonlinear excitations
has shed further light on the intricacies of pendulum dynamics in this context
Fialko  _et al._ (2012). Most of the studies have focused on the dynamical
stabilization of pendulum-like equilibrium states in optical potentials that
are unstable Gilary  _et al._ (2003); Jiang  _et al._ (2023); He and Liu
(2023). In a classical context, the pendulum dynamics of a microscopic
colloidal particle in a ring built with optical tweezers has been recently
reported Richards  _et al._ (2018).

The tilted ring. We assume a quasi-one-dimensional character of the system,
which can be ensured by imposing tight transverse confinement to the atoms.
While the role of quantum fluctuations is enhanced in one dimension, still, in
the mean-field regime, the gas is sufficiently coherent to be correctly
described by the Gross-Pitaevskii equation. Although strictly speaking Bose-
Einstein condensation is absent in one dimension, properties of a
quasicondensate are correctly captured by equations describing the evolution
of the condensate wave function. A tilting angle $\tau\in[0,\,\pi/2]$ produces
on the particles of mass $m$ the gravitational potential
$U(\tau,\theta)=-m\,g\,R\sin\tau\,\cos\theta$ along the azimuthal coordinate
$\theta\in[-\pi,\,\pi]$ of the ring, or alternatively
$U(\tau,\theta)=-m\omega^{2}R^{2}\cos\theta$, where
$\omega=\sqrt{g\sin\tau/R}$ is the pendulum angular frequency, so
$T=2\pi/\omega$ is the corresponding time period, see the depiction in Fig. 1.
For non-interacting particles, the time evolution is described by the
Schrödinger equation, i.e. Eq. 1 with $\gamma=0$ (see the early discussion of
pendulums in quantum mechanics by Condon Condon (1928), or for a more recent
account, for instance, Ref. Doncheski and Robinett (2003)). It admits general
solutions that can be written as superpositions of Mathieu functions DLMF
${\rm ce}_{2n}(\theta,q)$ and ${\rm se}_{2n+1}(\theta,q)$, with
$n=0,1,2,\dots$, and $2q=-(2R/a_{\rm ho})^{4}$. Here we have introduced the
harmonic oscillator length, $a_{\rm ho}=\sqrt{\hbar/m\omega}$, associated with
small-amplitude oscillations around the equilibrium point $\theta=0$. For
$q\rightarrow 0$, that is for a small tilting or $\omega\rightarrow 0$, these
functions tend to the trigonometric functions $\cos(2n\theta)$ and
$\sin[(2n+1)\theta]$.

Deeper insights can be drawn by considering the particle positions along the
circumference of the rings. The potential felt by a particle located at
position $x$ can be expressed by substituting the azimuthal angle by
$\theta=2\pi x/d$, where $d=2\pi R$ is the circumference of the ring, leading
to $U(\tau,x)=-gmR\sin\tau+2gmR\sin\tau\sin^{2}(\pi x/d)$. The resulting
tilting potential can be conceptualized Muñoz Mateo  _et al._ (2019) as a
periodic lattice potential $V_{0}\sin^{2}(\pi x/d)$ of amplitude
$V_{0}=2gmR\sin\tau=2m\omega^{2}R^{2}$. In this context, the natural energy
scale in a lattice is given by
$E_{L}=\hbar^{2}\pi^{2}/(2md^{2})=\hbar^{2}/(8mR^{2})$, referred to as the
recoil energy Bloch  _et al._ (2008). The dimensionless ratio between lattice
amplitude and the recoil energy, denoted as
$s=(1/2)\;V_{0}/E_{L}=8m^{2}\omega^{2}R^{4}/\hbar^{2}$ (with an extra $1/2$
factor as compared to the usual convention Morsch and Oberthaler (2006) ),
quantifies the lattice depth and matches in absolute value the Mathieu-
function parameter $s=|q|$. A shallow lattice is characterized by $s\lesssim
1$ while a deep lattice corresponds to $s\gg 1$. As in the tilted ring only a
single lattice site is available, i.e. $0\leq x\leq d$, the ring displays a
minimal band structure composed of a single Bloch state, with zero
quasimomentum, per energy band Kittel and McEuen (2018). In Figure
LABEL:fig:linear(a) two opposite limits are shown, corresponding to shallow
($s=0.5$) and deep ($s=32$) lattices. The three lowest-energy eigenstates
$\psi_{n,0}$ are shown with $n=0,1,2,\dots,$ standing for the band index and
reflecting the number of nodes. As expected, the density becomes increasingly
localized as the lattice amplitude is augmented. In the limit of very large
$s$, the lowest eigenstates are similar to those of a harmonic oscillator.

For the interacting gas, the system dynamics obeys the time-dependent Gross-
Pitaevskii equation

$i\hbar{\partial_{t}\psi}=\left(-\frac{\hbar^{2}}{2mR^{2}}\partial_{\theta}^{2}-m\omega^{2}R^{2}\cos\theta+\gamma|\psi|^{2}\right)\,\psi,$
(1)

where $\gamma=-2\hbar^{2}/(ma_{1D})<0$ denotes the coupling constant of the
attractive interparticle interaction and $a_{1D}$ is the $s$-wave scattering
length Olshanii (1998). The number of particles fixes the wave function
normalization $N=R\int d\theta\,|\psi|^{2}$, while the energy $E$ is provided
by the functional ${E[\psi]}=R\int
d\theta\left({\hbar^{2}|\partial_{\theta}\psi|^{2}}/{2mR^{2}}-m\omega^{2}R^{2}\cos\theta\,|\psi|^{2}+{\gamma}|\psi|^{4}/2\right)$.
Expectation values are computed as quantum-mechanical averages, $\langle
A\rangle=(1/N)R\int d\theta\,\psi^{*}A\psi$, where the $1/N$ appears due to
the wave function normalization. Figure 1(b) illustrates the effect of
increasing attractive interactions on the ground-state density profile, which
becomes increasingly localized. When the interparticle interaction dominates
over the external potential, the ground state approaches a free bright soliton
centered at $\theta_{0}=0$,

$\displaystyle\psi_{BS}(\theta)=\sqrt{\frac{N}{2\xi}}\,\mbox{sech}\left(R\frac{\theta-\theta_{0}}{\xi}\right),$
(2)

whose width $\xi=2\hbar^{2}/(m|\gamma|N)=a_{1D}/N$ scales inversely with the
number of particles, while the peak density does quadratically.

The close analogy between this system and the classical pendulum is evidenced
by the time evolution of expectation values, as obtained from the Ehrenfest
equations (see Appendix A). An illustrative example is represented in Fig.
1(c), where the evolution of the center of mass, $\langle\theta\rangle$, and
the angular momentum, $\langle\hat{p}_{\theta}\rangle$, is shown for a soliton
of width $\xi=R/4$ in a deep lattice with $s=32$. The time evolution follows
the equations of motion
$mR^{2}\,d\langle\theta\rangle/dt=\langle\hat{p}_{\theta}\rangle$ and
$d\langle\hat{p}_{\theta}\rangle/dt=-\langle\partial_{\theta}U\rangle$ and, in
agreement with Kohn’s theorem, is independent of the interaction strength.
Hence, the nonlinear equation of the pendulum is
$d^{2}\langle\theta\rangle/dt^{2}+\omega^{2}\langle\sin\theta\rangle=0$.
Nevertheless, the analogy with the classical system is not always fulfilled.
Had we chosen a wider wave packet (say $\xi=2R/5$) centered at the potential
maximum as the initial state, it would spread over the ring, since the wave
packet falls onto both sides of the maximum, as illustrated in Fig.
LABEL:fig:linear(b). We show next that this observation reflects the lack of
the classical equilibrium point at $\theta_{0}=\pi$, whose existence depends
on the system parameters.

Nonlinear spectrum. The lattice viewpoint offers valuable insight into the
structure of the spectrum of the nonlinear system. By switching on the
interaction, the linear Bloch states Kittel and McEuen (2018) characterized by
a given quasimomentum, extend into the nonlinear regime while maintaining the
number of nodes, which stands as a distinctive topological feature. In
addition, nonlinear lattice systems allow for new stationary states of
different nature, namely the gap solitons. Contrary to Bloch waves in an
infinite lattice, which are extended states over the whole system, gap
solitons are localized states that occupy a few lattice sites (see Ref. Morsch
and Oberthaler (2006), and references herein). In static systems, such states
emerge within the energy gaps of the underlying linear system beyond an
interaction threshold Louis  _et al._ (2003); Muñoz Mateo  _et al._ (2019). In
the tilted ring, gap solitons emerge as states centered at the potential
maximum.

Figure 3 illustrates the characteristic differences in the spectrum between
shallow ($s=0.5$) and deep ($s=32$) lattices. The chemical potential is
reported as a function of the interaction strength. The states labeled by
$\psi_{n,0}$ represent the nonlinear continuation of the linear solutions
shown in Fig. LABEL:fig:linear. These states have the center of mass located
at the potential minimum $\theta=0$. Beyond a certain interaction threshold, a
new nonlinear state emerges (shown by the continuous line with symbols and
labeled as $\psi_{0,\pi}$) with location at the potential maximum. Thus, these
states are centered in the classically-unstable equilibrium position of the
inverted pendulum with $\theta=\pi$.

Figure 3:  Chemical potential of stationary states
$\tilde{\mu}=\mu/(\hbar^{2}/mR^{2})$ versus interaction strength
$\tilde{\gamma}=|\gamma|mNR/\hbar^{2}$, for (a) deep potential $s=32$, and (b)
shallow potential $s=0.5$. The linear energy bands of an infinite lattice with
the same spacing (yellow shaded regions) are shown as a reference. The top
horizontal axis of each graph indicates the width of a free soliton $\xi$ in
units of the harmonic oscillator length $a_{\rm ho}$.

The emergence of equilibrium states at the potential maximum can be understood
by expanding the potential around $\theta=\pi$, up to second order in the
position. This results in a expulsive harmonic potential
$U(\pi+\delta\theta)-U(\pi)\approx-m\omega^{2}R^{2}\delta\theta^{2}/2$. Here,
the soliton Eq. (2) centered in $\theta_{0}=\pi$, with the soliton width $\xi$
as variational parameter, provides an ansatz. The energy in this system is
given by the standard Gross-Pitaevskii energy functional (see, for example,
Ref. Pitaevskii and Stringari (2003) for details). By performing the energy
minimization we arrive to the following quartic equation
$1-(a_{ho}/\xi_{f})\xi/a_{ho}+(\pi/4)(\xi/a_{ho})^{4}=0$, where
$\xi_{f}=2\hbar^{2}/(m|\gamma|N)$ is the free soliton width. This expression
is valid when both $\xi\ll R$ and $\xi_{f}\ll R$. For a wide soliton with
$\xi\gg a_{ho}$, no real solution exists, indicating an absence of minimum in
the energy functional and a lack of a stationary state in the vicinity of the
potential maximum. On the contrary, for a narrow soliton with $\xi\ll a_{ho}$,
the stationary solution exists with its width being roughly the same as that
of a free soliton $\xi\approx\xi_{f}$. The stationary solution exists when the
soliton width is smaller than the threshold value
$\xi_{f}/a_{ho}=(0.75/\pi^{1/3})^{0.75}\approx 0.6$, or equivalently, for
interaction strengths larger than $|\gamma|N\approx
3.3\sqrt{\hbar^{3}\omega/m}$. These values provide a reasonable estimate of
the threshold value, $\xi_{f}/a_{ho}\approx 0.45$, shown in Fig. 3.

Dynamical stabilization. The classical pendulum possesses only a single stable
position, corresponding to the pendulum stabilized at the bottom,
$\theta_{1}=0$, while the equilibrium of the pendulum in the top position,
$\theta_{2}=\pi$, is dynamically unstable. Figure 1(c) illustrates that a
similar phenomenon occurs in the tilted ring: a weak perturbation, generated
by adding a sinusoidal wave on the initial state that shifts the center of
mass to $\theta_{0}=-3.139$, makes the soliton roll down. Fascinatingly, the
inverted position of the pendulum can be stabilized by inducing fast vertical
vibrations (along the gravitational field direction) of the pendulum pivot as
experimentally shown and mathematically proved by introducing the method of
averaging of the fast variables by Pyotr Kapitza in his influential paper
Kapitza (1951). Commonly, the driven pendulum is referred to as a Kapitza
pendulumAstrakharchik and Astrakharchik (2011). In contrast, the horizontal
vibration moves the equilibrium point to a new position $\theta_{2}\neq\pi$
(see Ref. Landau and Lifshitz (1960), $\S 5,\,\S 30$). The potential energy of
the pivot vibration of frequency $\Omega\gg\omega$ along a generic angle
$\alpha$ can be calculated as the work $W=-{\bf f}\cdot{\bf r}$ done by the 2D
oscillating force ${\bf
f}=(\cos\alpha,\,\sin\alpha)\,m\Omega^{2}\ell\,\sin(\Omega\,t)$ acting within
the ring plane, where $\ell$ $\ll R$ is the characteristic amplitude of the
vibration, and ${\bf r}=R(\cos\theta,\,\sin\theta)$ are the ring coordinates.
As a result, the total potential energy becomes

$\displaystyle
U(\theta,t)=-m\omega^{2}R^{2}\left\\{\left[1-\frac{\Omega^{2}\ell}{\omega^{2}R}\sin(\Omega
t)\cos\alpha\right]\cos\theta\right.$
$\displaystyle\left.-\frac{\Omega^{2}\ell}{\omega^{2}R}\sin(\Omega
t)\sin\alpha\,\sin\theta\right\\}.$ (3)

The idea of Kapitza is based on separating processes occurring at slow and
fast time scales, i.e. corresponding to low ($\omega$) and high ($\Omega$)
frequencies. Then, the slow motion is governed by the effective time-
independent potential obtained by time averaging over the fast oscillations,

$\displaystyle U_{\rm
eff}(\theta)=-m\omega^{2}R^{2}\left\\{\left[\cos\theta-\left(\frac{\beta\cos\alpha\sin\theta}{2}\right)^{2}\right]\right.$
$\displaystyle\left.-\left(\frac{\beta\sin\alpha\cos\theta}{2}\right)^{2}\right\\}.$
(4)

where $\beta={\Omega\ell}/({\omega R})$ is the ratio between the velocity that
characterizes the vibration, $\Omega\ell$, and a velocity associated with the
unperturbed pendulum, $\omega R$. The transition to (classical) stable states
takes place at $\beta=\sqrt{2}$. The stability criterion can be stated in
terms of energetic quantities, i.e. that the kinetic energy, induced by the
driven oscillations, should be large compared to the potential energy above
the pivotal point. Figure LABEL:fig:Ueff(a) presents characteristic examples
illustrating how the effective potential depends on the driving strength
$\beta$ in a deep lattice (specifically, we use $s=32$). In Fig.
LABEL:fig:Ueff(a) we focus on the case of vertical driving, i.e. $\alpha=0$.
While for weak driving vibrations ($\beta=0$ and $\beta=1$) only a single
minimum exists, corresponding to the usual bottom position of the pendulum,
$\theta=0$, for strong driving ($\beta=2$) the pendulum might flip and a
second stable minimum is formed for $\theta=\pm\pi$, corresponding to an
inverted pendulum. The case of horizontal driving, $\alpha=\pi/2$, is
illustrated in Fig. LABEL:fig:Ueff(b). For weak driving oscillations, the
usual single minimum is observed. Instead, for sufficiently strong driving,
the new (classical) minima appear at $\theta_{1,2}=\pm\pi/3$, while $\theta=0$
becomes a local maximum.

We have tested the stability of quantum states [a wave packet similar to the
one described in Fig. 1(c)] subject to the time-dependent potential Eq. (3)
with vibration amplitude $\ell=0.1\,R$, and centered at the classical
equilibrium points of Eq. (4); the outcome is presented in Fig.
LABEL:fig:Ueff(b). In all considered cases, oscillations of frequency $\Omega$
can be observed as fast beating in the evolution of average momentum. For
$\alpha=0$ and $\beta=2$ (that is $\Omega=20\,\omega$) the equilibrium point
$\theta_{2}=\pi$ is stable (solid blue curve), whereas for $\beta=1$ (not
shown) it is not, in agreement with the classical prediction. The situation is
not that clear for $\alpha=\pi/2$, since $\beta=2$ (and $\phi=0$, see below)
does not lead to a static state; instead, the state is induced to tunnel
through the local maximum separating the two minima (dashed line), and higher
vibrations (for instance $\beta=3$, not shown) produce only a partial self
trapping. However, contrary to the classical case, the phase of the vibration
is key for the quantum pendulum; introducing a phase factor $\phi$ in Eq. (3)
by the substitution $\Omega t\rightarrow\Omega t+\phi$ has a crucial influence
on stability: although $\phi=0$ does not reproduce the classical result,
$\phi=\pi/2$ does (dash-dotted line). The cause resides in the velocity of the
initial state (or equivalently its phase) induced by the fast oscillations and
controlled by this phase factor, by virtue of which $\phi=\pi/2$ produces a
zero velocity. This mechanism is not qualitatively distinct from the usual
phase imprinting technique employed in ultracold-gas experiments Dobrek  _et
al._ (1999), and stands as an additional, key stabilization effect in quantum
systems, along with the effective classical potential Eq. (4), of Kapitza’s
procedure.

Measuring gravity. Historically, pendulums were the first devices used to
measure gravity, both for absolute and relative measurements, and also the
most accurate ones until the second half of the twentieth century, reaching
values of $\delta g\approx 10^{-6}g$; afterwards, they were replaced by free-
fall apparatus Zumberge (2021). Modern versions of the latter are made of
ultracold atoms Debs  _et al._ (2011), and have reached top-level performance
$\delta g\approx 10^{-9}g$ Ménoret  _et al._ (2018); Stray  _et al._ (2022).
By means of atomic interferometry, based on a sequence of Raman pulses that
split and reunite the falling atomic clouds Kasevich and Chu (1991), high-
contrast fringes are produced that provide the acceleration of gravity.

Differently, pendulum gravimetry relies on measuring the oscillation period,
hence length, positions and corresponding times have to be tracked over
repeated, small oscillations. In the present quantum pendulum, where
dissipation can be ruled out, the atomic cloud can be strongly localized; for
instance, an atomic cloud of $N=2.5\times 10^{4}$ particles of 7Li with
interaction strength $\gamma=-2\pi\nu\hbar a_{0}$, where $a_{0}$ is the Bohr
radius, and harmonic transverse confinement of frequency $\nu=350$ Hz [similar
parameters as in the experiment of Ref. Nguyen  _et al._ (2017), and away from
the critical particle number for collapse conditions $N_{c}\approx 4.5\times
10^{4}$ see Refs. Pérez-García  _et al._ (1997, 1998)], has a typical size of
$\xi\approx 0.6\,\mu$m (doubling $\nu$ reduces $\xi$ by half and $N_{c}$ by
$\sqrt{2}$). Even within a small ring of radius $R=12\,\mu$m, as those of Ref.
Beattie  _et al._ (2013), the relevant parameter $R/\xi=20$ is large enough to
closely reproduce the dynamics of a classical pendulum (with expected relative
differences of $\delta\sim 5\times 10^{-4}$ in the motional periods, see
Appendix A). In this case, assuming that shining far-from-resonance light
perpendicular to the ring on the atomic cloud does not modify the dynamics, as
for classical particles, laser photo-detection could be used to track the
motion with expected accuracy of, at least, the order or a few percent (as
typically reported in oscillation measurements Fang  _et al._ (2014); Huang
_et al._ (2019)); in-trap non-destructive absorption imaging Nguyen  _et al._
(2017) could also be performed. This uncertainty in measuring the motional
period might be reduced by resorting to similar procedures as for the
classical reversible pendulums of Kater and Bessel Poynting and Thomson
(1908): the measurement of two close but different periods. In the ring, they
are associated with different tilts, or equivalently to different lengths
$L_{j}=R/\sin\tau_{j}$, for $j=1,\,2$, so that gravity is obtained from
$g=(2\pi)^{2}\,(L_{1}^{2}-L_{2}^{2})/(L_{1}T_{1}^{2}-L_{2}T_{2}^{2})$.

A precise determination of the period would also allow one to use the tilted
ring for performing sensitive measurements of gravity Stray  _et al._ (2022).
By initially preparing the BEC at the energy maximum $\theta=\pi$, where small
displacements produce large differences in the period, the presence of a
sensible mass at the surface, horizontally located with respect to the ring
axis such that $\delta g$ and $g$ could act orthogonally, would induce an
angular displacement $\delta\theta\approx\tan\delta\theta=\delta g/g$. This
translates into pendulum periods [see Appendix A, Eq. (7)] that can be
approximated by $T_{k}(\delta\theta)\approx(2/\pi)(3\ln
2-\ln\delta\theta)\,T$; for instance, for a void of mass of 100 kg at 1-m
apart one obtains $\delta g/g\approx 7\times 10^{-10}$ and $T_{k}=145.5\,T$,
while for twice that mass, $\delta g/g\approx 14\times 10^{-10}$, the period
is $T_{k}=141.2\,T$.

## Appendix A

###  A.1  Numerical solutions

Equation (1) has been numerically solved through FFT techniques for the
spatial discretization, and standard, high-order (typically 5 to 9) time
integrators of Julia programming language. For time-independent solutions,
both imaginary time evolution, for the ground state, and Newton method, for
excited states as those presented in Fig. 3, have been used.

###  A.2  Adiabatic switch of fast vibrations

Instead of suddenly turning on the fast vibrations, as described in the main
text, the dynamical stabilization of equilibrium points can also be performed
in an adiabatic way. To this end, we have chosen a setup with similar
parameters as the ones used in Fig. LABEL:fig:Ueff(b) and $\alpha=\pi/2$, but
with an initial soliton state of varying amplitude (according to
$\tilde{\gamma}=8$ or 16), centered at the equilibrium point $\theta=0$ of the
non-vibrated system. Subsequently, the amplitude and high frequency of the
vibrating potential are adiabatically ramped up,
$U(\Omega)=f(t)\,m\Omega^{2}\,\ell R\,\sin[f(t)\,\Omega t+\phi]\,\sin\theta$,
where $f(t)=\tanh(0.1\,\omega\,t)$. Figure 5 shows some characteristic
examples for the time evolution of the center of mass in the adiabatic case.
As it can be inferred from this figure, the phase factor $\phi$ still plays a
relevant role, and eventually, for $\phi=\pi/2$, the initial state rolls
downhill in the effective potential generated by the fast vibration towards
one of its energy minima at $\theta=\pi/3$; higher and less regular
oscillations around the energy minimum are observed for the more localized
soliton (at $\tilde{\gamma}=16$) due to the higher momentum uncertainty.

Figure 5: Time evolution of the center of mass position in the case of an
adiabatic switch of the fast vibration along the direction $\alpha=\pi/2$. The
tilted ring is characterized by a potential depth of $s=32$. The initial
soliton states, featured by the interaction $\tilde{\gamma}$ and situated at
$\theta=0$, move towards the effective energy minima (indicated by horizontal,
thin dashed lines) depending on the vibration phase $\phi$ (see text for
additional details).

###  A.3  Classical vs quantum periods of oscillation

For oscillatory motion, the classical equation of the pendulum
$\ddot{\theta}+\omega^{2}\sin\theta=0$ has general solutions in terms of
elliptic functions $\theta(t)\sim\arcsin\left[k\,\mbox{sn}(\omega t,k)\right]$
(see, for instance Refs. Symon (1971); Tabor (1989)), where $\mbox{sn}(\omega
t,k)$ is the Jacobi elliptic sine function DLMF  of modulus
$k=\sin(\theta_{m}/2)$, with $\theta_{m}$ being the maximum, turning-point
angle. For generic initial conditions, $\theta(0)=\theta_{0}$ and
$\dot{\theta}(0)=\dot{\theta}_{0}$, the solution reads

$\displaystyle\theta(t)=\theta_{0}+$ $\displaystyle
2\left\\{\,\arcsin\left[k\,\mbox{sn}(\omega t+\varphi_{0},k)\right]\right.$
$\displaystyle\left.\qquad-\arcsin\left[k\,\mbox{sn}(\varphi_{0},k)\right]\,\right\\},$
(5)

from which, it follows the angular velocity

$\displaystyle\dot{\theta}(t)=2k\omega\,\mbox{cn}(\omega t+\varphi_{0},k),$
(6)

where cn is the Jacobi elliptic cosine function, hence
$\varphi_{0}=\mbox{arc\,cn}[\dot{\theta}_{0}/(2k\,\omega),\,k]$. The pendulum
period is given in terms of the complete elliptic integral of the first kind
$K(k)$ DLMF  as

$\displaystyle T_{k}=\frac{K(k)}{\pi/2}\,T,$ (7)

such that it differs from the period $T=2\pi/\omega$, achieved in the
approximation of small oscillations, by the factor
$K(k)/(\pi/2)\in[1,\,\infty]$ (a monotonically increasing function of $k$) for
$k\in[0,\,1]$, that is for $\theta_{m}\in[0,\,\pi]$.

Figure LABEL:fig:Tcomp shows a comparison between the classical predictions of
Eqs. (5-7) and the mean values of the motion of free-soliton wave packets 2 in
a tilted ring. As one could expect, the period predicted by Eq. (7),
$T_{k}=1.0732\,T$, is better approached at higher interactions
$\tilde{\gamma}=16$. By averaging over 25 periods of oscillation in order to
minimize the uncertainty in the measured times, we obtain the differences
$\delta(\tilde{\gamma})=[T_{k}(\tilde{\gamma})-T_{k}]/T_{k}$ to be
$\delta(\tilde{\gamma}=16)=3.36\times 10^{-3}$, and
$\delta(\tilde{\gamma}=8)=1.23\times 10^{-2}$; a slight reduction is found for
initial conditions closer to equilibrium positions, for instance,
$\delta(\tilde{\gamma}=16)=3.14\times 10^{-3}$ for $\theta_{0}=3\pi/4$.

These differences can be understood by considering the Ehrenfest equation
$d^{2}\langle\theta\rangle/dt^{2}+\omega^{2}\langle\sin\theta\rangle=0$, which
exactly matches the functional form of the classical equation, as
$d^{2}\langle\theta\rangle/dt^{2}+\omega^{2}\sin\langle\theta\rangle=0$, only
when the density profile of the soliton wave function becomes a Dirac delta
function
$|\psi_{BS}(\theta-\theta_{0})|^{2}\rightarrow\delta(\theta-\theta_{0})$;
otherwise, the difference between both equations can be written as the power
series

$\displaystyle\Delta=\,$
$\displaystyle\langle\sin\theta\rangle-\sin\langle\theta\rangle=$
$\displaystyle\sum^{n\geq 3}_{n\,\mbox{\tiny
odd}}\frac{(-1)^{\frac{n-1}{2}}}{n!}\,\left(\langle\theta^{n}\rangle-\left\langle\theta\right\rangle^{n}\right).$
(8)

For a free soliton state Eq. (2) of width $\xi$ and moving center
$\theta_{0}(t)$, the above series can be approximated, for large ratios
$(R/\xi=\tilde{\gamma}/2)\gg 1$, up to second order, by

$\displaystyle\Delta(\tilde{\gamma})\approx-\frac{\pi^{2}}{6\,\tilde{\gamma}^{2}}\,\sin\langle\theta\rangle,$
(9)

where we made use of ${\langle\theta\rangle}=\theta_{0}(t)$. This result
amounts to having an effective classical equation with a slightly reduced
angular frequency
$\omega(\tilde{\gamma})=\omega\,\sqrt{1-{\pi^{2}}/({6\,\tilde{\gamma}^{2}})}$,
or, equivalently, a slightly increased small-oscillation period

$\displaystyle
T(\tilde{\gamma})\approx\left(1+\frac{\pi^{2}}{12\,\tilde{\gamma}^{2}}\right)\,{T}\,.$
(10)

Evaluated at the values $\tilde{\gamma}=8,\,16$, used in Fig. LABEL:fig:Tcomp,
this estimate produces $\delta(\tilde{\gamma}=8)=1.31\times 10^{-2}$, and
$\delta(\tilde{\gamma}=16)=3.23\times 10^{-3}$, in good agreement with the
measured results previously reported.

Ring roundness.- Small azimuthal variations in the ring potential can also
give rise to alterations in the period of the motion. In order to estimate the
effect of such variation, we have introduced a perturbed potential
$U_{n}(\theta)=-m\omega^{2}R^{2}[\cos\theta-\epsilon_{n}\cos(n\,\theta)]$,
with a small integer $n>1$, and $\epsilon_{n}\ll 1$. It translates into an
extra force term
$F_{n}(\theta)=\epsilon_{n}m\omega^{2}R\,n\langle\sin(n\theta)\rangle$ in the
corresponding Ehrenfest equation that can also be calculated as a power series
up to second order in $\tilde{\gamma}^{-1}$ as

$\displaystyle
F_{n}(\tilde{\gamma})\approx\epsilon_{n}\,m\omega^{2}R\,n\,\left(1-\frac{n^{2}\pi^{2}}{6\,\tilde{\gamma}^{2}}\right)\,\sin\langle
n\,\theta\rangle\,.$ (11)

In the limit of small oscillations, a perturbed single frequency,
$\omega_{n}(\tilde{\gamma})=\omega\sqrt{1-\epsilon_{n}[n^{2}-{n^{4}\pi^{2}}/({6\,\tilde{\gamma}^{2}})]}$,
can be obtained; for instance, by setting $n=2,\,3$, and $\epsilon_{n}=0.01$,
this estimate provides us with the order of magnitude of the perturbed period
observed in our numerical results, as shown in Fig. LABEL:fig:Tcomp(c), where
we have measured $\delta_{2}(\tilde{\gamma}=16)=1.57\times 10^{-2}$ and
$\delta_{3}(\tilde{\gamma}=16)=1.16\times 10^{-2}$. From the reported data, it
becomes clear that the azimuthal variations of the ring potential can
introduce a significant source of uncertainty in the measured observables, at
least in the search for high precision measurements (see e.g. Ref. de Goër de
Herve  _et al._ (2021)); in this regard, time-averaging optical ring
potentials Bell  _et al._ (2016) could be useful to achieve improved trap
smoothness.

## Acknowledgements

The authors are grateful to Jean Dalibard for careful reading of a draft of
this paper and offering insightful comments. We are indebted to Joan Martorell
for multiple discussions, careful reading of drafts, and comprehensive help.
This work has been funded by Grants No. PID2020-114626GB-I00 and
PID2020-113565GB-C21 by MCIN/AEI/10.13039/5011 00011033 and ”Unit of
Excellence María de Maeztu 2020-2023” award to the Institute of Cosmos
Sciences, Grant CEX2019-000918-M funded by MCIN/AEI/10.13039/501100011033. We
acknowledge financial support from the Generalitat de Catalunya (Grants
2021SGR01411 and 2021SGR01095).

## References

  * Van Camp  _et al._ (2017) M. Van Camp, O. de Viron, A. Watlet, B. Meurers, O. Francis, and C. Caudron, Reviews of Geophysics 55, 938 (2017).
  * Ryu  _et al._ (2007) C. Ryu, M. F. Andersen, P. Cladé, V. Natarajan, K. Helmerson, and W. D. Phillips, Phys. Rev. Lett. 99, 260401 (2007).
  * Heathcote  _et al._ (2008) W. Heathcote, E. Nugent, B. Sheard, and C. Foot, New Journal of Physics 10, 043012 (2008).
  * Ramanathan  _et al._ (2011) A. Ramanathan, K. C. Wright, S. R. Muniz, M. Zelan, W. T. Hill, C. J. Lobb, K. Helmerson, W. D. Phillips, and G. K. Campbell, Phys. Rev. Lett. 106, 130401 (2011).
  * Moulder  _et al._ (2012) S. Moulder, S. Beattie, R. P. Smith, N. Tammuz, and Z. Hadzibabic, Phys. Rev. A 86, 013629 (2012).
  * Pandey  _et al._ (2019) S. Pandey, H. Mas, G. Drougakis, P. Thekkeppatt, V. Bolpasi, G. Vasilakis, K. Poulios, and W. von Klitzing, Nature 570, 205 (2019).
  * Smerzi  _et al._ (1997) A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy, Phys. Rev. Lett. 79, 4950 (1997).
  * Fialko  _et al._ (2015) O. Fialko, B. Opanchuk, A. I. Sidorov, P. D. Drummond, and J. Brand, Europhysics Letters 110, 56001 (2015).
  * Pigneur and Schmiedmayer (2018) M. Pigneur and J. Schmiedmayer, Physical Review A 98, 063632 (2018).
  * Fialko  _et al._ (2012) O. Fialko, M.-C. Delattre, J. Brand, and A. R. Kolovsky, Phys. Rev. Lett. 108, 250402 (2012).
  * Gilary  _et al._ (2003) I. Gilary, N. Moiseyev, S. Rahav, and S. Fishman, Journal of Physics A: Mathematical and General 36, L409 (2003).
  * Jiang  _et al._ (2023) J. Jiang, E. Bernhart, M. Röhrle, J. Benary, M. Beck, C. Baals, and H. Ott, Phys. Rev. Lett. 131, 033401 (2023).
  * He and Liu (2023) W. He and C.-Y. Liu, Annals of Physics , 169218 (2023).
  * Richards  _et al._ (2018) C. J. Richards, T. J. Smart, P. H. Jones, and D. Cubero, Scientific reports 8, 13107 (2018).
  * Condon (1928) E. U. Condon, Physical Review 31, 891 (1928).
  * Doncheski and Robinett (2003) M. Doncheski and R. Robinett, Annals of Physics 308, 578 (2003).
  * (17) DLMF, “NIST Digital Library of Mathematical Functions,” http://dlmf.nist.gov/, Release 1.0.28 of 2020-09-15 (2023), f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders,H. S.
  * Muñoz Mateo  _et al._ (2019) A. Muñoz Mateo, V. Delgado, M. Guilleumas, R. Mayol, and J. Brand, Phys. Rev. A 99, 023630 (2019).
  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Morsch and Oberthaler (2006) O. Morsch and M. Oberthaler, Rev. Mod. Phys. 78, 179 (2006).
  * Kittel and McEuen (2018) C. Kittel and P. McEuen,  _Introduction to solid state physics_ (John Wiley & Sons, 2018).
  * Olshanii (1998) M. Olshanii, Phys. Rev. Lett. 81, 938 (1998).
  * Louis  _et al._ (2003) P. J. Louis, E. A. Ostrovskaya, C. M. Savage, and Y. S. Kivshar, Physical Review A 67, 013602 (2003).
  * Pitaevskii and Stringari (2003) L. P. Pitaevskii and S. Stringari,  _Bose-Einstein Condensation_ , International Series of Monographs on Physics (Oxford University Press, Oxford, New York, 2003).
  * Kapitza (1951) P. Kapitza, JETP 21, 588 (1951).
  * Astrakharchik and Astrakharchik (2011) G. E. Astrakharchik and N. A. Astrakharchik, “Numerical study of kapitza pendulum,” (2011), arXiv:1103.5981 [cond-mat.other] .
  * Landau and Lifshitz (1960) L. D. Landau and E. M. Lifshitz,  _Mechanics_ , Vol. 1 (CUP Archive, 1960).
  * Dobrek  _et al._ (1999) L. Dobrek, M. Gajda, M. Lewenstein, K. Sengstock, G. Birkl, and W. Ertmer, Phys. Rev. A 60, R3381 (1999).
  * Zumberge (2021) M. A. Zumberge, Encyclopedia of Solid Earth Geophysics , 633 (2021).
  * Debs  _et al._ (2011) J. Debs, P. Altin, T. Barter, D. Doering, G. Dennis, G. McDonald, R. Anderson, J. Close, and N. Robins, Physical Review A 84, 033610 (2011).
  * Ménoret  _et al._ (2018) V. Ménoret, P. Vermeulen, N. Le Moigne, S. Bonvalot, P. Bouyer, A. Landragin, and B. Desruelle, Scientific reports 8, 12300 (2018).
  * Stray  _et al._ (2022) B. Stray, A. Lamb, A. Kaushik, J. Vovrosh, A. Rodgers, J. Winch, F. Hayati, D. Boddice, A. Stabrawa, A. Niggebaum, _et al._ , Nature 602, 590 (2022).
  * Kasevich and Chu (1991) M. Kasevich and S. Chu, Phys. Rev. Lett. 67, 181 (1991).
  * Nguyen  _et al._ (2017) J. H. Nguyen, D. Luo, and R. G. Hulet, Science 356, 422 (2017).
  * Pérez-García  _et al._ (1997) V. M. Pérez-García, H. Michinel, J. I. Cirac, M. Lewenstein, and P. Zoller, Phys. Rev. A 56, 1424 (1997).
  * Pérez-García  _et al._ (1998) V. M. Pérez-García, H. Michinel, and H. Herrero, Phys. Rev. A 57, 3837 (1998).
  * Beattie  _et al._ (2013) S. Beattie, S. Moulder, R. J. Fletcher, and Z. Hadzibabic, Physical review letters 110, 025301 (2013).
  * Fang  _et al._ (2014) B. Fang, G. Carleo, A. Johnson, and I. Bouchoule, Phys. Rev. Lett. 113, 035301 (2014).
  * Huang  _et al._ (2019) B. Huang, I. Fritsche, R. S. Lous, C. Baroni, J. T. M. Walraven, E. Kirilov, and R. Grimm, Phys. Rev. A 99, 041602 (2019).
  * Poynting and Thomson (1908) J. H. Poynting and J. J. Thomson,  _A Text-book of Physics_ (C. Griffin, 1908).
  * Symon (1971) K. Symon,  _Mechanics_, Addison-Wesley series in physics (Addison-Wesley Publishing Company, 1971).
  * Tabor (1989) M. Tabor,  _Chaos and Integrability in Nonlinear Dynamics: An Introduction_ (John Wiley and Sons, 1989).
  * de Goër de Herve  _et al._ (2021) M. de Goër de Herve, Y. Guo, C. D. Rossi, A. Kumar, T. Badr, R. Dubessy, L. Longchambon, and H. Perrin, Journal of Physics B: Atomic, Molecular and Optical Physics 54, 125302 (2021).
  * Bell  _et al._ (2016) T. A. Bell, J. A. P. Glidden, L. Humbert, M. W. J. Bromley, S. A. Haine, M. J. Davis, T. W. Neely, M. A. Baker, and H. Rubinsztein-Dunlop, New Journal of Physics 18, 035003 (2016).