30

11institutetext: Artur Palha, Pedro Pinto Rebelo and Marc Gerritsma
22institutetext: Delft University of Technology, Faculty of Aerospace
Engineering, Aerodynamics Group, 2600GB - Delft,The Netherlands, 22email:
{A.Palha, P.J.PintoRebelo<EMAIL_ADDRESS>
# Mimetic Spectral Element Advection

Artur Palha  Pedro Pinto Rebelo and Marc Gerritsma

###### Abstract

We present a discretization of the linear advection of differential forms on
bounded domains. The framework established in kreeft2011mimetic  is extended
to incorporate the Lie derivative, $\mathcal{L}$, by means of Cartan’s
homotopy formula. The method is based on a physics-compatible discretization
with spectral accuracy . It will be shown that the derived scheme has spectral
convergence with local mass conservation. Artificial dispersion depends on the
order of time integration.

##  1 INTRODUCTION

Consider the classical advection problem for a scalar function in conservation
form,

$\frac{\partial\rho}{\partial t}+\nabla\cdot\left(\vec{v}\rho\right)=0,$ (1)

where $\vec{v}$ is a prescribed uniformly Lipschitz continuous vector field
and $\rho$ the advected scalar function. The method presented in this work is
based on the approximation of the differential operators with the focus on the
spatial discretization and a time-stepping scheme that distinguishes between
quantities evaluated at time instants and quantities evaluated over time
intervals.

The mimetic framework, presented in kreeft2011mimetic , showed that using the
differential geometric approach for the representation of physical laws
clarifies the underlying structures. One clearly identifies to what kind of
geometrical object a certain physical quantity is associated and this
determines how its discretization must be done (e.g.: evaluation at points,
integration over lines, surfaces or volumes). Additionally, a well defined,
metric free, representation of differential operators is obtained, together
with their metric dependent Hilbert adjoints. For these reasons, the authors
followed this approach for the advection equation. It is known, see
(abraham_diff_geom, , pp. 317), that (1) is a particular case of the
generalized advection equation which can be written in terms of differential
geometry as,

$\frac{\partial\alpha^{(k)}}{\partial
t}+\mathcal{L}_{\vec{v}}\,\alpha^{(k)}=0.$ (2)

The advection operator, $\mathcal{L}_{\vec{v}}$, is the Lie derivative for the
prescribed velocity field $\vec{v}$ and the advected quantity is given by the
$k$-differential form $\alpha^{(k)}$. Depending on the index $k$ the quantity
$\alpha^{(k)}$ can represent scalar, vector and higher dimensional quantities.

##  2 DIFFERENTIAL GEOMETRY

In this section a brief introduction to differential geometry is given. For a
more detailed introduction the reader is directed to abraham_diff_geom . Given
an $n$-dimensional smooth orientable manifold $\Omega$ it is possible to
define in each point a tangent vector space $E$ of dimension $n$. The space of
smooth vector fields on a manifold is the space, $\Gamma$, of smooth
assignments of elements of $E$ to each point of the manifold. We denote by
$\Lambda^{k}$, $k$ an integer $0\leq k\leq n$, the space of differentiable
$k$-forms, i.e. all smooth $k$-linear, antisymmetric maps
$\omega^{(k)}:E\times\cdots\times E\rightarrow\mathbb{R}$, at every point of
the manifold. We recall the wedge product
$\wedge:\Lambda^{k}\times\Lambda^{l}\rightarrow\Lambda^{k+l}$ for $k+l\leq n$
with the property that
$\alpha^{(k)}\wedge\beta^{(l)}=(-1)^{kl}\beta^{(l)}\wedge\alpha^{(k)}$. The
inner product $\left(\cdot,\cdot\right)$ on $E$ induces at each point of the
manifold an inner product $\left(\cdot,\cdot\right)$ on $\Lambda^{1}$. In
turn, this can be extended to a local inner product on $\Lambda^{k}$,
(morita::geometry_diff_forms, , pp. 149). The local inner product gives rise
to a unique metric operator, Hodge-$\star$,
$\star:\Lambda^{k}\rightarrow\Lambda^{n-k}$, defined by
$\alpha^{(k)}\wedge\star\beta^{(k)}=\left(\alpha^{(k)},\beta^{(k)}\right)\omega^{(n)}$,
where $\omega^{(n)}=\star 1$ is the standard volume form. By integration, one
can define an inner product on $\Omega$ as
$\left(\cdot,\cdot\right)_{L^{2}}:=\int_{\Omega}\left(\cdot,\cdot\right)\omega^{(n)}$.
The exterior derivative $\mathrm{d}:\Lambda^{k}\rightarrow\Lambda^{k+1}$
satisfies the following rule,
$\mathrm{d}\left(\alpha^{(k)}\wedge\beta^{(l)}\right)=\mathrm{d}\alpha^{(k)}\wedge\beta^{(l)}+(-1)^{k}\alpha^{(k)}\wedge\mathrm{d}\beta^{(l)}$
and by definition $\mathrm{d}\alpha^{(n)}=0$. The flat operator, $\flat$, is a
mapping $\flat:\Gamma\mapsto\Lambda^{1}$.

The Lie derivative along a tangent vector field, $\vec{v}$, is denoted by
$\mathcal{L}_{\vec{v}}$ and represents the advection operator in differential
geometry. It is a mapping
$\mathcal{L}_{\vec{v}}:\Lambda^{k}\mapsto\Lambda^{k}$. From Cartan’s homotopy
formula the Lie derivative can be written as

$\mathcal{L}_{\vec{v}}\,\alpha^{(k)}():=\mathrm{d}\iota_{\vec{v}}\alpha^{(k)}+\iota_{\vec{v}}\mathrm{d}\alpha^{(k)}\;,$

where the interior product of a tangent vector field, $\vec{v}$, with a
$k$-form, $\alpha^{(k)}$, is a mapping
$\iota_{\vec{v}}\alpha^{(k)}:\Lambda^{k}\rightarrow\Lambda^{k-1}$ given by:

$\iota_{\vec{v}}\alpha^{(k)}(\vec{X}_{2},\cdots,\vec{X}_{k}):=\alpha^{(k)}\left(\vec{v},\vec{X}_{2},\cdots,\vec{X}_{k}\right),\quad\forall\vec{X}_{i}\in\Gamma\quad\mathrm{and}\quad\iota_{\vec{v}}\alpha^{(0)}=0,\quad\forall\vec{v}\in\Gamma\;.$

The interior product is the adjoint of the wedge product, made explicit by:

$\left(\iota_{\vec{v}}\,\alpha^{(k)},\beta^{(k-1)}\right)_{L^{2}\Lambda^{k-1}}=\left(\alpha^{(k)},\vec{v}^{\flat}\wedge\beta^{(k-1)}\right)_{L^{2}\Lambda^{k}},\quad\forall\beta^{(k-1)}\in\Lambda^{k-1}$
(3)

where $\vec{v}^{\,\flat}=\nu^{(1)}\in\Lambda^{1}$ and
$\alpha^{(k)}_{h}\in\Lambda^{k}$.

The relevance of this adjoint relation between the interior product and the
wedge product lies in the fact that it shows how a physical quantity
represented by an interior product with a vector field can be represented by
its dual differential 1-form.

For a volume form $\rho^{(n)}$ the Lie derivative is simply
$\mathcal{L}_{\vec{v}}\,\rho^{(n)}=\mathrm{d}\iota_{\vec{v}}\rho^{(n)}$ and
for a $0$-form,
$\mathcal{L}_{\vec{v}}\,\alpha^{(0)}=\iota_{\vec{v}}\mathrm{d}\alpha^{(0)}$.

##  3 MIMETIC DISCRETIZATION

In this section a brief introduction to the discretization of physical
quantities and to the discretization of the exterior derivative is presented.
For a more detailed presentation the reader is directed to kreeft2011mimetic ;
palhaLaplaceDualGrid ; gerritsma::geometricBasis

Consider a three dimensional domain $\Omega$ and an associated grid consisting
of a collection of points, $\tau_{(0),i}$, line segments connecting the
points, $\tau_{(1),i}$, surfaces bounded by these line segments,
$\tau_{(2),i}$, and volumes bounded by these surfaces, $\tau_{(3),i}$.

Let $\Lambda^{k}$ be the space of smooth differentiable $k$-forms.
Additionally, let the finite dimensional space of differentiable forms be
defined as
$\Lambda^{k}_{h}=\mathrm{span}(\\{\epsilon_{i}^{(k)}\\}),\,\,i=1,\cdots,\dim(\Lambda^{k}_{h})$,
where $\epsilon_{i}^{(k)}\in\Lambda^{k}$ are basis $k$-forms. Under these
conditions it is possible, see kreeft2011mimetic ; palhaLaplaceDualGrid , to
define a projection operator $\pi_{h}$ which projects elements of
$\Lambda^{k}$ onto elements of $\Lambda^{k}_{h}$ which satisfies:

$\pi_{h}\mathrm{d}=\mathrm{d}\pi_{h}\;.$ (4)

It is possible to write:

$\pi_{h}\alpha^{(k)}=\alpha_{h}^{(k)}=\sum_{i}\alpha_{i}\epsilon_{i}^{(k)}\;,$

where

$\alpha_{i}=\int_{\tau_{(k),i}}\alpha^{(k)}\quad\mathrm{and}\quad\int_{\tau_{(k),i}}\epsilon_{j}^{(k)}=\delta_{ij},\quad
k=0,1,\cdots,n\;.$

A set of basis functions yielding a projection operator $\pi_{h}$ that
satisfies (4) can be constructed using piecewise polynomial expansions on the
quadrilateral elements using tensor products. Thus, it suffices to derive the
basis forms in one dimension on a reference interval and generalize them in
$n$ dimensions.

In one dimension take a 0-form,
$\alpha^{(0)}\in\Lambda^{0}\left(Q_{ref}\right)$, where
$Q_{ref}:=\left[-1,1\right]$. Define on $Q_{ref}$ a cell complex $D$ of order
$p$ consisting of $(p+1)$ nodes $\tau_{(0),i}=\xi_{i}$ with $i=0,\cdots,p$,
where $-1\leq\xi_{0}<\cdots<\xi_{i}<\cdots\xi_{p}\leq 1$ are the Gauss-Lobatto
quadrature nodes, and $p$ edges, $\tau_{(1),i}=\left[\xi_{i-1},\xi_{i}\right]$
with $i=1,\cdots,p$. The projection operator $\pi_{h}$ reads:

$\displaystyle\pi_{h}\alpha^{(0)}\left(\xi\right)=\sum_{i=0}^{p}\alpha_{i}\epsilon^{(0)}_{i}(\xi)\;,$
(5)

where $\epsilon^{(0)}_{i}(\xi)=l_{i}\left(\xi\right)$ are the $p^{th}$ order
_Lagrange polynomials_ and $\alpha_{i}=\alpha^{0}(\xi_{i})$. Similarly in one
dimension for the projection of 1-forms Gerritsma gerritsma::edge_basis  and
Robidoux robidoux-polynomial  derived 1-form polynomials called _edge
polynomials_ , $\epsilon^{(1)}_{i}\in\Lambda^{1}_{h}\left(Q_{ref}\right)$,

$\displaystyle\epsilon^{(1)}_{i}(\xi)=e_{i}\left(\xi\right)\mathrm{d}\xi,\quad\text{with}\quad
e_{i}(\xi)=-\sum_{k=0}^{i-1}\frac{dl_{k}}{d\xi}\;.$ (6)

Note that in this way we have:

$\displaystyle\int_{\xi_{j-1}}^{\xi_{j}}\epsilon^{(1)}_{i}=\int_{\xi_{j-1}}^{\xi_{j}}e_{i}\left(\xi\right)\mathrm{d}\xi=\delta_{ij}\;.$
(7)

Moreover, the exterior derivative of the basis 0-forms is given by:

$\mathrm{d}\epsilon_{i}^{(0)}=\frac{\mathrm{d}l_{i}}{\mathrm{d}\xi}\,\mathrm{d}\xi=-\sum_{k=0}^{i-1}\frac{\mathrm{d}l_{k}}{\mathrm{d}\xi}\,\mathrm{d}\xi-\left(-\sum_{k=0}^{i}\frac{\mathrm{d}l_{k}}{\mathrm{d}\xi}\,\mathrm{d}\xi\right)=\epsilon_{i}^{(1)}-\epsilon_{i+1}^{(1)},\quad
i=1,\cdots,p-1\;.$ (8)

In this way, the exterior derivative of a discrete 0-form can be written as:

$\mathrm{d}\alpha_{h}^{(0)}=\mathrm{d}\sum_{i=0}^{p}\alpha_{i}\epsilon_{i}^{(0)}=\sum_{i=0,j=1}^{p}\mathsf{E}^{(1,0)}_{ij}\alpha_{j}\epsilon_{i}^{(1)}\;,$
(9)

where, $\mathsf{E}^{(1,0)}_{ij}$ is the incidence matrix containing only the
values, 0, 1 and -1, see kreeft2011mimetic ; palhaLaplaceDualGrid ;
gerritsma::geometricBasis  for more details. This idea can be extended to
higher dimensions, giving rise to $k$-incidence matrices,
$\mathsf{E}^{(k+1,k)}_{ij}$, which represent the discrete exterior derivative
on discrete $k$-forms, see palhaLaplaceDualGrid ; gerritsma::geometricBasis .

##  4 Mimetic spectral advection: an application to 1D advection

In this section we want to illustrate how to discretize the advection
equation. Take the Lie advection of a 1-form,

$\frac{\partial\rho^{(1)}}{\partial
t}+\mathrm{d}\iota_{\vec{v}}\rho^{(1)}=0\Leftrightarrow\left\\{\begin{array}[]{l}\frac{\partial\rho^{(1)}}{\partial
t}=-\mathrm{d}\varsigma^{(0)}\\\ \\\
\iota_{\vec{v}}\,\rho^{(1)}=\varsigma^{(0)}\end{array}\;.\right.$ (10)

Here $\rho^{(1)}$ is the advected quantity, say mass density, and
$\varsigma^{(0)}$ represents the instantaneous fluxes of the advected quantity
under the vector field $\vec{v}$, which are discretized in space as,

$\varsigma_{h}^{(0)}=\sum_{i=0}^{p}\varsigma_{i}(t)\epsilon_{i}^{(0)}\quad\mathrm{and}\quad\rho_{h}^{(1)}=\sum_{i=1}^{p}\rho_{i}(t)\epsilon_{i}^{(1)}\;.$
(11)

For the sake of clarity in the method presentation we first introduce the time
treatment, then the interior product discretization and finally their
combination for a numerical solution of the advection problem.

###  4.1 Time integration

The time integrator used for solving the time evolution part of the advection
equation is the canonical mimetic one, an arbitrary order symplectic operator
derived in palha::phd , which is connected to canonical Gauss collocation
integrators. Take an ordinary differential equation of the unknown function
$y(t)$:

$\frac{\mathrm{d}y}{\mathrm{d}t}=h(y,t),\quad t\in I\subset\mathbb{R}\;.$ (12)

Discretizing $y(t)$ as $y_{h}=\sum_{k=0}^{p}y^{k}l_{k}(t)$, one gets:

$\frac{\mathrm{d}y_{h}}{\mathrm{d}t}=\sum_{k=0}^{p}y^{k}\frac{\mathrm{d}l_{k}(t)}{\mathrm{d}t}=\sum_{k=1}^{p}(y^{k}-y^{k-1})\,e_{k}(t)\;,$

where the superscript $k$ denotes the time level.

The approximated solution, $y_{h}(t)$, is a polynomial of order $p$ determined
by means of $(p+1)$ degrees of freedom such as its values at the Gauss-Lobatto
nodes, red dots in Figure 1. On the other hand,
$\frac{\mathrm{d}y_{h}}{\mathrm{d}t}$ is a polynomial of order $(p-1)$ defined
by only $p$ degrees of freedom. One can set these degrees of freedom to be the
values of the derivative in one point inside each of the $p$ intervals
$[t^{k},t^{k+1}]$. A choice that results in a symplectic integrator of order
$2p$ is to select these points as the Gauss nodes of order $(p-1)$, the blue
nodes of Figure 1. Notice that along the trajectory these nodes will not show
the usual Gauss-Lobatto and Gauss distribution patterns, since in general the
velocity field is not constant. In this way the discrete integrator becomes:

$\sum_{k=1}^{p}(y^{k}-y^{k-1})\,e_{k}(\tilde{t}^{q})=h\left(\sum_{k=0}^{p}y^{k}l^{k}(\tilde{t}^{q}),\tilde{t}^{q}\right),\quad
q=1,2,\cdots,p\;,$ (13)

with $\tilde{t}^{j}$ the $p$ nodes of a Gauss quadrature formula. The fact
that the instants in time, $t^{k}$, where the $y^{k}_{i}$ are defined
alternate with the instants in time, $\tilde{t}^{j}$, where the $h_{i}$ are
evaluated (see Figure 1), corresponds to a staggering in time. This staggering
also appears in leap-frog methods and in the implicit midpoint rule, for
instance.

Figure 1: Geometric interpretation of the solution of (12) as given by (13):
$(t,y^{(0)}(t))$. In red the Gauss-Lobatto nodes where the trajectory is
discretized. In blue, the Gauss nodes where its derivative is discretized. The
flow field, represented by arrows, is tangent to the curve at the Gauss nodes.
That is, the derivative of the approximate trajectory is exactly equal to the
flow field at the Gauss nodes.

The first equation in (10) using the discretization (11) and (9) can be
written as:

$\dfrac{\sum_{i}\mathrm{d}\rho_{i}(t)\epsilon_{i}^{(1)}}{\mathrm{d}t}=-\sum_{il}\mathsf{E}^{(1,0)}_{il}\varsigma_{l}(t)\epsilon_{i}^{(1)}\quad\Rightarrow\quad\dfrac{\mathrm{d}\rho_{i}(t)}{\mathrm{d}t}=-\sum_{l}\mathsf{E}^{(1,0)}_{il}\varsigma_{l}(t)\;.$
(14)

This equation has a similar form as (12), but now as a system of equations,
therefore one can apply the mimetic integrator, yielding:

$\sum_{k}(\rho_{i}^{k+1}-\rho_{i}^{k})\,e_{k}(\tilde{t}^{q})=-\sum_{l}\mathsf{E}^{(1,0)}_{il}\varsigma_{l}^{q}\;.$
(15)

Recall that $\rho_{i}^{k}$ is the discrete degree of freedom of the advected
quantity at the $t^{k}$ instants of time associated to Gauss-Lobatto nodes and
$\varsigma_{l}^{q}$ is the discrete degree of freedom of the fluxes of the
advected quantity at the $\tilde{t}^{q}$ instants of time associated to the
Gauss nodes, just as stated for the systems of ordinary differential
equations.

###  4.2 Interior product

The discretization of the interior product is done using (3), in the following
way:

######  Definition 1 (Discrete interior product)

In one dimension, the discrete interior product
$\iota_{\vec{v},h}:\Lambda^{1}_{h}\rightarrow\Lambda^{0}_{h}$ is such that:

$\left(\iota_{\vec{v},h}\,\alpha^{(1)}_{h},\epsilon^{(0)}_{i}\right)_{L^{2}}=\left(\alpha^{(1)}_{h},\vec{v}^{\flat}\wedge\epsilon^{(0)}_{i}\right)_{L^{2}},\quad\forall\epsilon_{i}^{(0)}\in\Lambda^{0}_{h}$
(16)

where $\vec{v}^{\,\flat}=\nu^{(1)}\in\Lambda^{1}$ and
$\alpha^{(1)}_{h}\in\Lambda^{1}_{h}$.

In this way one satisfies the duality pairing between the interior product and
the wedge product in the discrete setting.

Partitioning the domain $\Omega$ in a spectral element cell complex one can
apply the discretization of the interior product in each spectral element,
obtaining:

$\sum_{i}\rho_{i}(t)\left(\epsilon_{i}^{(1)},\nu^{(1)}\wedge\epsilon_{j}^{(0)}\right)_{L^{2}}=\sum_{i}\varsigma_{i}(t)\left(\epsilon_{i}^{(0)},\epsilon_{j}^{(0)}\right)_{L^{2}}\;,\quad\forall\epsilon_{j}^{(0)}\in\Lambda^{0}_{h}\;.$
(17)

###  4.3 Putting things together: advection

The complete discrete systems becomes:

$\begin{cases}\sum_{k}(\rho_{i}^{k+1}-\rho_{i}^{k})\tilde{e}_{k}(\tilde{t}^{q})=-\sum_{l}\mathsf{E}^{(1,0)}_{il}\varsigma_{l}^{q}\\\
\\\
\sum_{i,k}\rho_{i}^{k}\epsilon_{k}^{(0)}(\tilde{t}^{q})\left(\epsilon_{i}^{(1)},\star\nu^{(0)}\wedge\epsilon_{j}^{(1)}\right)_{L^{2}\Lambda^{1}(\Omega_{m})}=\sum_{i}\varsigma_{i}^{q}\left(\epsilon_{i}^{(0)},\epsilon_{j}^{(0)}\right)_{L^{2}\Lambda^{0}(\Omega_{m})},\,\forall\epsilon_{j}^{(0)}\in\Lambda^{0}(\Omega_{m})\end{cases}$
(18)

##  5 Numerical results

This approach was applied to the two dimensional solution of an advected sine
wave and a sine bell in a constant velocity field $\vec{v}=\vec{e}_{x}$:
$\rho^{(2)}(x,y)=\sin(\pi x)\sin(\pi y)\mathrm{d}x\mathrm{d}y$ (sine wave) and
$\rho^{(2)}(x,y)=\sin(2\pi x)\sin(2\pi y)\mathrm{d}x\mathrm{d}y$ if
$(x,y)\in[0,0.5]\times[0,0.5]$ and $\rho^{(2)}(x,y)=0$ in
$(x,y)\in\mathbb{R}^{2}\backslash[0,0.5]\times[0,0.5]$ (sine bell), on a
domain with periodic boundary conditions.

In Figure 2 the error in time of the numerical solution of (10) for a mesh of
$4\times 4$ elements with a $\Delta t=0.1s$ and various polynomial orders in
space, $p$, and time, $p_{t}$, is presented. The initial error, due to the
discretization, is conserved, as long as the time integration is sufficiently
accurate.

Figure 2: Error in time of the numerical solution of (10) with
$\vec{v}=\vec{e}_{x}$ and $4\times 4$ elements of order $p=3$, $p=6$, $p=10$
and $p=12$ (from left to right and top to bottom) and $\Delta t=0.1s$, for the
sine wave $\rho^{(2)}(x,y)=\sin(\pi x)\sin(\pi y)\mathrm{d}x\mathrm{d}y$. As
shown, the error in the solution increases with time due to the inaccuracy of
time integration. When time integration is accurate enough the error in the
initial state is preserved. Here $p_{t}$ denotes the polynomial degree in
time.

In Figure 3, the $h$\- and $p$-convergence plots are shown for different
values of the order of the time integration scheme, $p_{t}$, and $\Delta
t=0.1s$. It is possible to see that the method presents algebraic
$h$-convergence rates of order $(p+1)$ as long as the time integration error
does not dominate the spatial one. The method shows a spectral $p$-convergence
as soon as the time integration is accurate enough.

Figure 3: Left: $h$ convergence in space for the advection of a sine wave with
$\Delta t=0.1s$. Right: $p$ convergence in space for the advection of a sine
wave, $4\times 4$ elements and $\Delta t=0.1s$.

In Figure 4 the error on the velocity is presented as a function of the
advected sine wave frequency. This figure shows that the numerical method
introduces an artificial dispersion if the time scheme is not accurate enough.

Figure 4: Error in velocity as a function of the frequency of the advected
sine wave: numerical dispersion. $p=10$, $\Delta t=0.1s$ and $n=4\times 4$
elements.

Another fundamental aspect is the conservation of the advected quantity.
Figure 5 shows the mass error in time, that is:
$\int_{\Omega}\rho^{(2)}_{t}-\int_{\Omega}\rho^{(2)}_{t_{0}}$. The error goes
from the zero machine in the first $10^{3}$ time steps while thereafter it
steadily increases. Notice that even after $2\times 10^{4}$ time steps the
error is still below $10^{-12}$.

Figure 5: Sum of the local errors of the advected 2-form for a sine bell in a
velocity field $\vec{v}=\vec{e}_{x}$, with $50\times 50$ elements of order
$p=0$ (blue) and $4\times 4$ elements of order $p=9$ (red), $\Delta t=0.01$
and $p_{t}=2$.

In Figure 6 one can see the advection of a sine wave of frequency $\omega=\pi$
in a Rudman vortex for 100 time steps after which the direction of the flow is
reversed and the calculation is continued for another 100 time steps, with
$4\times 4$ curved elements of order $p=9$, $\Delta t=0.1s$ and time
integration of order $p_{t}=2$, on a distorted mesh. The mimetic advection
enables one to recover the initial solution, thus demonstrating that the
integration method is reversible.

Figure 6: From left to right and from top to bottom: advection of a sine wave
of frequency $\omega=\pi$ on a Rudman vortex with $4\times 4$ curved elements
of order $p=9$, $\Delta t=0.1s$ and time integration of order $p_{t}=2$, on a
distorted mesh. At time $t=10$, the flow field is reversed. At times $t=8.0s$
and $t=12.0s$ the mesh is visible in the solution. See
http://www.youtube.com/watch?v=QmoJyqtk9YA for animation.

###### Acknowledgements.

The authors would like to thank the valuable comments of both reviewers and
the funding received by FCT - Foundation for science and technology Portugal
through SRF/BD/36093/2007 and SFRH/BD/79866/2011.

## References

  * [1] R. Abraham, J. E. Marsden, and T. Ratiu.  Manifolds, Tensor Analysis, and Applications, volume 75 of Applied Mathematical Sciences.  Springer, 2001.
  * [2] M. I. Gerritsma, R. Hiemstra, J. J. Kreeft, A. Palha, P. Rebelo, and D. Toshniwal.  The geometric basis of numerical methods.  Proceedings of ICOSAHOM 2012 (this issue).
  * [3] Marc Gerritsma.  Edge functions for spectral element methods.  In Jan S. Hesthaven and Einar M. Rønquist, editors, Spectral and High Order Methods for Partial Differential Equations, volume 76 of Lecture Notes in Computational Science and Engineering, pages 199–207. Springer Berlin Heidelberg, 2011.
  * [4] J. Kreeft, A. Palha, and M. I. Gerritsma.  Mimetic framework on curvilinear quadrilaterals of arbitrary order.  Arxiv preprint arXiv:1111.4304, 2011.
  * [5] A. Palha.  High order mimetic discretization: Development and application to Laplace and convection-diffusion problems in arbitrary quadrilaterals.  PhD thesis, TUDelft, 2013.
  * [6] A. Palha, P. Pinto Rebelo, R. Hiemstra, J. Kreeft, and M. I. Gerritsma.  Physics-compatible discretization techniques on single and dual grids, with application to the poisson equation of volume forms.  Submitted to J. Comp Phys., 2012.
  * [7] N. Robidoux.  Polynomial Histopolation, Superconvergent Degrees Of Freedom, And Pseudospectral Discrete Hodge Operators.  Unpublished: http://www.cs.laurentian.ca/nrobidoux/prints/super/histogram.pdf, 2008.
  * [8] Morita Shigeyuki.  Geometry of differential forms, volume 201 of Translations of mathematical monographs.  American mathematical society, 2001.


# A Fast Online Omnidirectional Quadrupedal Jumping Framework Via Virtual-
Model Control and Minimum Jerk Trajectory Generation

Linzhu Yue1, Lingwei Zhang1, Zhitao Song1, Hongbo Zhang1,
Jinhu Dong1, Xuanqi Zeng1, and Yun-Hui Liu1 1 L. Z. Yue, L. W. Zhang, Z. T.
Song, H. B. Zhang, J. H. Dong, X. Q. Zeng, and Y.-H. Liu is with the
Department of Mechanical and Automation Engineering at the Chinese University
of Hong Kong<EMAIL_ADDRESS>Corresponding author: Y.-H. Liu
<EMAIL_ADDRESS>work is supported by the InnoHK of the Government of
Hong Kong via the Hong Kong Centre for Logistics Robotics, the CUHK T Stone
Robotics Institute, and the Shenzhen Portion of Shenzhen-Hong Kong Science and
Technology Innovation Cooperation Zone under HZQB-KCZYB-20200089.

###### Abstract

Exploring the limits of quadruped robot agility, particularly in the context
of rapid and real-time planning and execution of omnidirectional jump
trajectories, presents significant challenges due to the complex dynamics
involved, especially when considering significant impulse contacts. This paper
introduces a new framework to enable fast, omnidirectional jumping
capabilities for quadruped robots. Utilizing minimum jerk technology, the
proposed framework efficiently generates jump trajectories that exploit its
analytical solutions, ensuring numerical stability and dynamic compatibility
with minimal computational resources. The virtual model control is employed to
formulate a Quadratic Programming (QP) optimization problem to accurately
track the Center of Mass (CoM) trajectories during the jump phase. The whole-
body control strategies facilitate precise and compliant landing motion.
Moreover, the different jumping phase is triggered by time-schedule. The
framework’s efficacy is demonstrated through its implementation on an enhanced
version of the open-source Mini Cheetah robot. Omnidirectional jumps—including
forward, backward, and other directional—were successfully executed,
showcasing the robot’s capability to perform rapid and consecutive jumps with
an average trajectory generation and tracking solution time of merely 50
microseconds.

##  I INTRODUCTION

Robotics jumping is critical for quadruped robots to traverse between
complicated unstructured terrains. Although quasi-static jumping algorithms
empower quadrupedal robots to cross mild irregular terrain, when aimed at
avoiding unknown obstacles, it is necessary for robots to perform
omnidirectional jumping promptly, which requires fast, skillful planning and
dynamical executing motions. We strive to propose a framework with minimal
computational resources that enables quadrupedal robots to make rapid
omnidirectional jump responses to avoid sudden dangers while exploring unknown
environments. This could be significantly helpful in inspection scenarios
where quadrupedal robots are equipped with large, heavy instruments and the
onboard computer has other heavy-load tasks. The real-time omnidirectional
jumping framework could improve the robot’s ability to traverse rough terrains
and ensure the safety of robots.

Figure 1: Various jumping motion experiments to validate the proposed
omnidirectional jumping framework. (a) Front-left jumping consecutively. (b)
Single rear-right jumping. (c) Single rear jumping. (d) Single front jumping.
The yellow line shows the trajectory of CoM, and the red dashed line
represents the trajectory of the selected foot. The red dots indicate the
foot’s contact with the ground, and the orange dots indicate the foot in the
flight phase.

Various experimental conductions demonstrate the great potential of optimal
control to provoke robots to achieve robust jumping motions. Robots such as
ANYmal[1], SALTO-1P[2] and Mini Cheetah[3] can perform remarkable jumping
motions and navigating in the local environment. A traditional method
simplifies the robot model as a 2-D planar model by combining left and right
legs to accelerate reference trajectory optimization speed. Model Predictive
Control (MPC)[4], Whole-Body Impulse Control (WBIC)[5], and Virtual Model
Control (VMC)[6] are implemented to track the reference trajectory. The MIT
Cheetah 2 is able to autonomously jump over obstacles up to 40 cm in height
during bounding gait through Nonlinear MPC[7]. Except for single rigid body
models, with kino-dynamic and novel torque-speed limitation constraints, the
Mini Cheetah can reliably produce successful aerial motions such as flips and
barrel rolls[8]. The cost on the order of seconds rather than milliseconds
confines the performance of the kino-dynamic planners to generate a motion
plan, which implies the real-time limitation of running jumping that starts
from the non-static initial position. Though our previous work[9],[10] based
on offline Differential Evolution (DE) and online evolutionary algorithms with
pre-motion library also endow quadrupedal robots with the ability to execute
complicated jumping motions, the trajectory optimization result is not
reliable in some extreme cases.

Reinforcement Learning (RL) is another broad way to achieve accurate and
aggressive quadrupedal jumping motions. In [11], with policies considering the
robot’s total power limits and torque-speed relationships, Unitree A1[12]
achieves aggressive and accurate jumping motions. Cat-like jumping exploits
the possibility of adjusting robot body gestures by jumping motions in a low
gravity environment [13]. However, few researchers attach great importance to
using a single policy to make robots perform complicated or omnidirectional
jumping motions.

For real-time omnidirectional jumping frameworks controlled by a
microcontroller unit, the robot Moobot[14] is capable of traversing different
platforms from all directions. However, the insect-scale robot cannot jump
consecutively and must be reloaded by hand after each jump. As for quadruped
robots, a hierarchical planning and control framework[15] is proposed to
enable Mini Cheetah to traverse complex multi-layered terrain. A novel high-
level jump selection controller is implemented to ensure the most robustness
guarantees. In our work, enabling quadruped robots to perform omnidirectional
jumping consecutively, we proposed a framework that minimizes the computing
time of the trajectory planner and the reference tracking controller to
improve the agility of quadruped robots significantly. Due to analytic
solutions, the minimum-jerk trajectory planner efficiently generates a CoM
reference path in the jumping phase with numerical stability and dynamic
compatibility. Aiming to track reference motions, an intuitive virtual model
controller is deployed to compute the force of feet with low computational
resources. Then, in the landing phase, a whole-body controller guarantees the
CoM inside the support polygon of four feet, ensuring the stability and
compatibility of the robot’s body.

This paper makes the following contributions:

  * •

We propose an omnidirectional jumping framework that generates and tracks
aerial motions for quadrupedal robots, made up of a minimum jerk CoM
trajectory planner, a VMC reference tracking controller, and a WBC landing
controller.

  * •

The average cost of generating trajectory and computing tracking motor torque
commands to perform omnidirectional jumping motion is within 50 us, implying
the remarkable real-time performance of the framework.

  * •

The efficacy of the framework is verified on open-source Mini Cheetah[24], and
the robot succeeds in generating multiple reference CoM trajectories and
performing omnidirectional jumping behaviors consecutively (see Fig. 1)

The remaining content of this paper is structured as follows. Section II
briefly introduces the models of robot and omnidirectional jumping. Section
III details the formulation of the minimum jerk trajectory planner and the VMC
controller, and section IV is the implementation details and experiments,
including the hardware setup, real-time performance, and, verification of the
omnidirectional jumping.

##  II Models and dynamics

###  II.A Robot Model

The reduced-order dynamic model of jumping motion treats the robot as a single
rigid body (SRB) with a specified moment of inertia. The robot state $\bm{x}$
can be written as:

$\displaystyle\bm{x}:=[{\bm{P}_{com}^{T}}\quad{\bm{\Theta}^{T}}\quad{\bm{V}_{C}^{T}}\quad{\bm{\omega}^{T}_{B}}]^{T}\in\mathbb{R}^{12}$
(1a)
$\displaystyle{\bm{Q}}:=[{\bm{q}_{i}}\quad{\dot{\bm{q}}_{i}}]\in\mathbb{R}^{24}$
(1b)

Where ${\bm{P}_{com}}\in\mathbb{R}^{3}$ is the position of the robot’s body
center of mass (CoM) with respect to (w.r.t.) inertial frame (see Fig.2);
${\bm{\Theta}=[\psi\ \phi\ \theta]}$ represents the Euler angles of the
robot’s body; ${\bm{V}_{C}\in\mathbb{R}^{3}}$ is the velocity of the CoM.
${\bm{\omega}}_{B}\in\mathbb{R}^{3}$ is the angular velocity of CoM w.r.t. the
robot frame $\\{B\\}$. $\bm{q}_{i}\in\mathbb{R}^{3}$ and
$\bm{\dot{}}{\bm{q}}_{i}\in\mathbb{R}^{3}$ are the joint angles and velocities
of each leg. $i$ is the number of feet. The ground reaction force (GRF)
$\bm{f}_{i}\in\mathbb{R}^{3}$ at each contact point consists of the dynamic
system control input ${\bm{u}}:=[{\bm{f}_{i}}]$. ${\bm{r}_{i}}$ is the vector
from CoM to the robot foot. Then, The linear acceleration of CoM and the
angular acceleration of the base are shown:

$\displaystyle m\ddot{\bm{P}}_{com}$
$\displaystyle={\sum_{i=1}^{n_{c}}\bm{f}_{i}}-{\bm{g}}$ (2a)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}({\bm{I}_{B}\bm{\omega}_{B}})$
$\displaystyle=\sum_{1}^{n_{c}}\bm{r}_{i}\times\bm{f}_{i}+{\bm{\omega}_{B}}\times{\bm{I}}_{B}{\bm{\omega}}_{B},$
(2b)

where ${\bm{g}}\in\mathbb{R}^{3}$ represents gravitational acceleration.
$n_{c}$ represents the number of contacts.
${{}^{B}\bm{I}}\in\mathbb{R}^{3\times 3}$ is the robot’s rotation inertial
tensor, which is assumed as a constant in this work,
$\textrm{diag}(^{B}\bm{I})=[0.07,0.26,0.242]^{T}$.

Figure 2: A model of a single rigid body (SRB) utilized in the framework for
VMC. The blue arrow represents the CoM to the plantar position vector, while
the red arrow represents the Ground Reaction Forces (GRFs).

###  II.B Omnidirectional Jumping Model

Omnidirectional jumping significantly increases quadruped robots’ capability
to access terrains instead of traditional assumptions that robots jump in 2D.
The motion of omnidirectional jumping comprises three phases: preparing,
flight, and landing (see Fig. 3). In the preparing phase, given desired
$[\bm{P}_{CoM}\ \bm{\dot{P}}_{CoM}\ \bm{\ddot{P}}_{CoM}]$, mini jerk planner
and motor torques generate a smooth CoM trajectory are calculated the motor
torque commands $\tau_{cmd}$ to track the reference path. In the flight phase,
a proportional and derivative controller ensures the preparation of the robot
landing. In the landing phase, given reference GRFs and desired base angles, a
whole-body controller maintains the stability of the robot’s base regardless
of the impact of the ground on the robot’s legs.

Figure 3: The Jumping phase, Flight phase, and Landing phase consist of the
omnidirectional jumping motion.

##  III Framework

This section introduces the online omnidirectional jumping framework with low
computational requirements, comprising three primary components: a minimum
jerk trajectory generator, a virtual model controller, and a whole-body
landing controller. With an analytic solution of relative coefficients, the
five-order polynomial trajectory generator rapidly generates a smooth and
dynamic feasible path for the CoM of robots based on its current states.
Through optimizing output GRFs, the VMC enables the robot to track the CoM
reference trajectory using a quadratic program and a Cartesian PD controller.
Finally, given reference GRFs and reference base Euler angles, the whole-body
landing controller ensures the stability of the robot’s landing, and a PD
motor torque controller is implemented to improve the robot’s performance of
tracking references in the landing phase. An overview is shown in Fig 4. The
first subsection shows details of the mini-jerk trajectory planner, the second
section reveals the process of generating motor torque commands by VMC, and
the last subsection briefly introduces the WBC and motor torque control
methods in the landing phase.

Figure 4: The overview of the fast online omnidirectional quadrupedal jumping
framework, made of 3 primary components.

###  III.A Minimum Jerk Trajectory Planner

The minimum jerk/snap algorithm has been extensively utilized in UAV
systems[16]. UAVs possess unique dynamic features that satisfy the criteria of
differential flatness, demonstrating that UAV control can be used for
trajectory generation. When the UAV system employs mini-jerk trajectory
creation, the UAV experiences minimal thrust, while using mini-snap results in
the minimum differential thrust. For the quadrupedal robot, considering the
SRB model, when all four feet are in contact with the ground, only the force
in the Z direction and the acceleration of the CoM are taken into account,
similar to trajectory generation in a UAV system. The system with a fixed
inclination satisfies the differential flatness property even if the force
direction does not align perfectly with the vector leading to the center of
mass along the sole of the quadruped robot[17].

Currently, the trajectory produced by the mini-jerk allows for the generation
of relatively reasonable GRFs in a short amount of time. (2a) the connection
between the CoM acceleration and the contact force helps create a smooth
trajectory for the quadruped robot’s jumping motion. This trajectory is then
sent to the VMC for tracking and optimization to achieve a 12-dimensional
contact force that meets kino-dynamics requirements. In this work, reference
Euler angles and angular velocities are selected as hyperparameters, and the
trajectory outputs are chosen as:

$\displaystyle s(t)=[x\quad y\quad z]=\bm{\psi}$ (3)

To make a quadruped robot well track the space the flatness outputs, the
polynomial order of the smooth trajectory is selected as five, and it is
convenient to formulate it as three-piecewise segments:

$\displaystyle s(t)$ $\displaystyle=\left\\{\begin{aligned}
s_{1}(t)=\sum&{}_{i=0}^{5}a_{1,i}t^{i}\quad T_{0}\leq t\leq T_{1}\\\
s_{2}(t)=\sum&{}_{i=0}^{5}a_{2,i}t^{i}\quad T_{1}\leq t\leq T_{2}\\\
s_{3}(t)=\sum&{}_{i=0}^{5}a_{3,i}t^{i}\quad T_{2}\leq t\leq
T_{3}\end{aligned}\right.$ (4a) $\displaystyle
s_{j}(t)^{{}^{\prime\prime\prime}}$ $\displaystyle=\sum_{i\geq
3}^{5}i(i-1)(i-2)t^{i-3}a_{j,i}$ (4b)

Where $s(t)$ represents the whole trajectory of CoM, $s_{j}(t)$ is each
segment’s trajectory. $a_{j,i}$ is the coefficient of each 5-order polynomial
segment path. $s_{j}(t)^{{}^{\prime\prime\prime}}$ is the jerk of each
segment.

Under the precondition of tracking reference trajectory, smaller reference GRF
values in (2a) ensure motors’ accuracy of torque and speed. Therefore, the
optimization program of figuring out a trajectory equipped with the minimal
acceleration reference can be formulated as follows:

$\displaystyle J_{j}(T)$ $\displaystyle=$
$\displaystyle\int_{T_{j-1}}^{T_{j}}(s_{j}(t)^{{}^{\prime\prime\prime}})^{2}dt=\bm{p}_{j}^{T}\bm{Q}_{j}\bm{p}_{j}$
(5a) $\displaystyle s.t.\quad s_{j}^{k}(T_{j})$ $\displaystyle=$
$\displaystyle\bm{\psi}_{T,j}^{k}$ (5b) $\displaystyle s_{j}^{k}(T_{j})$
$\displaystyle=$ $\displaystyle s_{j+1}^{k}(T_{j})$ (5c)

where (5b) is the derivative constraint for one polynomial segment, (5c) is
the continuity constraint between two segments, and k is the derivative order.
Then the solution of this problem can be solved in closed form ref to [18]:

$\displaystyle J$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{c}\bm{a}_{f}\\\
\bm{a}_{p}\end{array}\right]^{T}\underbrace{\left[\begin{array}[]{cc}\bm{R}_{ff}&\bm{R}_{fp}\\\
\bm{R}_{pf}&\bm{R}_{pp}\end{array}\right]}_{\bm{R}}\left[\begin{array}[]{c}\bm{a}_{f}\\\
\bm{a}_{p}\end{array}\right]$ (6g) $\displaystyle=$
$\displaystyle\bm{a}_{f}^{T}\bm{R}_{ff}\bm{a}_{f}+\bm{a}_{f}^{T}\bm{R}_{fp}\bm{d}_{p}+\bm{a}_{p}^{T}\bm{R}_{pf}\bm{a}_{f}+\bm{a}_{p}^{T}\bm{R}_{pp}\bm{a}_{p}$
(6h)

where $\bm{a}=\bm{C}^{T}\left[\bm{a}_{f}\quad\bm{a}_{p}\right]^{T}$, $\bm{C}$
is the selecting matrix that divides $\bm{a}$ to free $\bm{a}_{f}$ and
constrained $\bm{a}_{p}$.
$\bm{R}=\bm{C}\bm{M}^{-T}\bm{Q}\bm{M}^{-1}\bm{C}^{T}$, $\bm{M}$ is the
decision variable mapping matrix. The derivative of (6h) gives the solution:
$\bm{a}_{p}^{*}=-\bm{R}_{pp}^{-1}\bm{R}^{T}_{fp}\bm{a}_{f}$.

Given a final position target w.r.t body frame, piecewise trajectories depend
on the piecewise time allocation. We assume the robot body’s trapezoidal
velocity to simplify the time allocation process. Every piece accelerates to
max velocity, keeps its velocity, and decreases to $0m/s^{2}$.

###  III.B Virtual Model Trajectory Tracking Controller

The last subsection generates the robot’s base’s reference position, velocity,
and acceleration. Given the reference acceleration and reference angular
acceleration, the GRFs can be solved rapidly with a virtual model controller.
The virtual components producing virtual forces do not exist in the real
world. Choose spring and damper as the robot’s virtual components, and the
demonstration is shown as Fig. 4 and the relationship between
$\ddot{\bm{P}}_{com}$, $\dot{\bm{\omega}}_{B}$ and the reference trajectory
from minimum jerk is:

$\displaystyle\ddot{\bm{P}}_{com}^{r}$
$\displaystyle=\bm{k}_{p}^{p}(\bm{P}_{com}^{r}-\bm{P}_{com})+\bm{k}_{d}^{p}(\dot{\bm{P}}_{com}^{r}-\dot{\bm{P}}_{com})$
(7a) $\displaystyle\dot{\bm{\omega}}^{r}_{B}$
$\displaystyle=\bm{k}_{p}^{\omega}\bm{e}_{R}+\bm{k}_{d}^{\omega}(\bm{\omega}^{r}_{B}-\bm{\omega}_{B})$
(7b)

where $\ddot{\bm{P}}_{com}^{r}$ is the reference acceleration, and
$\dot{\bm{P}}_{com}^{r}$, $\bm{P}_{com}^{r}$ are reference base’s velocity and
position separately from the trajectory planner.
$\bm{k}_{p}^{p}\in\mathbb{R}^{3\times 3}$ and
$\bm{k}_{d}^{p}\in\mathbb{R}^{3\times 3}$ are the proportional and derivative
gains matrices. $\bm{\omega}_{B}^{r}$ is a constant, and $\bm{e}_{R}$ is the
error function for the rotation matrix, given by [19] as
$e_{R}=log(\bm{R}_{r}^{T}\cdot\bm{R})^{v}$. $\bm{k}_{p}^{\omega}$ and
${\bm{k}_{d}^{\omega}}$ are proportional and derivative gain matrices. The
reason for not considering the acceleration generated from the trajectory
planner is to get a better trajectory tracking performance while robots
prepare to jump. The accelerations produced by the virtual model promote the
robot’s dynamic efficiency and tracking performance with rarely any additional
computational cost.

Rewrite (2a) and (2b) in matrix formulation[20]:

$\underbrace{\left[\begin{array}[]{ccc}\mathbb{I}&\cdots&\mathbb{I}\\\
\bm{r}_{1}\times&\cdots&\bm{r}_{4}\times\end{array}\right]}_{\bm{M}\in\mathbb{R}^{6\times
12}}\underbrace{\left[\begin{array}[]{c}\bm{f}_{1}\\\ \vdots\\\
\bm{f}_{4}\end{array}\right]}_{\bm{f}\in\mathbb{R}^{12\times
1}}=\underbrace{\left[\begin{array}[]{c}m\ddot{\bm{P}}_{com}+\bm{g}\\\
\bm{I}_{B}\
\bm{\dot{}}{\omega}_{B}\end{array}\right]}_{\bm{N}\in\mathbb{R}^{6\times 1}}$
(8)

With known Cartesian and angular accelerations of the base, the constraint of
(8) is 6-dimension, but $f$ is 12-dimension. Thus, an optimization program can
be formed as follows:

$\displaystyle\bm{f}=\mathop{\arg\min}\limits_{f}(\bm{M}\bm{f}-$
$\displaystyle\bm{N})^{T}\bm{Q}(\bm{M}\bm{f}-\bm{N})+(\bm{J}_{c}\bm{f})^{T}\bm{R}\bm{J}_{c}\bm{f}$
(9a) $\displaystyle s.t.\quad\quad$
$\displaystyle\bm{c}_{l}<\bm{Gf}<\bm{c}_{u}$ (9b)

Where (9b) is the friction cone constraint of contact feet, whose matrices are
shown as:

$\displaystyle\bm{G}_{i}=\left[\begin{array}[]{c}-(\mu_{i}\bm{n}_{i})^{T}\\\
(\mu_{i}\bm{n}_{i})^{T}\\\
\bm{n}_{i}^{T}\end{array}\right],\bm{c}_{l}=\left[\begin{array}[]{c}-\infty\\\
0\\\
\bm{f}_{min_{i}}\end{array}\right],\bm{c}_{u}=\left[\begin{array}[]{c}0\\\
\infty\\\ \bm{f}_{max_{i}}\end{array}\right]$ (19)

$\mu_{i}\in\mathbb{R}$ is the friction parameter between the contact point and
the ground. $\bm{n}_{i}\in\mathbb{R}^{3}$ is the directional normal vector to
the ground. $\bm{Q}\in\mathbb{R}^{6\times 6}$ and
$\bm{R}\in\mathbb{R}^{6\times 6}$ is weight matrices and the term
$(\bm{J}_{c}\bm{f}^{T})\bm{R}\bm{J}_{c}\bm{f}$ penalizes large motor torque
output when robot’s base tracking the references.
$\bm{J}_{c}\in\mathbb{R}^{6\times 12}$ is the contact Jacobian Matrix.
$\bm{c}_{l}\in\mathbb{R}^{20\times 1}$ is the lower bound of friction cone
constrain and $\bm{c}_{u}\in\mathbb{R}^{20\times 1}$ is the upper bound.
$\bm{G}\in\mathbb{R}^{20\times 12}$ is the mapping matrix.

The 9a is a quadratic optimized problem and can be solved efficiently by
QuadProg++[21, 22]. The motor torque command is calculated by:

$\displaystyle\bm{\tau}_{ff}$ $\displaystyle=-\bm{S}_{f}\bm{J}_{c}\bm{f}$ (20)
$\displaystyle\bm{\tau}_{cmd}$
$\displaystyle=\bm{\tau_{ff}}+\bm{k}_{cp}(\bm{p}^{r}_{f}-\bm{p}_{f})+\bm{k}_{cd}({\bm{\dot{p}}^{r}_{f}-\bm{\dot{}}{\bm{p}}_{f}})$

where $\tau_{ff}\in\mathbb{R}^{12}$ is the motor feed-forward torque vector
and $\bm{S}_{f}\in\mathbb{R}^{12\times 12}$ is the selecting matrix.
$\bm{\tau}_{cmd}\in\mathbb{R}^{12}$ is the motor command torque.
$\bm{k}_{cp}\in\mathbb{R}^{12\times 12}$ and
$\bm{k}_{cd}\in\mathbb{R}^{12\times 12}$ are the Cartesian gain matrices
separately. $\bm{p}_{f}^{r}$ and $\bm{\dot{p}}_{f}^{r}$ are the reference foot
position and velocity vectors. $\bm{p}_{f}$ and $\bm{\dot{p}}_{f}$ are the
actual foot position and velocity vectors.

###  III.C WBC Landing Controller

As a quadratic program[23], the whole body controller with a controlling
frequency of 500Hz ensures the robot’s base stability in the landing phase.
When the robot touches the ground, the summation of reference GRFs equals its
gravity, and the reference base angles are all zeros. In this way, the robot
can land steadily and safely. The whole framework algorithm is shown as
algorithm 1. The motor torque commands from the whole-body controller come
from as follows:

$\bm{\tau}_{cmd}=\bm{\tau_{ff}}+\bm{k}_{p}(\bm{q}^{r}-\bm{q})+\bm{k}_{d}({\bm{\dot{q}}^{r}-\bm{\dot{}}{\bm{q}}})$
(21)

where $\bm{k}_{p}$ and $\bm{k}_{d}$ are the PD gains matrices. $\bm{q}^{r}$
and $\bm{\dot{q}}^{r}$ are the reference motor position and velocity vectors.
$\bm{q}$ and $\bm{\dot{q}}$ are the actual motor position and velocity
vectors.

input :
$\bm{P}_{\textrm{com}}^{\textrm{e}},\bm{\dot{P}}_{\textrm{com}}^{\textrm{e}},\bm{\ddot{P}}_{\textrm{com}}^{\textrm{e}},\bm{\omega}^{\textrm{e}},\psi^{\textrm{e}},\phi^{\textrm{e}},k$

output : ${\bm{\tau}_{\textrm{cmd}}}$

1 while  _$j <k$_ do

2 $s_{j}(t)=\sum_{i=0}^{5}a_{j,i}t^{i}\quad T_{j}\leq t\leq T_{j+1}$;

3 $s_{j}(t)^{{}^{\prime\prime\prime}}=\sum_{i\geq
3}^{5}i(i-1)(i-2)t^{i-3}a_{j,i}$;

4

5 end while

6$\bm{a}_{p}^{*}=-\bm{R}_{pp}^{-1}\bm{R}_{fp}^{T}\bm{a}_{f}$;

7

8$\ddot{\bm{P}}_{com}^{r}=\bm{k}_{p}^{p}(\bm{P}_{com}^{r}-\bm{P}_{com})+\bm{k}_{d}^{p}(\dot{\bm{P}}_{com}^{r}-\dot{\bm{P}}_{com})$;

9

10$\dot{\bm{\omega}}^{r}_{B}=\bm{k}_{p}^{\omega}e_{R}+\bm{k}_{d}^{\omega}(\bm{\omega}^{r}_{B}-\bm{\omega}_{B})$;

11

12while  _$n <maximum\ iterations$_ do

13 $\min_{f}\
(\bm{M}\bm{f}-\bm{N})^{T}\bm{Q}(\bm{M}\bm{f}-\bm{N})+(\bm{J}_{c}\bm{f})^{T}\bm{R}\bm{J}_{c}\bm{f}$

14 if  _$\bm{f}$ $<$ Tolerance_ then

15 Terminate

16 end if

17

18 end while

$\bm{\tau}_{cmd}=-\bm{S}_{f}\bm{J}_{c}\bm{f}$

Algorithm 1 The Framework Algorithm

##  IV Implementation Details and Experiments

This part focuses on the hardware implementation specifics of the online
omnidirectional jumping framework and experiments conducted using the improved
open-source Mini Cheetah[24]. Specifying hardware setup, we first offer
solving time data of a particular omnidirectional leaping framework to
illustrate the framework’s capacity to function in real time on the robot. We
showcase the framework’s capacity to facilitate Mini-Cheetah executing
omnidirectional jumps, such as front-left and rear-right jumps. Finally,
comparing the robot performance between the minimum jerk trajectory planner
and minimum snap trajectory planner illustrates the disparities between them.

###  IV.A Hardware Setup and Real-Time Performance

The framework developed with C++ is deployed in an enhanced Mini-Cheetah
equipped with an Intel NUC<EMAIL_ADDRESS>Hz. The trajectory planning, VMC, and
whole-body controllers run on the NUC in real-time. Motor data is collected at
1000 Hz through a New-design USB to Can board. The Lightweight Communications
and Marshalling package (LCM) is used for the asynchronous communication
between the simulation, hardware bridge, and data collecting.

Hyper-parameters are shown in Tab. I, including hyperparameters of mini jerk
trajectory planner and hyperparameters in VMC.

To verify the real-time performance of the fast omnidirectional jumping
framework, with the random input reference end CoM position, velocity, and
acceleration w.r.t the inertial frame, we record the time-consuming summation
of the minimum jerk trajectory planner and VMC controller, which starts from
the reference inputting and comes to an end with the generation of motor
torque commands. Several time measurements are repeated on actual robots, and
the average consuming time is compared with our previous work, the
evolutionary-based jumping framework[10] and the offline jumping framework[9].

TABLE I: HYPERPARAMETERS

Parameters | Symbol | Values
---|---|---
Minimum Jerk
End Position Range | $\bm{P}_{CoM}[m]$ | $[\pm 0.15,\pm 0.1,\pm 0.05]$
End Velocity Range | $\bm{\dot{}}{P}_{CoM}[m/s]$ | $[\pm 0.3,\pm 0.16,[1.5,3.5]]$
End Acceleration Range | $\bm{\ddot{P}}_{CoM}^{z}[m/s^{2}]$ | $[\pm 20,\pm 20,[10,40]]$
VMC
Friction Coefficient | $\mu$ | 0.5
Weights for CoM Wrench | $\bm{Q}$ | $diag(1,1,10,20,10,25)$
Weights for Torque Wrench | $\bm{R}$ | $diag(5,50,2)10^{-5}$
CoM Proportional Gains | $\bm{k}_{p}^{p}$ | $diag(1070,1070,1070)$
CoM Derivative Gains | $\bm{k}_{d}^{p}$ | $diag(12,12,10)$
Attitude Proportional Gains | $\bm{k}_{p}^{\omega}$ | $diag(800,800,800)$
Attitude Derivative Gains | $\bm{k}_{d}^{\omega}$ | $diag(20,10,20)$
GRFs Minimum | $\bm{f}_{min}[N]$ | 5
GRFs Maximum | $\bm{f}_{max}[N]$ | 250
Cartesian Proportional Gains | $\bm{k}_{cp}$ | $diag(0,0,0)$
Cartesian Derivative Gains | $\bm{k}_{cd}$ | $diag(15,15,15)$

As Tab. II demonstrated, the average time resource required for trajectory
planner and VMC is around 50 microseconds on NUC, depicting well the real-time
efficacy of the framework.

TABLE II: JUMPING ALGORITHM COMPARISON Algorithm | Average Time(s)
---|---
This work | 0.00005
Evolution-Based Jumping Framework[10] | 0.14
Offline Jumping Framework[9] | 65

###  IV.B Ominidirectional Jumping

Made consisting of the jumping phase, flight phase, and landing phase, the
omnidirectional jumping framework is capable of enabling Mini-Cheetah to jump
in any direction w.r.t its base frame, such as the front side, the left side,
the rear side, and angles between them. Fig. 5 illustrates the robot’s ability
to perform omnidirectional jumping well.

Figure 5: Mini-Cheetah deployed with the omnidirectional jumping framework
performs two jumping motions. Data charts in Fig. (a) and Fig. (b) are data
records of the rear-right leg in front-left jumping and rear-right jumping,
respectively. The blue line represents the start of the jumping phase, the
dark red line indicates the start of the flight phase, the green line is the
beginning of the landing phase, and the red dot lines represent the upper and
lower bound of the torque limit with an absolute value of 24 Nm.

As the figures show, the period between the start of the jumping phase and the
beginning of the flight phase is within 0.5 seconds. The maximum of motor
torques in the jumping phase is around 22 Nm. Because of the imprecision of
the SRB model in dynamics, motor torques generated by the VMC are smaller than
the actual desired torques. Thus, a Cartesian PD controller is implemented to
improve the performance of foot position and speed tracking in the jumping
phase, conducting a better jumping motion. The front-left jumping direction
and the rear-right jumping direction are selected intentionally to show the
robot’s potential capability of performing omnidirectional jumping, which
matches the result of the experiments.

##  V CONCLUSIONS

In conclusion, this paper introduces a new omnidirectional jumping framework
for quadruped robots, leveraging Virtual Model Control (VMC) and Minimum Jerk
principles. This framework is distinguished by its rapid execution, achieving
microsecond levels with impressive speed. By employing a 5th-order polynomial
for trajectory generation alongside defined constraints for continuity, the
minimum-jerk trajectory method proves highly effective for this application.
Implementing a VMC-based trajectory tracking controller, utilizing PD control
for centroid trajectory tracking and GRFs optimization, has shown remarkable
efficiency, and in real-world jumping experiments, the solution time for
trajectory tracking reached approximately 30 us, with an additional 20 us for
trajectory generation. This significantly enhances the optimization time to 50
us, marking a substantial improvement over previous methods that relied on the
differential evolution algorithm. Such computing speed is critical for
fulfilling the real-time requirements of emergency obstacle avoidance tasks in
robotic applications.

Looking ahead, our future work will explore improved time allocation
techniques and enhanced precision in jump control, particularly in conjunction
with external localization systems. We also aim to address the challenges
observed with the current impulse-based WBC landing controller[5], which
sometimes obtained failure results in landing phases. A primary focus will be
developing a new landing controller that effectively manages horizontal and
vertical speeds to preserve landing posture.

## References

  * [1] M. Hutter, C. Gehring, A. Lauber, F. Gunther, C. D. Bellicoso, V. Tsounis, P. Fankhauser, R. Diethelm, S. Bachmann, M. Blösch, et al., “Anymal-toward legged robots for harsh environments,” Advanced Robotics, vol. 31, no. 17, pp. 918–931, 2017
  * [2] J. K. Yim, B. R. P. Singh, E. K. Wang, R. Featherstone and R. S. Fearing, ”Precision Robotic Leaping and Landing Using Stance-Phase Balance,” Robotics and Automation Letters, vol. 5, no. 2, pp. 3422-3429, 2020
  * [3] H.-W. Park, M. Y. Chuah, and S. Kim, “Quadruped bounding control with variable duty cycle via vertical impulse scaling,” in 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems, 2014, pp. 3245–3252.
  * [4] J. Di Carlo, P. M. Wensing, B. Katz, G. Bledt and S. Kim, ”Dynamic Locomotion in the MIT Cheetah 3 Through Convex Model-Predictive Control,” 2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2018, pp. 1-9
  * [5] Kim, D., Carlo, J.D., Katz, B., Bledt, G., Kim, S. (2019). ”Highly Dynamic Quadruped Locomotion via Whole-Body Impulse Control and Model Predictive Control”. ArXiv, abs/1909.06586.
  * [6] Xie H, Ahmadi M, Shang J, Luo Z. An intuitive approach for quadruped robot trotting based on virtual model control. Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering. 2015;229(4):342-355
  * [7] Park H W, Wensing P M, Kim S. Jumping over obstacles with MIT Cheetah 2[J]. Robotics and Autonomous Systems, 2021, 136: 103703.
  * [8] M. Chignoli and S. Kim, ”Online Trajectory Optimization for Dynamic Aerial Motions of a Quadruped Robot,” 2021 IEEE International Conference on Robotics and Automation (ICRA), 2021, pp. 7693-7699
  * [9] Z. Song et al., ”An Optimal Motion Planning Framework for Quadruped Jumping,” 2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2022, pp. 11366-11373
  * [10] L. Yue, Z. Song, H. Zhang, X. Zeng, L. Zhang, and Y. -H. Liu, ”Evolutionary-Based Online Motion Planning Framework for Quadruped Robot Jumping,” 2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2023, pp. 767-773
  * [11] Bellegarda G, Nguyen C, Nguyen Q. Robust quadruped jumping via deep reinforcement learning[J]. arXiv preprint arXiv:2011.07089, 2020.
  * [12] “Unitree a1 robot,” https://www.unitree.com/products/ a1/.
  * [13] Rudin N, Kolvenbach H, Tsounis V, et al. Cat-like jumping and landing of legged robots in low gravity using deep reinforcement learning[J]. IEEE Transactions on Robotics, 2021, 38(1): 317-328.
  * [14] L. Tang, Y. Li and B. Li, ”Moobot: A Miniature Origami Omnidirectional Jumping Robot With High Trajectory Accuracy,” IEEE Transactions on Industrial Electronics, vol. 71, no. 6, pp. 6032-6040
  * [15] Chignoli M, Morozov S, Kim S. ”Rapid and reliable quadruped motion planning with omnidirectional jumping,” 2022 International Conference on Robotics and Automation (ICRA). 2022: 6621-6627.
  * [16] Mellinger D, Kumar V. Minimum snap trajectory generation and control for quadrotors[C]//2011 IEEE international conference on robotics and automation. IEEE, 2011: 2520-2525.
  * [17] Mu B, Chirarattananon P. Trajectory generation for underactuated multirotor vehicles with tilted propellers via a flatness-based method[C]//2019 IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM). IEEE, 2019: 1365-1370.
  * [18] Richter C, Bry A, Roy N. Polynomial trajectory planning for aggressive quadrotor flight in dense indoor environments[C]//Robotics Research: The 16th International Symposium ISRR. Springer International Publishing, 2016: 649-666.
  * [19] F. Bullo and A. D. Lewis, Geometric Control of Mechanical Systems: Modeling, Analysis, and Design for Simple Mechanical Control Systems, vol. 49. Cham, Switzerland: Springer Science & Business Media, 2004.
  * [20] Focchi M, Del Prete A, Havoutis I, et al. High-slope terrain locomotion for torque-controlled quadruped robots[J]. Autonomous Robots, 2017, 41: 259-272.
  * [21] Goldfarb D, Idnani A. A numerically stable dual method for solving strictly convex quadratic programs[J]. Mathematical programming, 1983, 27(1): 1-33.
  * [22] https://github.com/liuq/QuadProgpp
  * [23] S. Fahmi, C. Mastalli, M. Focchi and C. Semini, ”Passive Whole-Body Control for Quadruped Robots: Experimental Validation Over Challenging Terrain,” in IEEE Robotics and Automation Letters, vol. 4, no. 3, pp. 2553-2560
  * [24] Katz B, Di Carlo J, Kim S. Mini cheetah: A platform for pushing the limits of dynamic quadruped control[C]. 2019 International Conference on Robotics and Automation (ICRA). , IEEE, 2019: 6295-6301.


# FSSH-2: Fewest Switches Surface Hopping with robust switching probability

Leonardo Araujo<EMAIL_ADDRESS>Caroline Lasser<EMAIL_ADDRESS>Department of Mathematics, TUM School of Computation, Information and
Technology, Technische Universität München, Boltzmannstr. 3, 85748 Garching
bei München  Burkhard Schmidt<EMAIL_ADDRESS>Weierstrass
Institute for Applied Analysis and Stochastics, Mohrenstr. 39, 10117 Berlin

###### Abstract

This study introduces the FSSH-2 scheme, a redefined and numerically stable
adiabatic Fewest Switches Surface Hopping (FSSH) method for mixed quantum-
classical dynamics. It reformulates the standard FSSH hopping probability
without non-adiabatic coupling vectors and allows for numerical time
integration with larger step sizes. The advantages of FSSH-2 are demonstrated
by numerical experiments for five different model systems in one and two
spatial dimensions with up to three electronic states.

…

##  I Introduction

The widely–used Fewest Switches Surface Hopping (FSSH) method is a
cost–effective approach for the numerical simulation of non–adiabatic effects
in quantum–classical dynamics. Even though originally developed more than 30
years ago [1], this concept is still widely used, especially in the context of
photochemistry and excited-state molecular dynamics [2, 3, 4, 5]. The FSSH
scheme adopts a mixed quantum–classical paradigm where typically nuclei are
treated classically and electrons quantum-mechanically. The method’s principal
efficiency stems from its use of classical trajectories moving on adiabatic
potential energy surfaces with the possibility of stochastic ’hopping’ between
states according to calculated probabilities. Because of this conceptual
simplicity and the straightforward implementation, FSSH is still among the
most commonly used methods for non–adiabatic molecular dynamics.

Despite its efficiency, FSSH has the drawback that it relies on the notorious
Non–Adiabatic Couplings (NACs). While there are methods to compute NACs
accurately and efficiently in regions of large gaps between neighboring
potential energy surfaces [6, 7, 4], they may become increasingly large near
energy crossings, if not divergent. For trajectories approaching genuine or
avoided crossings, this necessitates the use of exceptionally small time step
sizes to achieve numerical convergence, if at all possible. The use of small
time steps may result in crucial problems when combined with "on the fly"
electronic structure calculations of complex chemical systems. There, even a
single time step may become computationally expensive, in particular because
of the calculations of forces and NAC vectors. This underscores the necessity
of robust time integration approaches that allow for large time steps.

The present study introduces FSSH-2, an improved numerical time integration
scheme for the FSSH method. FSSH-2 incorporates a modified fewest switches
hopping probability that utilizes time–overlap integrals for electronic
amplitude propagation, a technique known as local diabatization [8]. An
essential feature of this approach is its numerical stability, as it operates
independently of the potentially divergent NACs. It is similar to the global
flux hopping method [9], but stands apart in both its derivation and its
efficient formulation.

We begin with a brief derivation of the conventional FSSH method in Sec. II,
highlighting its assumptions, numerical aspects, and established challenges.
In Sec. III, the FSSH-2 method is presented, along with a discussion of its
theoretical advantages over FSSH. In Sec. IV, we present numerical experiments
demonstrating the advantages of FSSH-2.

##  II FSSH Method

The traditional FSSH method propagates swarms of classical trajectories, that
represent the nuclei, with associated electronic amplitudes that trigger
hopping events to model non–adiabatic transitions. We employ a fixed time step
size denoted by $\Delta_{t}>0$, with time steps represented as
$t_{s}=s\Delta_{t}$ for $s\in{\mathbb{N}}_{0}$. We discuss the three
fundamental tasks for advancing a given classical trajectory with position
$q(t)\in{\mathbb{R}}^{d_{\text{nuc}}}$ and momentum
$p(t)\in{\mathbb{R}}^{d_{\text{nuc}}}$ one time step further: classical
propagation, evolution of the corresponding electronic amplitudes, and surface
hopping.

###  II.1 Classical Dynamics

For each state, let the function $\theta_{m}(q,x)$ be the $m$-th adiabatic
state, where $q\in{\mathbb{R}}^{d_{\text{nuc}}}$ represents the nuclear
coordinates and $x\in{\mathbb{R}}^{d_{\text{el}}}$ the electronic coordinates.
Each adiabatic state is an eigenfunction corresponding to the $m$-th
eigenvalue of the self–adjoint electronic Hamiltonian $H^{\text{el}}(q)$
parametrically dependent on the coordinates $q$.. We assume that the adiabatic
states are real–valued. The eigenvalue equation for the electronic Hamiltonian
can be written as

$\displaystyle H^{\text{el}}(q)\theta_{m}(q,x)=E_{m}(q)\theta_{m}(q,x),$

where $E_{m}(q)$ defines the $m$-th adiabatic potential energy
(hyper-)surface. The classical dynamics of a trajectory occupying the
adiabatic state $\theta_{m}$ is described by Hamilton’s equations of motion:

$\displaystyle\begin{split}\dot{q}(t)&=\frac{p(t)}{M}\\\
\dot{p}(t)&=F_{m}(q(t)).\end{split}$ (1)

Here, $M>0$ stands for the nuclear mass and
$F_{m}(q)=-\nabla_{q}E_{m}(q)\in{\mathbb{R}}^{d_{\text{nuc}}}$ represents the
$m$-th adiabatic force which can be computed by means of the Hellmann–Feynman
theorem [10, 11, 12],

$\displaystyle
F_{m}(q)=-\braket{\theta_{m}(q,\cdot),\nabla_{q}H^{\text{el}}(q)\theta_{m}(q,\cdot)},$
(2)

where $\braket{\cdot,\cdot}$ represents the inner product of the electronic
Hilbert space.

For propagating Eq. (1), we used the Störmer–Verlet method [13] since it is an
explicit, second–order, symplectic integrator. Due to the singularity of the
adiabatic force near eigenvalue crossings, we refrain from using a higher
order method.

###  II.2 Electronic Amplitudes

Let us denote by $L\in{\mathbb{N}}$ the total number of electronic states
being considered. To compute hopping probabilities, FSSH propagates adiabatic
amplitudes, denoted as
$c^{\text{adi}}(t)=\left(c^{\text{adi}}_{m}(t)\right)_{m=1}^{L}\in{\mathbb{C}}^{L}$.
These amplitudes obey the following ordinary differential equation:

$\displaystyle i\dot{c}^{\text{adi}}(t)=H^{\text{adi}}(t)c^{\text{adi}}(t).$
(3)

In this equation, the elements of the Hamiltonian matrix
$H^{\text{adi}}(t)=\left(H^{\text{adi}}_{m,n}(t)\right)_{m,n=1}^{L}\in{\mathbb{C}}^{L\times
L}$ are defined as

$\displaystyle
H^{\text{adi}}_{m,n}(t)=\delta_{m,n}E_{m}(q(t))-\frac{i}{M}d^{m,n}(q(t))^{T}p(t),$

where $d^{m,n}(q)=\braket{\theta_{m}(q,\cdot),\nabla_{q}\theta_{n}(q,\cdot)}$
represents the non–adiabatic coupling (NAC) vector. Initially normalized
amplitudes stay normalized for all times,$\|c^{\text{adi}}(t)\|=1$. Similarly
to the adiabatic force in Eq. (2), the (extended) Hellmann–Feynman theorem
provides a straightforward formula for the computation of the NAC vector for
$E^{m}(q)\neq E^{n}(q)$:

$\displaystyle
d^{m,n}(q)=-\frac{\braket{\theta_{m}(q,\cdot),\nabla_{q}H^{\text{el}}(q)\theta_{n}(q,\cdot)}}{E_{m}(q)-E_{n}(q)}.$
(4)

However, propagating the adiabatic electronic amplitudes is significantly
hampered by the divergence of the NAC vector, and the literature [14, 15]
recommends to introduce a second, smaller time step size for this task.

The direct propagation of Eq. (3) is typically carried out using the first
term of the Magnus expansion [16] [17, Chapter IV.7] in conjunction with the
trapezoidal rule:

$\displaystyle
c^{\text{adi}}(t_{s+1})\approx\exp\left({-i\int\limits_{t_{s}}^{t_{s+1}}H^{\text{adi}}(t)dt}\right)c^{\text{adi}}(t_{s})\approx\exp\left({-i\Delta_{t}\frac{H^{\text{adi}}(t_{s})+H^{\text{adi}}(t_{s+1})}{2}}\right)c^{\text{adi}}(t_{s}).$
(5)

The utilization of the trapezoidal rule here is advantageous as it circumvents
the need to calculate adiabatic states at locations not associated with the
defined time steps. As for the classical dynamics, a higher order
approximation could lead to faster divergences for larger time step sizes.

###  II.3 Surface Hopping

After propagating the electronic amplitudes, a trajectory may undergo a
transition from its current state $m$ to a new state $n$. The FSSH hopping
probability $P_{m,n}(t_{s},t_{s+1})$ is derived by considering the electronic
densities
$\rho(t)=c^{\text{adi}}(t)\overline{c^{\text{adi}}(t)}^{T}\in\mathbb{C}^{L\times
L}$ based on the electronic amplitudes. By means of Eq. (3), the rate of
change of these densities are characterized by

$\displaystyle\dot{\rho}_{m,m}(t)$
$\displaystyle=\frac{2}{M}\sum\limits_{\begin{subarray}{c}n=1\\\ n\neq
m\end{subarray}}^{L}\text{Re}\left(\rho_{m,n}(t)\right)\
d^{m,n}(q(t))^{T}p(t).$ (6)

The FSSH hopping probability is obtained by tracking the outflux from state
$m$ to all other states at time $t$ during a short time interval $\Delta_{t}$.
Assuming no influx into state $m$ during this brief period, the rate of
outflux is described as $\rho_{m,m}(t_{s})-\rho_{m,m}(t_{s+1})$. The
probabilities of staying at and leaving state $m$ are then given by:

$\displaystyle
P_{m,\text{stay}}(t_{s},t_{s+1})=\frac{\rho_{m,m}(t_{s+1})}{\rho_{m,m}(t_{s})}\quad\text{
and }\quad
P_{m,\text{out}}(t_{s},t_{s+1})=\frac{\rho_{m,m}(t_{s})-\rho_{m,m}(t_{s+1})}{\rho_{m,m}(t_{s})}.$
(7)

This is known as the fewest switches assumption after which the FSSH method is
named. To compute the probability of hopping from state $m$ to another state
$n\neq m$, we can express $P_{m,n}$ by inserting Eq. (6) into (7):

$\displaystyle
P_{m,\text{out}}(t_{s},t_{s+1})=-\frac{\int_{t_{s}}^{t_{s+1}}\dot{\rho}_{m,m}(t)dt}{\rho_{m,m}(t_{s})}=-\frac{2}{M\rho_{m,m}(t_{s})}\sum\limits_{\begin{subarray}{c}n=1\\\
n\neq
m\end{subarray}}^{L}\int_{t_{s}}^{t_{s+1}}\text{Re}\left(\rho_{m,n}(t)\right)\
d^{m,n}(q(t))^{T}p(t)dt.$

The last equation formulates a hopping probability from state $m$ to $n$ as:

$\displaystyle P_{m,n}(t_{s},t_{s+1})$
$\displaystyle=-\frac{2}{M\rho_{m,m}(t_{s})}\int_{t_{s}}^{t_{s+1}}\text{Re}\left(\rho_{m,n}(t)\right)\
d^{m,n}(q(t))^{T}p(t)dt.$ (8)

By applying an Euler approximation to the integral, the original FSSH formula
[1] is reconstructed

$\displaystyle P_{m,n}(t_{s},t_{s+1})$
$\displaystyle=-\frac{2\Delta_{t}\text{Re}\left(\rho_{m,n}(t_{s+1})\right)}{M\rho_{m,m}(t_{s})}d^{m,n}(q(t_{s+1}))^{T}p(t_{s+1}),$
(9)

where negative values are treated as zeros. Note that $P_{m,n}$ depends on the
NAC, thus inheriting its possible divergence issue.

For a hopping event from state $m$ to another state $n$, the process involves
a uniformly distributed random number $r\in[0,1]$. The state $n$ is determined
by the condition:

$\displaystyle\sum_{\begin{subarray}{c}k=1\\\ k\neq
m\end{subarray}}^{n-1}P_{m,k}(t_{s},t_{s+1})\leq
r<\sum_{\begin{subarray}{c}k=1\\\ k\neq
m\end{subarray}}^{n}P_{m,k}(t_{s},t_{s+1}).$

If no such $n$ satisfies this criterion, no hopping is performed and the
trajectory remains in state $m$. This random hopping criterion together with
initial random sampling constitute the stochastic nature of the FSSH method.

##  III FSSH-2 Method

The above numerical issues have led us to explore alternative and more
efficient approaches to refine the FSSH method. The novel FSSH-2 presents two
modifications to the original FSSH method: the propagation of adiabatic
amplitudes through local diabatization and the redefinition of FSSH hopping
probabilities.

###  III.1 Local Diabatization Method for Propagating Adiabatic Amplitudes

In the spirit of the local diabatization method [8], we use a diabatic
representation, which often offers smoother conditions. For the while being,
we assume that a diabatic potential energy matrix
$V^{\text{dia}}(q)\in{\mathbb{R}}^{L\times L}$ is given but later on it will
turn out that an adiabatic–to–diabatic transformation is not required. Instead
of considering the adiabatic amplitudes, we work in a diabatic representation
with diabatic amplitudes that evolve according to

$\displaystyle i\dot{c}^{\text{dia}}(t)=H^{\text{dia}}(t)c^{\text{dia}}(t),$
(10)

where the diabatic Hamiltonian consists solely of the diabatic potential
energy matrix, $H^{\text{dia}}(t)=V^{\text{dia}}(q(t))$. The diabatic
amplitudes are related to the adiabatic amplitudes through the diagonalization
of $V^{\text{dia}}(q)$

$\displaystyle V^{\text{dia}}(q)=\overline{U(q)}^{T}E(q)U(q),$ (11)

where
$E(q)=\operatorname{\text{diag}}\left(E_{1}(q),...,E_{L}(q)\right)\in{\mathbb{R}}^{L\times
L}$ denotes the adiabatic potential energy matrix, which is diagonal.
Consequently, the unitary matrix $U(q(t))\in{\mathbb{R}}^{L\times L}$ serves
as the transformation matrix from diabatic to adiabatic representation. This
implies that the electronic amplitudes in the diabatic and adiabatic
representations are related as follows:

$\displaystyle c^{\text{adi}}(t)=U(q(t))c^{\text{dia}}(t).$ (12)

Note that we will show below that the transformation matrices $U$ will not be
used in the FSSH-2 scheme. Analogously to Eq. (5), we can propagate Eq. (10)
using the first term of the Magnus expansion but this time in conjunction with
the Euler method

$\displaystyle
c^{\text{dia}}(t_{s+1})\approx\exp\left({-i\int\limits_{t_{s}}^{t_{s+1}}V^{\text{dia}}(q(t))dt}\right)c^{\text{dia}}(t_{s})\approx\exp\left({-i\
V^{\text{dia}}(q(t_{s+1}))\ \Delta_{t}}\right)c^{\text{dia}}(t_{s}).$

The Euler method is beneficial here since the diagonalization of
$V^{\text{dia}}(q(t))$ is inherently involved in the force calculation for
classical dynamics when using a diabatic model. Together with Eq. (11), this
leads to:

$\displaystyle c^{\text{dia}}(t_{s+1})$ $\displaystyle\approx\exp\left({-i\
V^{\text{dia}}(q(t_{s+1}))\ \Delta_{t}}\right)c^{\text{dia}}(t_{s})$
$\displaystyle=\overline{U(q(t_{s+1}))}^{T}\exp\left({-i\ E(q(t_{s+1}))\
\Delta_{t}}\right)U(q(t_{s+1}))\ c^{\text{dia}}(t_{s}).$

Transforming the diabatic amplitudes back into the adiabatic representation by
using Eq. (12), we obtain the following propagator for the adiabatic
amplitudes:

$\displaystyle c^{\text{adi}}(t_{s+1})$ $\displaystyle\approx\exp\left({-i\
E(q(t_{s+1}))\ \Delta_{t}}\right)\Gamma(t_{s},t_{s+1})\ c^{\text{adi}}(t_{s})$
(13)

with
$\Gamma(t_{s},t_{s+1})=U(q(t_{s+1}))\overline{U(q(t_{s}))}^{T}\in{\mathbb{R}}^{L\times
L}$. Evaluating the remaining matrix exponential is straightforward, as $E(q)$
is a diagonal matrix. Interestingly, this propagator remains valid even in the
absence of a known diabatic framework. The matrix elements of
$\Gamma(t_{s},t_{s+1})$ are defined as

$\displaystyle\Gamma_{m,n}(t_{s},t_{s+1})=\braket{\theta_{m}(q(t_{s+1}),\cdot),\theta_{n}(q(t_{s}),\cdot)},$

where the $\theta$ are the adiabatic states defined in Sec. II.1. In Ref.
[18], these matrix elements are referred to as time–overlaps between adiabatic
states. From the point of view of numerical calculations, they have the
advantage that their absolute value is bounded by one. Thus, Eq. (13) results
in a propagator of the adiabatic states that only depends on the adiabatic
potential energies, $E$, and time–overlaps, $\Gamma$, of adiabatic states.
Hence, this propagator offers several advantages in addressing the numerical
challenges discussed for Eq. (5). First, it eliminates the need to consider
the NAC vectors, thus ensuring there are no divergent quantities in the
calculations and computational costs do not scale with $\mathcal{O}(L^{4})$,
as demonstrated in our previous work [19]. In fact, the computational costs of
$\Gamma$ are comparable to those of matrix multiplication, scaling with
$\mathcal{O}(L^{2})$. Second, the implementation of FSSH-2 is straightforward
and all operations in Eq. (13) can be trivially parallelized across all
trajectories. Thus, the local diabatization propagator offers a simple,
cost–effective, and parallelizable method for propagating adiabatic
amplitudes.

###  III.2 FSSH-2 Hopping Formula

Another challenge is to derive a hopping formula that can harness the
advantages of the local diabatization propagator fully. The traditional FSSH
hopping formula exhibits the stability issues mentioned in Sec. II.2 and II.3.
To tackle this, we revisit the outflow probability in Eq. (7). By first
considering the case of only two electronic states ($L=2$), the outflow
probability characterizes the hopping probability from state $m$ to the other
state $n\neq m$. Utilizing the normalization of our electronic densities,
$\rho_{m,m}(t)+\rho_{n,n}(t)=1$, we can express the outflow probability in
terms of state $n$:

$\displaystyle
P_{m,\text{out}}(t_{s},t_{s+1})=\frac{\rho_{m,m}(t_{s})-\rho_{m,m}(t_{s+1})}{\rho_{m,m}(t_{s})}=\frac{\rho_{n,n}(t_{s+1})-\rho_{n,n}(t_{s})}{\rho_{m,m}(t_{s})}.$
(14)

This hopping formula is straightforward to evaluate, given that all its
elements are densities, bounded by 1. Notably, this formula is not only more
concise than the FSSH hopping formula in Eq. (9) but it is also independent of
NACs. Consequently, it provides a simple approach to circumvent potential
challenges of NACs.

We extend this approach to an arbitrary number of electronic states
$L\in\mathbb{N}$. The normalization condition for electron densities can be
expressed as

$\displaystyle\sum_{n=1}^{L}\rho_{n,n}(t)=1\quad\Leftrightarrow\quad\rho_{m,m}(t)=1-\sum_{\begin{subarray}{c}n=1\\\
n\neq m\end{subarray}}^{L}\rho_{n,n}(t).$

Substituting this into Eq. (7), we obtain the expression

$\displaystyle P_{m,\text{out}}(t_{s},t_{s+1})=\sum_{\begin{subarray}{c}n=1\\\
n\neq
m\end{subarray}}^{L}\frac{\rho_{n,n}(t_{s+1})-\rho_{n,n}(t_{s})}{\rho_{m,m}(t_{s})}.$

The individual terms in this sum,
$\frac{\rho_{n,n}(t_{s+1})-\rho_{n,n}(t_{s})}{\rho_{m,m}(t_{s})}$, indicate a
hopping probability between state $m$ and $n$, analogous to Eq. (14) for the
case of two electronic states. However, for $L>2$, the physical interpretation
of these terms might be incorrect in certain scenarios. For instance, when
$\rho_{n,n}(t_{s+1})-\rho_{n,n}(t_{s})>0$ and the outflow of state $m$ is
zero, it could incorrectly suggest a hopping from $m$ to $n$. Therefore, we
introduce $P_{m,\text{out}}$ as an upper bound and define the FSSH-2 hopping
probability as

$\displaystyle\tilde{P}_{m,n}(t_{s},t_{s+1})=\min\left\\{P_{m,\text{out}}(t_{s},t_{s+1}),\frac{\rho_{n,n}(t_{s+1})-\rho_{n,n}(t_{s})}{\rho_{m,m}(t_{s})}\right\\}.$
(15)

As mentioned for the case $L=2$, it addresses some of the issues encountered
with the traditional FSSH hopping formula. Based solely on the densities
derived from adiabatic amplitudes, $\tilde{P}_{m,n}$ circumvents the issues
related to the divergent NAC vectors. Additionally, the computation of this
refined hopping formula is simple, enhancing its practicality. With the
absence of approximations in its formulation, the precision of
$\tilde{P}_{m,n}$ is solely connected to the accuracy of the propagation of
the adiabatic amplitudes.

###  III.3 Limitations of the FSSH-2 hopping probability

The FSSH-2 hopping probability still presents certain shortcomings. Like FSSH,
the probabilities do not form a valid probability distribution. Another
concern is that the upper bound in Eq. (15) may not be sufficient to prevent
incorrect transfers between states that are not coupled but have transfers at
the same time. Imagine a scenario with four states where the first and last
two states cross each other simultaneously, but there is a substantial gap
between states 2 and 3. If there is population on all states, it may occur
that a hopping event from state 3 to 2 occurs.

##  IV Simulations

Model |  Difference between quantum solution and FSSH |  Difference between quantum solution and FSSH-2 |  Difference between FSSH
and FSSH-2
---|---|---|---
Tully 1 |  0.004 |  0.003 |  0.001
Tully 2 |  0.008 |  0.004 |  0.007
2D LVC |  0.015 |  0.012 |  0.014
2D Well |  0.023 |  0.021 |  0.004
Model X |  0.007 |  0.010 |  0.003

Table 1:  Comparison of the quantum reference solution to converged results of
FSSH and FSSH-2 for each model, evaluated at our highest chosen (256) substep
number. The numbers showcase the mean absolute population differences of the
second adiabatic state.

The FSSH-2 method offers a straightforward yet promising redefinition of the
traditional FSSH method, showing potential for enhanced numerical efficiency,
which needs further validation across various scenarios. To assess the
performance of the two methods, we explored a variety of diabatic models of
non–adiabatic dynamics in molecular systems: the Tully single and double
crossing models [1], the two–dimensional linear vibronic coupling (2D LVC)
model [20], the two–dimensional well [21], and Model X [22]. While the latter
one incorporates three electronic states, the other models feature two
electronic states. This exhaustive approach is in line with the recommended
test suite for non–adiabatic algorithms by Nelson et al. [23].

Our simulations were carried out using version 7.2.0 of the WavePacket MATLAB
software package [24, 25, 19], which is designed to treat fully
quantum–mechanical (grid based) and quantum–classical (trajectory based)
simulations on an equal footing. Given that the use of grid–based
representations of wave functions is restricted to low dimensional systems,
its quantum–classical counterpart offers an efficient alternative for
simulations involving higher dimensional systems. Within WavePacket, the main
time steps can be divided into an arbitrary number of substeps, which are then
used as our previously defined (propagation) time steps $\Delta_{t}$.
Therefore, by keeping the main time steps fixed for a model, we can vary
$\Delta_{t}$ by varying the number of substeps. This also enables us to easily
maintain consistency across the models by selecting the same substep number.

In our simulations, we chose the simulation time such that the classical
trajectories could traverse the (avoided) crossings up to two times. The
surface hopping methods used swarms of 10,000 trajectories. Following each
hopping event, the momentum was rescaled along the NAC vectors, the direction
of which can be efficiently computed without singularity issues by evaluating
only the numerator in Eq. (4). The main time step sizes were selected based on
the limits of energy conservation of the classical dynamics propagator. To
compare our results to a full quantum solution, we also performed a wave
function simulation using a second–order Strang-Marchuk operator splitting
scheme with the FFT grid representation [26, 27, 28]. Further details for each
simulation set–up can be found in the appendix.

Figure 1:  Population of the second adiabatic state for the Tully 1 model
using FSSH (left) and FSSH-2 (center) for different numbers of substeps. The
right panel shows the absolute population error after the initial hopping
event for both methods.

Figure 2:  Same as Fig. 1 but for the Tully 2 model.

Figure 3:  Same as Fig. 1 but for the 2D LVC model.

Figure 4:  Same as Fig. 1 but for the 2D Well model.

Figure 5:  Same as Fig. 1 but for Model X.

###  IV.1  Comparison between FSSH and FSSH-2

The primary aspect of our investigation involves the comparison of population
transfers. In the case of two electronic states and under the assumption of
continuous adiabatic states, FSSH-2 and FSSH exhibit a first–order
approximation relationship. Consequently, for models featuring only two
electronic states, we anticipate close accordance between FSSH and FSSH-2
results. The sole scenario where qualitative discrepancies might emerge is
Model X because it involves three electronic states. Nonetheless, as shown in
Tab. 1, across all models considered, the absolute differences in population
transfers averaged over the running time generated by these methods with their
maximum substep number of 256 were consistently below $0.025$. This level of
variance is acceptable for a stochastic simulation method, considering the
number of trajectories involved. Except for the 2D Well case which displayed
the largest deviations, no discernible advantages were observed when these
results were compared against the quantum reference solutions.

###  IV.2 Toward larger time step sizes

Next, we assessed how FSSH and FSSH-2 responded to enlarged time step sizes by
systematically halving the substep number, starting from 256 and descending to
1. Obviously, each decrease in the substep number directly cuts the total
simulation time in half. Such a gain of efficiency has profound implications
for realistic system simulations, especially when surface hopping simulations
are coupled to large scale electronic structure calculations.

In Figs. 1 to 5, the left and center panels illustrate the population dynamics
of the second adiabatic state for both the traditional FSSH and our novel
FSSH-2 method, respectively. For the standard FSSH method, results with 1, 2,
and 4 substeps often significantly deviate from the converged outcomes,
sometimes even presenting qualitatively incorrect results. In contrast, the
results for our FSSH-2 method demonstrate a consistently much better
performance - in most cases, the results with just a single substep per main
time step align closely with the converged outcomes. However, the 2D LVC model
is an exception. Initially, after the first hopping event, the FSSH-2 scheme
yields nearly perfect results, but following the second hopping event, 16
substeps are required to approach the converged population results. For fewer
time steps, significant issues with the conservation of energy were found in
addition to the discrepancies of the population dynamics. To further
investigate this anomaly, we re-ran the simulations for the 2D LVC model, this
time rescaling the momenta along the previous momenta after each hopping
event. While the results for the traditional FSSH method remain largely
unaffected, the FSSH-2 method demonstrates enhanced accuracy: with as few as 4
substeps, convergence is reached, similarly as observed for the other models.
This observation highlights the importance of smaller time step sizes for
precise handling of NACs in FSSH, particularly in multidimensional models.

To quantitatively support our observations, we assessed the maximum substep
number at which the mean population deviation of the second state, calculated
across all main time steps, exceeds a threshold of $0.015$, compared to the
converged results. This threshold indicates the beginning of significant
divergence. The traditional FSSH method typically surpassed this threshold
only at relatively high substep numbers, i. e., 4 for the Tully models and
between 20 and 100 for the other models, see Tab. 2. In marked contrast, the
FSSH-2 method demonstrated superior performance, exceeding the threshold
already for a single substep in all models. The only exception was the 2D LVC
model, where the FSSH-2 method surpassed this threshold only when the number
of substeps was below 8. Hence, the FSSH-2 method showed a relative
improvement over the traditional FSSH method by a factor ranging from
approximately 4 to 32. This substantial increase in efficiency suggests a
potential reduction in running time by a similar factor.

The superior stability of FSSH-2 comes to the forefront especially during the
first hopping event. This is underscored by the graphs in the right panels of
Figs. 1 to 5 which show the population after the first passage of the
crossings or intersections of the lowest two potential energy curves or
surfaces. Especially for the 2D LVC model and Model X, the FSSH-2 method
offers significant advantages, allowing for time step sizes that are around 50
times greater than what FSSH requires.

Overall, these results demonstrate the superior numerical efficiency of our
FSSH-2 method compared to the traditional FSSH method.

Model |  FSSH |  FSSH-2
---|---|---
Tully 1 |  4 |  1
Tully 2 |  4 |  1
2D LVC |  91 |  8
2D Well |  23 |  1
Model X |  32 |  1

Table 2:  Comparison of the FSSH and FSSH-2 methods at different substep
numbers. The second and third columns specify the maximum substep number where
the FSSH and FSSH-2 methods start to display population divergences exceeding
$0.015$.

##  V Conclusions

We presented the FSSH-2 scheme, an enhanced numerical time integration variant
of the original FSSH method. The main advantage of FSSH-2 is that it is
tailored to operate independently of the notorious NAC which may become
(nearly) divergent at (avoided) intersections of adiabatic potential energy
surfaces. This is known to be one of the main causes for the numerical
problems hampering the original FSSH scheme, in particular the restriction to
rather small time steps. Even though our theoretical analysis indicates that
the two methods do not necessarily produce identical outcomes, our simulations
– primarily based on the test suite for non–adiabatic algorithms by Nelson et
al. [23] – revealed congruent results between the FSSH and FSSH-2 methods when
using very small time steps. However, the most important finding in our
analysis is that the FSSH-2 scheme exhibits swifter convergence with respect
to the number of time steps, which makes it approximately 1 to 2 orders of
magnitude faster than conventional FSSH, especially during the initial hopping
event.

## Data Availability Statement

The MATLAB scripts used to generate the results shown in the figures and
tables are openly available in the ZENODO repository at
https://doi.org/10.5281/zenodo.10497421. While we performed the simulations
using version 7.2.0 of WavePacket [24, 25, 19], it is worth mentioning that
FSSH-2 is also available in the Libra package v5.5.0 [29].

###### Acknowledgements.

This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) through project TRR 352 – Project-ID 470903074. Moreover,
L. A. acknowledges support for a research stay at the WIAS in Berlin.

##  Appendix A Simulation Details

Throughout all of our simulations, the initial state is modeled as a Gaussian
wave packet parameterized as follows

$\displaystyle
g(q)=\left(\det(W_{G})\right)^{-\frac{1}{2}}\left(\frac{2}{\pi}\right)^{\frac{d_{\text{nuc}}}{4}}\exp\left(-\frac{1}{4}(q-q_{G})^{T}W_{G}^{-2}(q-q_{G})+ip_{G}^{T}(q-q_{G})\right),$

where $W_{G}\in\mathbb{R}_{>0}^{d_{\text{nuc}}\times d_{\text{nuc}}}$ is a
positive diagonal width matrix, determining the spread of the wave packet in
each dimension. The packet is centered at position
$q_{G}\in\mathbb{R}^{d_{\text{nuc}}}$, with a momentum
$p_{G}\in\mathbb{R}^{d_{\text{nuc}}}$. The initial conditions for the
trajectories are derived from the Wigner transform of the Gaussian wave
packet, as expressed by

$\displaystyle\mathcal{W}_{g}(q,p)=\exp\left(-2(q-q_{G})^{T}W_{G}^{-2}(q-q_{G})-\frac{1}{2}(p-p_{G})^{T}W_{G}^{2}(p-p_{G})\right),$

implying that the sampling follows a normal distribution in $q$ and $p$.

###  A.1 Tully 1

A single crossing model consisting of two electronic states with diabatic
potential matrix given by [1]

$\displaystyle V^{\text{dia}}(q)$
$\displaystyle=\begin{pmatrix}\text{sign}\left(q\right)A\left(1-\exp(-\text{sign}\left(q\right)Bq)\right)&C\exp(-Dq^{2})\\\
C\exp(-Dq^{2})&-\text{sign}\left(q\right)A\left(1-\exp(-\text{sign}\left(q\right)Bq)\right)\end{pmatrix}$

with $A=0.01$, $B=1.6$, $C=0.005$, $D=1.0$, and mass $M=2000$. The initial
Gaussian wave packet has the parameters $W_{G}=0.75$, $q_{G}=-6$, $p_{G}=15$.
Initially, the electronic system occupies the lower adiabatic state. The main
time step size is $200$ a.u. for a total of $9$ main time steps. The main time
step size is $100$ a.u. for a total of $10$ main time steps. The grid for the
quantum solution consisted of $2^{8}$ equally spaced points in the range
$[-10,10]$.

###  A.2 Tully 2

A double crossing model consisting of two electronic states with diabatic
potential matrix given by [1]

$\displaystyle V^{\text{dia}}(q)$
$\displaystyle=\begin{pmatrix}0&C\exp(-Dq^{2})\\\
C\exp(-Dq^{2})&-A\exp(-Bq^{2})+E_{0}\end{pmatrix}$

with $A=0.1$, $B=0.28$, $C=0.015$, $D=0.06$, $E=0.05$, and mass $M=2000$. The
initial Gaussian wave packet has the parameters $W_{G}=0.75$, $q_{G}=-7$,
$p_{G}=30$. Initially, the electronic system occupies the lower adiabatic
state. The main time step size is $100$ a.u. for a total of $10$ main time
steps. The grid for the quantum solution consisted of $2^{9}$ equally spaced
points in the range $[-12,12]$.

###  A.3 2D LVC

A conical intersection model consisting of two electronic states with diabatic
potential matrix given by [20]

$\displaystyle V^{\text{dia}}(q_{1},q_{2})$
$\displaystyle=\frac{1}{2}\begin{pmatrix}\omega_{1}^{2}(q_{1}+\frac{a}{2})^{2}+\omega_{2}^{2}q_{2}^{2}+\sigma&2cq_{2}\\\
2cq_{2}&\omega_{1}^{2}(q_{1}-\frac{a}{2})^{2}+\omega_{2}^{2}q_{2}^{2}-\sigma\end{pmatrix}$

with the parameters used for bis(methylene) adamantyl cation:
$\omega_{1}=7.743\cdot 10^{-3}$, $\omega_{2}=6.680\cdot 10^{-3}$, $a=31.05$,
$c=8.092\cdot 10^{-5}$, $\sigma=0$, and mass $M=1$. The initial Gaussian wave
packet has the parameters
$W_{G}=\operatorname{\text{diag}}(2\omega_{1},2\omega_{2})^{-\frac{1}{2}}$,
$q_{G}=(a/2,0)^{T}$, $p_{G}=(0,0)^{T}$. Initially, the electronic system
occupies the upper adiabatic state. The main time step size is $100$ a.u. for
a total of $11$ main time steps. The grid for the quantum solution consisted
of $2^{8}\times 2^{7}$ equally spaced points in the range
$[-80,80]\times[-40,40]$.

###  A.4 2D Well

A 2D extension of the dual avoided crossing model (Tully 2, see above)
consisting of two electronic states with diabatic potential matrix given by
[21]

$\displaystyle V^{\text{dia}}(q_{1},q_{2})$
$\displaystyle=\begin{pmatrix}\scriptstyle-E_{0}&\scriptstyle
C\exp\left(-D\left(0.25(q_{1}+q_{2})^{2}+0.75(q_{1}-q_{2})^{2}\right)\right)\\\
\scriptstyle
C\exp\left(-D\left(0.25(q_{1}+q_{2})^{2}+0.75(q_{1}-q_{2})^{2}\right)\right)&\scriptstyle-A\exp\left(-B\left(0.75(q_{1}+q_{2})^{2}+0.25(q_{1}-q_{2})^{2}\right)\right)\end{pmatrix}$

with $A=0.15$, $B=0.14$, $C=0.015$, $D=0.06$, $E_{0}=0.05$, and mass $M=2000$.
The initial Gaussian wave packet has the parameters
$W_{G}=\frac{1}{\sqrt{2}}{\rm Id}_{2\times 2}$, $q_{G}=(-8,0)^{T}$,
$p_{G}=(20,0)^{T}$. Initially, the electronic system occupies the lower
adiabatic state. The main time step size is $100$ a.u. for a total of $13$
main time steps. The grid for the quantum solution consisted of $2^{8}\times
2^{7}$ equally spaced points in the range $[-15,15]\times[-10,10]$.

###  A.5 Model X

A triple avoided crossing model consisting of three electronic states with
diabatic potential matrix given by [22]

$\displaystyle\scriptstyle V^{\text{dia}}(q)$
$\displaystyle=\begin{pmatrix}\scriptstyle
A\left(\tanh\left(Bq\right)+\tanh\left(B(q+7)\right)\right)&\scriptstyle
C\exp\left(-q^{2}\right)&\scriptstyle C\exp\left(-(q+7)^{2}\right)\\\
\scriptstyle
C\exp\left(-q^{2}\right)&\scriptstyle-A\left(\tanh\left(Bq\right)+\tanh\left(B(q-7)\right)\right)&\scriptstyle
C\exp\left(-(q-7)^{2}\right)\\\ \scriptstyle
C\exp\left(-(q+7)^{2}\right)&\scriptstyle
C\exp\left(-(q-7)^{2}\right)&\scriptstyle-A\left(\tanh\left(B(q+7)\right)-\tanh\left(B(q-7)\right)\right)\end{pmatrix}$

with $A=0.03$, $B=1.6$, $C=0.005$, and mass $M=2000$. The initial Gaussian
wave packet has the parameters $W_{G}=0.75$, $q_{G}=-12$, $p_{G}=15$.
Initially, the electronic system occupies the central adiabatic state. The
main time step size is $125$ a.u. for a total of $20$ main time steps. The
grid for the quantum solution consisted of $2^{9}$ equally spaced points in
the range $[-15,25]$.

## References

  * Tully [1990] J. C. Tully, Molecular dynamics with electronic transitions, J. Chem. Phys. 93, 1061 (1990).
  * Crespo-Otero and Barbatti [2018] R. Crespo-Otero and M. Barbatti, Recent advances and perspectives on nonadiabatic mixed quantum–classical dynamics, Chemical Reviews 118, 7026 (2018).
  * Ibele and Curchod [2020] L. M. Ibele and B. F. E. Curchod, A molecular perspective on Tully models for nonadiabatic dynamics, Physical Chemistry Chemical Physics 22, 15183 (2020).
  * Jain and Sindhu [2022] A. Jain and A. Sindhu, Pedagogical overview of the fewest switches surface hopping method, ACS Omega 10.1021/acsomega.2c04843 (2022).
  * Freixas  _et al._ [2021] V. M. Freixas, A. J. White, T. Nelson, H. Song, D. V. Makhov, D. V. Shalashilin, S. Fernandez-Alberti, and S. Tretiak, Nonadiabatic excited-state molecular dynamics methodologies: Comparison and convergence, Journal of Physical Chemistry Letters 12, 2970 (2021).
  * Hammes-Schiffer and Tully [1994] S. Hammes-Schiffer and J. C. Tully, Proton transfer in solution: Molecular dynamics with quantum transitions, The Journal of Chemical Physics 101, 4657 (1994).
  * Jain  _et al._ [2016a] A. Jain, E. Alguire, and J. E. Subotnik, An Efficient, Augmented Surface Hopping Algorithm That Includes Decoherence for Use in Large-Scale Simulations, J. Chem. Theory Comput. 12, 5256 (2016a).
  * Granucci  _et al._ [2001] G. Granucci, M. Persico, and A. Toniolo, Direct semiclassical simulation of photochemical processes with semiempirical wave functions, J. Chem. Phys. 114, 10608 (2001).
  * Wang  _et al._ [2014] L. Wang, D. Trivedi, and O. V. Prezhdo, Global flux surface hopping approach for mixed quantum-classical dynamics, J. Chem. Theory Comput. 10, 3598 (2014).
  * Hellmann [1937] H. Hellmann,  _Einführung in die Quantenchemie_ (Deuticke, Leipzig, 1937).
  * Feynman [1939] R. P. Feynman, Forces in molecules, Phys. Rev. 56, 340 (1939).
  * Epstein [1948] S. T. Epstein, A Note on Perturbation Theory, Phys. Rev. 74, 206 (1948).
  * Lasser and Lubich [2020] C. Lasser and C. Lubich, Computing quantum dynamics in the semiclassical regime, Acta Numer. 29, 229 (2020).
  * Hack  _et al._ [1999] M. D. Hack, A. W. Jasper, Y. L. Volobuev, D. W. Schwenke, and D. G. Truhlar, Quantum Mechanical and Quasiclassical Trajectory Surface Hopping Studies of the Electronically Nonadiabatic Predissociation of the Ã State of $\text{NaH}_{2}$, J. Phys. Chem. A 103, 6309 (1999).
  * Jain  _et al._ [2016b] A. Jain, E. Alguire, and J. E. Subotnik, An Efficient, Augmented Surface Hopping Algorithm That Includes Decoherence for Use in Large-Scale Simulations, J. Chem. Theory Comput. 12, 5256 (2016b).
  * Magnus [1954] W. Magnus, On the exponential solution of differential equations for a linear operator, Commun. Pure Appl. Math. 7, 649 (1954).
  * Hairer and Wanner [2006] E. Hairer and G. Wanner,  _Geometric numerical integration algorithms for ordinary differential equations_ (2006) p. 644.
  * Fabiano  _et al._ [2008] E. Fabiano, T. W. Keal, and W. Thiel, Implementation of surface hopping molecular dynamics using semiempirical methods, Chem. Phys. 349, 334 (2008).
  * Schmidt  _et al._ [2019] B. Schmidt, R. Klein, and L. Cancissu Araujo, WavePacket: A Matlab package for numerical quantum dynamics. III: Quantum-classical simulations and surface hopping trajectories, J. Comput. Chem. 40, 2677 (2019).
  * Gherib  _et al._ [2016] R. Gherib, L. Ye, I. G. Ryabinkin, and A. F. Izmaylov, On the inclusion of the diagonal Born-Oppenheimer correction in surface hopping methods, J. Chem. Phys. 144, 154103 (2016).
  * Shenvi  _et al._ [2011] N. Shenvi, J. E. Subotnik, and W. Yang, Phase-corrected surface hopping: Correcting the phase evolution of the electronic wavefunction, J. Chem. Phys. 135, 024101 (2011).
  * Subotnik and Shenvi [2011] J. E. Subotnik and N. Shenvi, Decoherence and surface hopping: When can averaging over initial conditions help capture the effects of wave packet separation?, J. Chem. Phys. 134, 244114 (2011).
  * Nelson  _et al._ [2020] T. R. Nelson, A. J. White, J. A. Bjorgaard, A. E. Sifain, Y. Zhang, B. Nebgen, S. Fernandez-Alberti, D. Mozyrsky, A. E. Roitberg, and S. Tretiak, Non-adiabatic Excited-State Molecular Dynamics: Theory and Applications for Modeling Photophysics in Extended Molecular Materials, Chem. Rev. 120, 2215 (2020).
  * Schmidt and Lorenz [2017] B. Schmidt and U. Lorenz, WavePacket: A Matlab package for numerical quantum dynamics. I: Closed quantum systems and discrete variable representations, Comput. Phys. Commun. 213, 223 (2017).
  * Schmidt and Hartmann [2017] B. Schmidt and C. Hartmann, WavePacket: A Matlab package for numerical quantum dynamics. II: Open quantum systems, optimal control, and model reduction, Comput. Phys. Commun. 228, 229 (2017).
  * Askar and Cakmak [1978] A. Askar and A. S. Cakmak, Explicit integration method for the time-dependent schrodinger equation for collision problems, The Journal of Chemical Physics 68, 2794 (1978).
  * Kosloff and Kosloff [1983] D. Kosloff and R. Kosloff, A Fourier method solution for the time dependent Schrödinger equation as a tool in molecular dynamics, Journal of Computational Physics 52, 35 (1983).
  * Leforestier  _et al._ [1991] C. Leforestier, R. Bisseling, C. Cerjan, M. D. Feit, R. Friesner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H.-D. Meyer, N. Lipkin, O. Roncero, and R. Kosloff, A comparison of different propagation schemes for the time dependent Schrödinger equation, Journal of Computational Physics 94, 59 (1991).
  * Akimov  _et al._ [2023] A. Akimov, M. Shakiba, B. Smith, M. Dutra, K. Sato, DaehoHan, S. Temen, W. Li, T. Khvorost, X. Sun, and L. Stippell, Quantum-Dynamics-Hub/libra-code: Libra with TC- NBRA, MASH, FSSH2, revDISH, revised TSH and Ehrenfest, phase-corrected exact calculations (2023).