# Learning ODEs via Diffeomorphisms for Fast and Robust Integration

Weiming Zhi1 &Tin Lai1 &Lionel Ott2 &Edwin V. Bonilla 3 &Fabio Ramos1,4
1 School of Computer Science, the University of Sydney, Australia
2 Autonomous Systems Lab, ETH Zurich, Switzerland
3 CSIRO’s Data61, Australia
4 NVIDIA, USA Correspondence to<EMAIL_ADDRESS>
###### Abstract

Advances in differentiable numerical integrators have enabled the use of
gradient descent techniques to learn ordinary differential equations (ODEs).
In the context of machine learning, differentiable solvers are central for
Neural ODEs (NODEs), a class of deep learning models with continuous depth,
rather than discrete layers. However, these integrators can be
unsatisfactorily slow and inaccurate when learning systems of ODEs from long
sequences, or when solutions of the system vary at widely different timescales
in each dimension. In this paper we propose an alternative approach to
learning ODEs from data: we represent the underlying ODE as a vector field
that is _related_ to another _base_ vector field by a differentiable
bijection, modelled by an invertible neural network. By restricting the base
ODE to be amenable to integration, we can drastically speed up and improve the
robustness of integration. We demonstrate the efficacy of our method in
training and evaluating continuous neural networks models, as well as in
learning benchmark ODE systems. We observe improvements of up to two orders of
magnitude when integrating learned ODEs with GPUs computation.

##  1 Introduction

The problem of fitting an ordinary differential equation (ODE) to observed
data appears in many disciplines such as biology, chemistry and physics [35,
41]. In the context of machine learning and, in particular, neural networks,
this problem arises within the framework of _Neural ODEs_ [7], a family of
continuous-depth deep learning models that parameterize the dynamics of hidden
states using neural networks. Recent developments in learning ODEs generally
use differentiable adaptive step-size numerical integrators to learn these
dynamics [9, 32].

To accurately integrate ODEs with rapidly varying solutions over time, the
step sizes taken by the numerical integrator can become exceedingly small.
Furthermore, to roll out long sequences of learned ODEs, the neural network
dynamics model is queried sequentially at each step. This can be
unsatisfactorily slow for time critical applications, such as those in robot
control, and is known to suffer from numerical instabilities [17, 8],
especially if the ODE exhibits stiffness. With these challenges in mind, we
propose an alternative approach to learning ODEs from data. We view the
desired target ODE as a vector field that is a “morphed” version of an
alternative base vector field via a _diffeomorphism_ , i.e., a bijective
mapping where both the mapping and its inverse are differentiable. Thus,
instead of directly modelling the desired dynamics with a neural network, we
use an invertible neural network to learn the diffeomorphism that relates the
target ODE to the base ODE. Crucially, the base ODE is much more amenable to
integration and, therefore, we obtain a solution to this simpler ODE and pass
it through the bijection to obtain a solution to the (more complex) target
ODE. Figure 1 shows an example of related vector fields.

Figure 1: Related vector fields can be thought of as a vector field that is
“morphed” into another. (Left) Integral curves (red) of the vector field of a
Linear ODE overlaid on grid points (blue); (Right) Morphed integral curves and
grid.

We investigate the use of base ODEs that are (1) linear or (2) modelled by a
neural network, allowing for integration with accurate numerical integrators
at reasonably large step-sizes. Restricting the base ODE to be linear allows
the use of closed-forms solutions, providing major speed benefits. In this
setup, when integrating trajectories of reasonable length, at no degradation
of performance, we can achieve a drastic speedup of up to two orders of
magnitude when employing GPUs, when compared against existing differentiable
integrators with standard settings. We can also restrict the learned ODE to be
provably asymptotically stable by adding simple constraints to the linear base
ODE. Alternatively, when additional flexibility is required, we remove the
restrictive assumption of a linear base ODE and model the dynamics using a
neural network. We empirically show that in this setup, we can improve the
performance of learning stiff ODEs compared to existing differentiable
integrators, even when we use a much smaller neural network for the base ODE.

In summary, our main contributions are:

  1. 1.

a novel paradigm to learn ODEs from data: invertible neural networks are
trained to morph the target ODE to an alternative base ODE that can be more
tractably integrated;

  2. 2.

analysis of the base as (i) a linear ODE and (ii) a non-linear ODE with
learned dynamics and restricted to take large steps. In the linear case, we
demonstrate how to obtain closed-form integrals, providing significant speed-
ups. In the non-linear case, we demonstrate learning challenging and stiff
ODEs, with smaller networks;

  3. 3.

a principled method to enforce asymptotic stability of learned ODEs, by adding
restrictions to the base ODE.

Proofs and additional details can be found in the supplement.

##  2 Related Work

Our work addresses the problem of learning ODEs that match observations which
is typically encountered when learning with a class of continuous-depth neural
network model, _neural ODEs_ [7]. Our approach also leverages learnable
invertible function approximators in a manner that resembles normalising
flows. Here, we briefly outline the related work around learning of ODEs,
Neural ODEs, invertible neural networks and normalising flows.

Neural ODEs and Learning of ODEs: Neural ODEs are neural network models which
model the hidden state as continuous ODEs rather than discrete layers [7, 30].
Chen et al. proposed to model the dynamics of this information flow by another
neural network, which is trained by differentiating integrators via the
adjoint method [36] for tractable memory usage. More memory efficient
alternatives have also been introduced [45]. In particular, neural network
models which incorporate an ODE, such as latent ODEs [40] and neural CDEs
[25], have found application in time-series tasks. Subsequent strategies have
been introduced to improve the training of these models, including augmenting
the ODE state-space [13], regularisation techniques [16], and hyper-network
extensions [44, 8]. Gaussian processes have also been used in earlier work to
model differential equations [42, 37]. At the core of neural ODE models is a
differentiable integrator which is used to learn the underlying ODE. Our
proposed approach improves the learning of the underlying ODE, and is
compatible with models that incorporate learnable ODEs. We shall empirically
show in section 5.1.3 our method used to learn dynamics within continuous-
depth models. The term “neural ODE” has typically been used in the literature
[44, 8] to refer to neural networks that incorporate ODEs, including the
original work in [7]. However, “neural ODE” has occasionally been used to
refer to an ODE with dynamics parameterised by a neural network [31]. We
follow the former convention.

Invertible neural networks and Normalising Flows: Invertible neural networks
(INN) are a class of function approximators that learn bijections where the
forward and inverse mapping and their Jacobians can be efficiently computed
[5]. INN are typically constructed by invertible building blocks, such as
those introduced in [26, 11, 14]. Advances in INNs are largely motivated by
normalising flows [38, 33], an approach to construct a flexible probability
distribution by finding a differentiable bijection, or diffeomorphism, between
the target distribution and a base distribution. Our approach is similar in
spirit to normalising flows, as we analogously aim learn a diffeomorphism that
relates the vector fields of the target ODE and some base ODE. However, unlike
normalising flows, we do not require the burdensome computation of Jacobian
determinants [23]. A separate line of work, broadly characterised as
_continuous normalising flows_ , use ODEs to build invertible approximators
[18, 7]. Our work proposes the opposite where invertible approximators are
used to learn ODEs.

##  3 Preliminaries

We shall introduce learning Ordinary differential equations (ODEs) with neural
networks. We briefly present background on invertible neural networks, a core
component of our method. We then describe the notions of tangent spaces and
pushforwards, which will be used to elaborate our method.

###  3.1 Learning ODEs with neural networks

ODEs are central in neural ODEs, a class of continuous-depth neural network
models. Many problems in science and engineering can also be described by
ordinary differential equations of the form:

$\displaystyle\mathbf{y}^{\prime}(t)=f(\mathbf{y}(t),t),$
$\displaystyle\mathbf{y}(0)=\mathbf{y}_{0},$ (1)

where $t$ is time, $\mathbf{y}(t)$ are the states at time $t$, and $f$
provides the dynamics. We use a neural network $f_{\omega}$ with parameters
$\omega$ to model the dynamics. We shall henceforth drop the explicit
dependence on time, and consider the autonomous ODEs given by
$\mathbf{y}^{\prime}(t)=f(\mathbf{y}(t))$. Non-autonomous ODEs, which
explicitly depend on time, can be equivalently expressed as autonomous ODEs by
adding a dimension to the states $\mathbf{y}$. Although there have been
attempts to explicitly condition the network weights on time [10], the
dominant approach in the neural ODE literature is simply to append time to the
network input. For an initial condition $\mathbf{y}_{t_{0}}$ at start $t_{0}$,
and some end time $t_{e}$, a solution of the ODE can be evaluated by a
numerical integrator (ODESolve), such as Euler’s or Runge-Kutta methods [6]:

$\displaystyle\mathbf{y}(t_{e})=\mathbf{y}_{t_{0}}+\int^{t_{e}}_{t_{0}}f_{\mathbf{\omega}}(\mathbf{y}(t))\mathrm{d}t=\mathrm{ODESolve}(f_{\mathbf{\omega}},\mathbf{y}_{t_{0}},t_{e}).$
(2)

The learning problem involves estimating, with $f_{\omega}$, the dynamics of
the ODE, provided $n$ observations
$\mathbf{y}^{obs}_{t_{1}}\ldots\mathbf{y}^{obs}_{t_{n}}$ at specified times.
We can learn the ODE by optimising the parameters $\omega$ to minimise a loss
between the observations at the given times and the integrated ODE,
$\ell(\mathbf{\omega})=\mathrm{Loss}(\mathbf{y}^{obs}_{t_{i}},\mathbf{y}(t_{i}))$
for $i=1,\ldots,n$. Advances in the neural ODE literature have introduced
differentiable numerical integrators, which allow gradient descent
optimisation to be applied. By using the adjoint sensitivity method as
outlined in [7], the gradients of adaptive integrators can be obtained in a
memory tractable manner, without differentiating through the operations of the
integrator.

###  3.2 Invertible Neural Networks

A key building block of our method is the invertible neural network (INN).
Invertible neural networks are a class of function approximators which learn
differentiable bijections. INNs can be trained on a forward mapping, and get
the inverse mapping with no further labour, owing to the definition of their
architecture. Throughout this paper, we use INNs of the type described in
[11]. The basic unit is a reversible block, where inputs are split into two
halves, $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$, and the outputs
$\mathbf{v}_{1}$ and $\mathbf{v}_{2}$ are:

$\displaystyle\mathbf{v}_{1}=\mathbf{u}_{1}\odot\exp(s_{2}(\mathbf{u}_{2}))+t_{2}(\mathbf{u}_{2}),$
$\displaystyle\mathbf{v}_{2}=\mathbf{u}_{2}\odot\exp(s_{1}(\mathbf{u}_{1}))+t_{1}(\mathbf{u}_{1}),$
(3)

where $\odot$ indicates element-wise multiplication, and $t_{1}$, $t_{2}$ and
$s_{1}$, $s_{2}$ are functions modelled by fully-connected neural networks
with non-linear activations. These expressions are clearly invertible:

$\displaystyle\mathbf{u}_{1}=(\mathbf{v}_{1}-t_{2}(\mathbf{u}_{2}))\odot\exp(-s_{2}(\mathbf{u}_{2})),$
$\displaystyle\mathbf{u}_{2}=(\mathbf{v}_{2}-t_{1}(\mathbf{u}_{1}))\odot\exp(-s_{1}(\mathbf{u}_{1})).$
(4)

The functions $t_{1}$, $t_{2}$ and $s_{1}$, $s_{2}$ themselves are not
required to be invertible, and can be modelled by regular neural networks.
Although further developments have given rise to methods of constructing
better performing invertible blocks, such as those in [14], the coupling-based
block introduced here is fast, sufficiently flexible, and has been shown to be
a universal diffeomorphism approximator [43].

###  3.3 Tangent Spaces and Pushforwards

Here we briefly introduce the differential geometry notions of tangent spaces
and pushforwards, which will be later used to elaborate on geometric concepts
of our method.

Tangent Spaces: A manifold is a space that locally resembles Euclidean space.
Throughout this paper, all manifolds will be assumed to be differentiable,
with defined _tangent spaces_. For an $n$-dimensional manifold $\mathcal{M}$,
at a point $\mathbf{p}\in\mathcal{M}$, the tangent space
$T_{\mathbf{p}}\mathcal{M}$ is an $n$-dimensional real vector space, where
each element passes $\mathbf{p}$ tangentially and is referred to as a _tangent
vector_. The tangent space provides a higher-dimensional analogue of a tangent
plane at a point on a surface. The collection of tangent spaces for all points
on $\mathcal{M}$ is known as the _tangent bundle_ denoted by $T\mathcal{M}$.

Pushfoward: For a mapping $F:\mathcal{M}\rightarrow\mathcal{N}$ between two
manifolds, $\mathcal{M}$ and $\mathcal{N}$, the _pushforward_ by $F$ is a
linear mapping between the tangent spaces of the manifolds,
$D_{\mathbf{p}}F:T_{\mathbf{p}}\mathcal{M}\rightarrow
T_{F(\mathbf{p})}\mathcal{N}$. Tangent vectors at $\mathbf{p}$ in the domain
$\mathcal{M}$ can be mapped to tangent vectors at the corresponding point
$F(\mathbf{p})$ in the codomain $\mathcal{N}$ via the pushforward. This can be
computed by the matrix product of the Jacobian of $F$ at $\mathbf{p}$ and a
tangent vector at $\mathbf{p}$.

##  4 Methodology

To elaborate on our method, we shall study ODEs as vector fields. The dynamics
of an ODE, $f$, can be viewed as a vector field, and solutions as integral
curves of the vector field [27]. First, we introduce the learning of vector
fields which are _related_ to another _base_ vector field. Second, we describe
possible choices for such _base_ vector fields.

(a) If $F$ maps points $\mathbf{p}\in\mathcal{M}$ to
$F(\mathbf{p})\in\mathcal{N}$, a single tangent vector at $\mathbf{p}$,
$X_{\mathbf{p}}\in T_{\mathbf{p}}\mathcal{M}$, can be mapped to
$T_{\mathbf{F(\mathbf{p})}}\mathcal{N}$. However, an entire vector field $X$
on $\mathcal{M}$ cannot in general be mapped to a valid vector field on
$\mathcal{N}$. The pushforward by a diffeomorphism is a special case where a
valid vector field can be obtained.

${T\mathcal{M}}$${T\mathcal{N}}$${\mathcal{M}}$${\mathcal{N}}$$DF$$X$$F$$F^{-1}$$Y$
(b) If vector fields $X$ and $Y$ on manifolds $\mathcal{M}$ and $\mathcal{N}$
respectively are _related_ by diffeomorphism $F$, then they are related via
the pushforward of $F$. If $Y$ is unknown, we have another path to evaluate
$Y$ by
$\mathcal{N}\xrightarrow[]{F^{-1}}\mathcal{M}\xrightarrow[]{X}T\mathcal{M}\xrightarrow{DF}T\mathcal{N}$.

Figure 2: Vector fields related by a diffeomorphism
$F:\mathcal{M}\rightarrow\mathcal{N}$ allow the pushforward of $F$ to act on
vector fields on $\mathcal{M}$ producing vector fields on $\mathcal{N}$.

###  4.1 Related Vector Fields for ODE Learning

A vector field $X$ defined on manifold $\mathcal{M}$ is a function that
assigns a tangent vector $X_{\mathbf{p}}\in T_{\mathbf{p}}\mathcal{M}$ to each
point $\mathbf{p}\in\mathcal{M}$. Intuitively, our aim is to construct a
mapping $F$ which shapes the manifold where a _base_ vector field $X$ is
defined, such that the pushforward of $X$ by $F$ extrinsically appears
“morphed” to match the data. A question that arises is: _what are the
requirements of these mappings, for the “pushed forward” vector field to be
valid?_

Provided a mapping between manifolds $F:\mathcal{M}\rightarrow\mathcal{N}$, we
can push a single vector, $X_{\mathbf{p}}\in T_{\mathbf{p}}\mathcal{M}$, to
the tangent space of $\mathcal{N}$ at $F(\mathbf{p})$,
$T_{F(\mathbf{p})}\mathcal{N}$, via the _pushforward_ ,
$D_{\mathbf{p}}F(X_{\mathbf{p}})$. Figure 2(a) sketches out an example of a
“pushing forward” a tangent vector between tangent spaces. However, this
notion does not extend in general to vector fields. If mapping $F$ is
injective and non-surjective, the pushforward of $X$ outside the image of $F$
is not defined. On the other hand, if mapping $F$ is surjective and non-
injective, there may be multiple differing pushforwards given for a point. In
special cases when the pushforward of mapping $F$ defines a valid vector field
on the codomain $\mathcal{N}$, the vector field and its pushforward are known
to be $F$-related.

######  Definition 4.1 (Related vector fields).

Let $F:\mathcal{M}\rightarrow\mathcal{N}$ be a smooth mapping of manifolds. A
vector field $X$ on $\mathcal{M}$ and a vector field $Y$ on $\mathcal{N}$ are
related by $F$, or _$F$ -related_, if for all $\mathbf{p}\in\mathcal{M}$,

$D_{\mathbf{p}}F(X_{\mathbf{p}})=Y_{F(\mathbf{p})}$ (5)

Related vector fields arise in particular when $F$ is a bijective mapping,
where both the mapping itself and its inverse are differentiable, i.e. a
_diffeomorphism_.

######  Proposition 4.1 (Proposition 8.19 in [27]).

Suppose $F:\mathcal{M}\rightarrow\mathcal{N}$ is a diffeomorphism between
smooth manifolds $\mathcal{M}$, $\mathcal{N}$. For every vector field $X$ on
$\mathcal{M}$, there is a unique vector field $Y$ on $\mathcal{N}$ that is
$F$-related to $X$.

By considering the $F$-related properties of vector fields, we have a pathway
to define unknown vector fields using the pushforward of $F$, as shown in fig.
2(b). If vector field $X$ on $\mathcal{M}$ is $F$-related to some vector field
$Y$ on $\mathcal{N}$, instead of directly evaluating the vector field $Y$, we
can instead obtain tangent values for any $\mathbf{q}\in\mathcal{N}$, via
$\mathcal{N}\xrightarrow[]{F^{-1}}\mathcal{M}\xrightarrow[]{X}T\mathcal{M}\xrightarrow{DF}T\mathcal{N}$.
Therefore, the vector attached by vector field $Y$ is,

$Y_{\mathbf{q}}=D_{F^{-1}(\mathbf{q})}F(X_{F^{-1}(\mathbf{q})})=J_{F}(F^{-1}(\mathbf{q}))X_{F^{-1}(\mathbf{q})}\text{,
for each }\mathbf{q}\in\mathcal{N},$ (6)

where $J_{F}$ is the Jacobian of diffeomorphism $F$. One would naturally ask:
_why would it be beneficial to construct a desired vector field $Y$ in the
form of eq. 6?_

We shall answer this by considering integral curves on $Y$, which represent
solutions to the ODE associated with $Y$. An integral curve of $Y$ on
$\mathcal{N}$ is a differentiable curve
$\mathbf{y}:\mathbb{R}\rightarrow\mathcal{N}$, whose velocity at each point is
equal to the value of $Y$ at that point, i.e.
$\mathbf{y}^{\prime}(t)=Y_{\mathbf{y}(t)}\in T_{\mathbf{y}(t)}\mathcal{M}$,
for all $t\in\mathbb{R}$. The integral curves of $F$-related vector fields are
also linked by $F$: integral curves on one vector field are mapped to the
other via a single pass through $F$, and evaluation of Jacobian $J_{F}$ is not
required.

######  Proposition 4.2 (Proposition 9.6 in [27]).

Suppose $X$ and $Y$ are vector fields on manifolds $\mathcal{M}$ and
$\mathcal{N}$ respectively. $X$ and $Y$ are related by mapping
$F:\mathcal{M}\rightarrow\mathcal{N}$ if and only if for each integral curve
$\mathbf{x}:\mathbb{R}\rightarrow\mathcal{M}$, $\mathbf{y}=F(\mathbf{x})$ is
an integral curve of $Y$.

In the ODE learning problem outlined in section 3.1, during both training and
inference, we need to obtain an integral curve $\mathbf{y}$ of the vector
field $Y$ either by numerical integration on $Y$, or by
$\mathbf{y}=F(\mathbf{x})$, where $\mathbf{x}$ denotes the corresponding
integral curve of $X$, related to $Y$ via the diffeomophism $F$. The Jacobian
of $F$ does not need to be evaluated when we are working with the integral
curves.

If integral curves of $X$ can be found in a more efficient, or less error-
prone manner, than by numerically integrating curves of $Y$, we can leverage
the relationship $\mathbf{y}=F(\mathbf{x})$ for ODE learning. This can be done
by an invertible neural network, with parameters $\theta$, to parameterise a
diffeomorphism, $F_{\theta}$, which maps between $\mathbb{R}^{n}$ and some
$n$-dimensional manifold that is diffeomorphic to Euclidean space. We denote
the base ODE associated with the $F$-related vector field of our target as
$\mathbf{x}^{\prime}(t)=g_{\varphi}(\mathbf{x}(t))$, where $\varphi$
parameterises the base ODE. We can then use the target ODE within some
learning problem, such as minimising the mean absolute error (MAE) between the
target ODE and observations:

$\ell(\theta,\varphi)=MAE\Big{(}\mathbf{y}_{t_{i}}^{obs},F_{\theta}(F_{\theta}^{-1}(\mathbf{y}_{0})+\int^{t_{i}}_{0}g_{\varphi}(\mathbf{x}(t))\mathrm{d}t)\Big{)},$
(7)

where $\mathbf{y}_{0}$ is an initial condition for the system, and
$\mathbf{y}_{t_{i}}^{obs}$ are observed data points at times $t_{i}$. By eq.
6, the learned ODE can then be written as
$\mathbf{y}^{\prime}(t)=J_{F_{\theta}}(F_{\theta}^{-1}(\mathbf{y}(t)))g_{\varphi}(F_{\theta}^{-1}(\mathbf{y}(t)))$.
Provided an initial solution $\mathbf{y}_{0}$, we can solve via
$\mathbf{y}(t)=F_{\theta}(F_{\theta}^{-1}(\mathbf{y}_{0})+\int^{t}_{0}g_{\varphi}(\mathbf{x}(t))\mathrm{d}t)$.
In practice, we are often required to evaluate an entire trajectory, ie.
$\mathbf{y}(t)$ at multiple times $t_{1},\ldots,t_{end}$ with one initial
$\mathbf{y}_{0}$, as outlined in Algorithm 1. This allows us to batch up the
pass through $F_{\theta}$, which makes this highly efficient when executed on
a GPU. The benefits of our method are apparent when it is advantageous to
integrate the base ODE and then pass the solution through the diffeomorphism,
$F_{\theta}$, rather than numerically integrate the target ODE. Next, we
investigate two choices of base ODE: (1) Linear ODE; (2) Non-linear ODE
amenable to numerical integration.

Input : $F_{\theta}$, $g_{\varphi}$, $\mathbf{y}_{0}$, $t_{1},\ldots,t_{end}$

Output : $\mathbf{y}(t_{1}),\ldots,\mathbf{y}(t_{end})$

1 $\mathbf{x}_{0}\leftarrow F^{-1}_{\theta}(\mathbf{y}_{0})$

$\mathbf{x}(t_{i})\leftarrow\mathbf{x}_{0}+\int^{t_{i}}_{0}g_{\varphi}(\mathbf{x}(t))\mathrm{d}t$,
for $i=1,\ldots,end$ ;

// The integral either admits a closed-form solution or is easily numerically
integrated.

$\mathbf{y}(t_{1}),\ldots,\mathbf{y}(t_{end})\leftarrow
F_{\theta}(\mathbf{x}(t_{1}),\ldots,\mathbf{x}(t_{end}))$ ;

// The pass through the INN can be batched, and efficiently computed on GPUs.

2

Algorithm 1 Efficient integration of learned ODEs

###  4.2 Linear ODE as Base: Fast Integration and Easy Enforcement of
Asymptotic Stability

We can speed-up integration significantly by modelling the base as a Linear
ODE. Linear ODEs are in the form of $\mathbf{x}^{\prime}(t)=A\mathbf{x}(t)$,
where $\mathbf{x}(t)\in\mathbb{R}^{n}$ are $n$-dimensional variables, and
$A\in\mathbb{R}^{n\times n}$. Linear ODEs can be solved very efficiently as
they admit closed-form solutions. Provided an initial solution
$\mathbf{x}_{0}$, the solution of $\mathbf{x}(t)$ and the target ODE
$\mathbf{y}(t)$ are then

$\displaystyle\mathbf{x}(t)=\sum_{k=1}^{n}(\mathbf{l}_{k}\cdot\mathbf{x}_{0})\mathbf{r}_{k}\exp(\lambda_{k}t),$
$\displaystyle\mathbf{y}(t)=F_{\theta}(\mathbf{x}(t)),$ (8)

where $\mathbf{l}_{k}$, $\mathbf{r}_{k}$ and $\lambda_{k}$ are the
corresponding left, right eigenvectors and eigenvalue of matrix $A$
respectively. We learn the eigenvalues and eigenvectors or matrix $A$ jointly
with diffeomorphism $F_{\theta}$.

Linear ODEs are also interesting because their long-term behaviour, which is
determined by their eigenvalues, is easy to analyse. We shall see how this
property allows us to craft the long-term behaviour of the desired target ODE.
In particular, in many applications, consideration is given to the asymptotic
properties of ODEs, namely what happens to the solutions after a long period
of time. Will the solution converge to equilibrium points, periodic orbits, or
diverge and fly off? Non-linear ODEs with neural network dynamics typically do
not restrict the learned ODE to be provably stable. Learning ODEs with our
method provides a straightforward way to restrict the ODE to be asymptotically
stable. We begin by defining equilibrium points and asymptotic stability of
first order ODEs.

######  Definition 4.2 (Equilibrium point).

An equilibrium point $\mathbf{y}^{*}$ of an ODE
$\mathbf{y}^{\prime}(t)=f(\mathbf{y}(t))$, is a point where
$f(\mathbf{y}^{*})=0$.

######  Definition 4.3 (Asymptotic stability).

An ODE $\mathbf{y}^{\prime}(t)=f(\mathbf{y}(t),t)$ is asymptotically stable if
for every solution $\mathbf{y}(t)$, there exists a $\delta>0$, such that
whenever $\lvert\lvert\mathbf{y}(t_{0})-\mathbf{y}^{*}\lvert\lvert<\delta$,
then $\mathbf{y}(t)\rightarrow\mathbf{y}^{*}$ as $t\rightarrow\infty$, where
$\mathbf{y}^{*}$ is some equilibrium point.

Intuitively, asymptotically stable systems of ODEs will always settle at some
equilibrium points after a long period of time. In the context of vector
fields related by a diffeomorphism, the asymptotic stability properties of the
ODEs are shared.

######  Theorem 4.1.

Suppose two ODEs $\mathbf{x}^{\prime}(t)=g(\mathbf{x}(t))$,
$\mathbf{y}^{\prime}(t)=f(\mathbf{y}(t))$ are related via
$\mathbf{y}(t)=F(\mathbf{x}(t))$, where $F$ is a diffeomorphism. If the former
ODE is asymptotically stable with $n_{e}$ equilibrium points
$\mathbf{x}^{*}_{1},\ldots,\mathbf{x}^{*}_{n_{e}}$, then the latter ODE is
also asymptotically stable, with equilibrium points
$F(\mathbf{x}^{*}_{1}),\ldots,F(\mathbf{x}^{*}_{n_{e}})$.

Therefore, if we can restrict the base ODE to be asymptotically stable, then
the target ODE learned by our method is also asymptotically stable. When the
base is an $n$ dimensional linear ODE, we can restrict it to be asymptotically
stable by directly learning the eigenvalues, $\lambda_{i}$ for $i=1,\ldots,n$,
and constraining them to be negative, i.e. $\lambda_{i}<0$ for $i=1,\dots,n$.
This can be done by $\lambda_{i}=-(s_{\lambda_{i}})^{2}-\varepsilon I$, where
$\varepsilon$ is a small positive constant, and learning $s_{\lambda_{i}}$
instead of learning the eigenvalues.

###  4.3 Non-linear ODE as Base: Improved Robustness for ‘Difficult’ ODEs

Using linear systems as base ODEs provides a dramatic increase in speed at the
cost of flexibility. We observe that the computation overhead of a single
backward pass $F^{-1}_{\theta}$ and a batched single forward pass $F_{\theta}$
is minimal when compared with numerical integration. When the ODE is difficult
to learn and only a moderate speed-up is required, we can also parameterise
the dynamics of the base ODE using a neural network. This is particularly
appealing for ODEs which are considered stiff.

Although there is no precise definition of stiffness, a common characteristic
of stiff ODEs is a rapid varying of the solution in time at different orders
of magnitude across state dimensions. Existing differentiable explicit
integrators are unable to learn these ODEs. Learning the target ODE will
require step-sizes that are exceedingly small at even large tolerances.
Adaptive step-size solvers will encounter arithmetic underflow, while fixed
step-size solvers are struggle to accurately integrate [19]. We restrict the
base ODE to be solved with an Euler method integrator, and learn it jointly
with the diffeomorphism. The burden of accurately representing the stiff
dynamics is shared by $F_{\theta}$, providing added flexibility. The
diffeomorphism is observed to learn to relate the target ODE to an ODE that is
amenable to integration, where the timescales across dimensions do not differ
greatly at the same time. This set-up allows us to learn stiff systems to an
accuracy that cannot be achieved by typical integrators, even when the neural
network model of the base ODE is much smaller than the neural network used to
directly learn the target ODE, providing a speed-up during integration.

##  5 Experimental Results

We empirically evaluate the ability of our method to speed up the integration
of learned ODEs, along with the robustness of integration when learning
potentially stiff ODEs. Throughout this section, we compare the error and
integration times of our method against a variety of solvers. Including fixed
step-size solvers: Euler’s, midpoint, and Runge-Kutta 45 (RK4), and the
adaptive step-size solvers Dormand–Prince 5 (DOPRI5) and Dormand–Prince 8
(DOPRI8). For all fixed step-size solvers, we set the step-size equal to the
smallest time increment for which we require outputs. We augment the ODE
states in accordance to [13] in all of the ODEs trained during the
experiments, except when recreating results from [40] in section 5.1.3, where
we use the implementation from the original authors. For all adaptive step-
size solvers, we set absolute and relative tolerances to $10^{-5}$. The
differentiable solvers are implemented in the _torchdiffeq_ library, with
neural networks implemented with Pytorch [34]. Additional details on
experimental setup are available in the supplement.

###  5.1 Substantial Integration Speed-up by Learning with a Linear ODE Base

We test our hypothesis that the availability of a closed-form expression for
the integral, when using a linear base ODE, provides a substantial speed-up of
integration. We evaluate on learning synthetic ODE systems, real-world robot
demonstrations, and within a _Latent ODE_ [40]. Here, we report performance
and integration times. Training times can be found in the supplement.

Figure 3: Learning the 3D Lotka–Volterra. (Left) Interpolating (red) data
(blue); (Right) Generating trajectories (red) at unseen initial conditions and
the ground truth (blue).

####  5.1.1 Learning 3D Lotka-Volterra

We train and evaluate models on data from the 3D Lokta-Volterra system, which
models the dynamics of predator-prey populations. The data is corrupted by
white noise with standard deviation of $0.05$. We train our model using a
linear base ODE, and assess the capability of our model in interpolating the
data points at $10$x the data resolution, and generalising to $16$ unseen
initial conditions to integrate trajectories, also at $10$x the data
resolution. The integrated trajectories of the same time duration for
interpolation and generalisation to novel initial conditions have similar
integration run-times. Hence, we only report the integration time for
generalisation. Figure 3 shows interpolation results and newly generated
trajectories, where we see that our model is able to capture the dynamics of
the system. Furthermore, table 1 provides a quantitative evaluation, where we
see that our method not only achieves better performance in terms of MSE but
also is significantly faster than the competing numerical integrators with
speed-ups of more than two orders of magnitude.

####  5.1.2 A Time Critical Application: Motion Generation After Learning from
Demonstrations

The ability to quickly roll-out trajectories is crucial in robotics settings.
In particular, we consider the application of generating new motion
trajectories from provided demonstrations. We evaluate the ability of our
method in learning end-effector trajectories of robot manipulators. The goal
is to learn an ODE system where trajectories integrated at different starting
points mimic the shown demonstrations. We use three sets of data from [24]:
(i) 7 demonstrations of drawing “S” shapes on a flat surface; (ii) 14
demonstrations of placing a cube on a shelf; (iii) 12 demonstrations of
drawing out large “C”-shapes. We use $70\%$ of the data for training, and test
our generalisation capabilities on the remaining demonstrations. In these
datasets, the motion converges to equilibrium points. Hence, we constrain the
learned ODE to be asymptotically stable. We report the performance and run-
times of generalising to new starting points in table 1. We see that our
approach is competitive in the quality of generalised trajectories, while
achieving speed-ups of more than two-orders of magnitude.

####  5.1.3 ODE Learning for Continuous Deep Learning Models

We evaluate our method as a component of _Latent ODEs_ [40], a continuous-
depth deep learning model. Latent ODEs embeds the time series observations as
hidden states via an encoder-decoder. An ODE, with dynamics parameterised by a
neural network, is fit on the hidden states, using a differentiable
integrator. By assuming the hidden states follow some continuous time
dynamics, the model is capable of handling irregularly sampled series. This
avoids the need to group observations into equally-timed bins. In our
experiments, we apply our method, with a linear ODE base, to learn the
dynamics governing the hidden states. We report results for reconstructing the
periodic curves and the human activity classification problem, which were used
in the original latent ODE paper [40], as well as an additional ECG
classification problem. We reconstruct the periodic curves for 100 and 1000
time-steps. The performance and times spent on integrating the hidden state
dynamics are reported in table 2. In our comparisons, the latent ODE is set up
according to the original paper. We see that by leveraging the closed-form
expression of integrals of the ODEs, we achieve integration times that are
hundreds of times faster. We also note that the main computation cost of the
integral in our method is the pass through the invertible neural network. The
parallel computation capabilities of GPUs allow us to batch the pass at
constant cost, whereas the sequential nature of numerical integrators result
in a linear increase in run-time, as demonstrated by the differences in
integration times of “Periodic 100” and “Periodic 1000”.

| 3D Lotka-Volterra | Imitation S | Imitation cube pick | Imitation C
---|---|---|---|---
| MSE (I) | MSE (G) | Time (ms) | MSE (G) | Time (ms) | MSE (G) | Time (ms) | MSE (G) | Time (ms)
Ours (Lin) | 0.143± 0.096 | 1.482± 0.11 | 9.32± 0.37 | 6.11±1.21 | 6.62± 0.15 | 18.61± 6.21 | 7.12 ± 1.58 | 8.12± 1.58 | 7.53± 0.76
Euler | 4.460± 0.264 | 4.621± 0.123 | 385.55± 14.42 | 10.30±2.90 | 724.65± 8.32 | 14.92± 1.38 | 728.40± 9.46 | 7.26± 1.99 | 753.94± 1.44
Midpoint | 0.382± 0.047 | 5.514± 0.120 | 670.37± 31.25 | 10.93±3.25 | 581.59± 13.27 | 12.89± 1.25 | 1267.16± 13.58 | 6.87± 2.23 | 1305.44± 14.74
RK4 | 0.354± 0.005 | 5.630± 0.149 | 1316.11± 30.84 | 10.27±2.99 | 2501.67± 18.91 | 15.89± 0.89 | 2522.74± 23.14 | 7.59± 2.64 | 1292.27± 21.95
DOPRI5 | 0.927± 0.050 | 5.186± 0.365 | 264.67± 17.02 | 10.83±2.76 | 1277.71± 14.24 | 14.88± 0.86 | 504.04± 12.33 | 7.13± 1.92 | 623.38± 15.57

Table 1:  The mean squared error for interpolation, MSE (I), and
generalisation, MSE (G), and mean execution times (± 1 standard deviations) on
the 3D Lotka-Volterra system and the time critical application of motion
trajectory generation from demonstrations for our method using a linear base
ODE and competing numerical integrators.

| Periodic 100 | Periodic 1000 | Human Activity | ECG
---|---|---|---|---
| MSE | Int. time (ms) | MSE | Int. time (ms) | Acc. | Int. time (ms) | Acc. | Int. time (ms)
Ours (Lin) | 0.030 | 2.7± 0.6 | 0.008 | 2.8± 0.8 | 0.864 | 4.2± 1.8 | 0.966 | 7.6± 2.5
Euler | 0.040 | 33.7± 2.6 | 0.043 | 326.6± 9.5 | 0.815 | 67.9± 2.9 | 0.963 | 100.0± 2.8
Midpoint | 0.032 | 54.5± 1.8 | 0.074 | 510.1± 15.5 | 0.865 | 114.2± 2.4 | 0.963 | 169.7± 3.8
RK4 | 0.039 | 95.6± 1.6 | 0.052 | 1020.0± 60.0 | 0.857 | 221.2± 4.2 | 0.963 | 325.5± 4.8
DOPRI5 | 0.045 | 83.4± 2.2 | 0.050 | 264.7± 4.6 | 0.869 | 67.9± 5.0 | 0.963 | 123.3± 2.6
DOPRI8 | 0.041 | 99.6± 2.3 | 0.049 | 282.7± 6.4 | 0.724 | 94.8± 1.6 | 0.963 | 171.6± 3.6

Table 2:  The mean squared error and mean integration times (± 1 standard
deviations) when using latent ODEs on the tasks of periodic curve
reconstruction using 100 and 1000 time-steps and the classification problems
of human activity and ECG for our method using a linear base ODE and competing
numerical integrators.  Figure 4: An example of integrating stiff ROBER (log-
scale). (Left) When learning stiff ROBER with RK4, we see that integrated
curves struggle to match the ground truth at dimensions with small scales.
(Right) Learning and integration by our method matches the ground truth
better.

###  5.2 Robust Integration by Learning with a Non-Linear Neural Network Base

We test our hypothesis that using a neural network base ODE allows us to learn
ODEs that are difficult to integrate or stiff. We learn and evaluate models
trained on the chaotic Lorenz system, and the stiff Robertson’s system [39]
(ROBER). A simple differentiable Euler integrator, with step-size equal to the
data time-step size, is used to learn our base ODE. Figure 5(a) illustrates a
generated trajectory from the Lorenz system at $10$x the data resolution,
which closely resembles the data points (blue). For the stiff ROBER system,
attempts to learn the system with adaptive step solvers lead to numerical
errors, while fixed-step size solvers struggle to accurately integrate as the
dimensions operate on very different scales. The results of RK4 and our method
in learning the ROBER system are shown in fig. 4. Note that the y-axis of the
plots are in log-scale. The trajectory obtained from RK4 can match the
dimension that operates in the $10^{4}$ scale, but cannot concurrently
adequately account for the trends in the other dimensions, which operate in
much smaller scales. In our method, we observe the learned base ODE has
dimensions that are much closer in scale, allowing for easier integration with
large steps. The integrated base trajectories can then be passed through
diffeomorphism $F$ to recover the integrated target trajectories. Figure 5
(right) provides the performance and integration times of learning with our
method and baseline numerical integrators, where we see that our method is
more accurate than competing approaches while also requiring smaller execution
times.

(a) Lorenz system. Data (blue) with generated trajectory (red).

| Lorenz | ROBER
---|---|---
| MAE | Time(ms) | MAE | Time(ms)
Ours (Non-Lin) | 0.131 | 403.62± 64.89 | 1.33 | 24.60± 1.83
Euler | 10.99 | 455.67± 13.11 | 2.03 | 26.35± 0.57
Midpoint | 6.60 | 805.05± 42.95 | 2.49 | 48.89± 3.52
RK4 | 6.81 | 1760.84± 205.54 | 1.60 | 98.68± 4.65
DOPRI5 | 7.55 | 632.22± 83.33 | N/A | N/A

Figure 5:  Results on chaotic (Lorenz) and stiff (ROBER) systems. Left: a
generated trajectory by our method along with the true underlying data. Right:
Quantitative evaluation of our method using a non-linear base ODE and
competing numerical integrators.

##  6 Conclusions, Limitations and Societal Impacts

We have proposed a novel approach to learning ODEs with neural networks,
speeding up integration times by up to two orders of magnitude. Our method
uses invertible neural networks to learn a diffeomorphism relating a desired
target ODE to a _base_ ODE that is easier to integrate. We have shown that the
base ODE can be linear or parameterised by a neural network. We have leveraged
the closed form solution of linear ODEs to provide remarkable speed-ups and
allow constraints to be applied on the asymptotic properties of the learned
ODE. We have also shown that, by using a base ODE parameterised by a neural
network, we can learn “difficult” ODEs, with smaller networks to model
dynamics. We have validated our approach by learning ODEs on synthetic and
real-world data, and within continuous-depth neural network models.

Limitations: We are restricted to learning ODE systems that evolve on some
manifold that is diffeomorphic globally to Euclidean space (this includes
Euclidean space). For example, we cannot learn an ODE that evolves on a torus.
Standard ODE integrators are typically restricted to Euclidean space, so our
method is still applicable to scenarios when differentiable numerical
integrators are used. However, recent works have extended differentiable
integrators to ODEs on more general manifolds [29, 15]. An extension of our
approach to address this could be to segment the space up and use multiple
diffeomorphisms.

Potential negative societal impact: The learning of ODEs appears in many
disciplines, including fields that guide decisions around social policy.
Potential negative impacts can be mitigated by ensuring that results from
models are validated by human judgement.

## References

  * [1] ECG heartbeat categorization dataset.  Kaggle Dataset.
  * [2] Pytorch implementation of differentiable ode solvers (torchdiffeq).  GitHub repository, https://github.com/rtqichen/torchdiffeq.
  * [3] Framework for easily invertible architectures (freia).  GitHub repository, https://github.com/VLL-HD/FrEIA.
  * [4] Latent odes for irregularly-sampled time series.  GitHub repository, https://github.com/YuliaRubanova/latent_ode.
  * [5] L. Ardizzone, J. Kruse, S. Wirkert, D. Rahner, E. W. Pellegrini, R. S. Klessen, L. Maier-Hein, C. Rother, and U. Köthe.  Analyzing inverse problems with invertible neural networks.  In ICLR, 2019.
  * [6] J. C. Butcher.  The Numerical Analysis of Ordinary Differential Equations: Runge-Kutta and General Linear Methods.  Wiley-Interscience, 1987.
  * [7] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud.  Neural ordinary differential equations.  Advances in Neural Information Processing Systems, 2018.
  * [8] K. M. Choromanski, J. Q. Davis, V. Likhosherstov, X. Song, J.-J. Slotine, J. Varley, H. Lee, A. Weller, and V. Sindhwani.  Ode to an ode.  In Advances in Neural Information Processing Systems, 2020.
  * [9] T. Daulbaev, A. Katrutsa, L. Markeeva, J. Gusak, A. Cichocki, and I. Oseledets.  Interpolation technique to speed up gradients propagation in neural odes.  In Advances in Neural Information Processing Systems, 2020.
  * [10] J. Q. Davis, K. Choromanski, J. Varley, H. Lee, J.-J. Slotine, V. Likhosterov, A. Weller, A. Makadia, and V. Sindhwani.  Time dependence in non-autonomous neural odes.  ICLR 2020 Workshop DeepDiffEq, 2020.
  * [11] L. Dinh, J. Sohl-Dickstein, and S. Bengio.  Density estimation using real NVP.  In International Conference on Learning Representations, ICLR, 2017.
  * [12] A. Dontchev and R. Rockafellar.  Implicit Functions and Solution Mappings: A View from Variational Analysis.  2009\.
  * [13] E. Dupont, A. Doucet, and Y. W. Teh.  Augmented neural odes.  In Advances in Neural Information Processing Systems, 2019.
  * [14] C. Durkan, A. Bekasov, I. Murray, and G. Papamakarios.  Neural spline flows.  In Advances in Neural Information Processing Systems, 2019.
  * [15] L. Falorsi and P. Forré.  Neural ordinary differential equations on manifolds, 2020.
  * [16] C. Finlay, J.-H. Jacobsen, L. Nurbekyan, and A. Oberman.  How to train your neural ODE: the world of Jacobian and kinetic regularization.  In Proceedings of the 37th International Conference on Machine Learning, 2020.
  * [17] A. Gholami, K. Keutzer, and G. Biros.  Anode: Unconditionally accurate memory-efficient gradients for neural odes.  In IJCAI, 2019.
  * [18] W. Grathwohl, R. T. Q. Chen, J. Bettencourt, I. Sutskever, and D. Duvenaud.  Ffjord: Free-form continuous dynamics for scalable reversible generative models.  ICLR, 2019.
  * [19] E. Hairer, S. P. Nørsett, and G. Wanner.  Solving Ordinary Differential Equations I (2nd Revised. Ed.): Nonstiff Problems.  Springer-Verlag, Berlin, Heidelberg, 1993.
  * [20] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant.  Array programming with NumPy.  Nature, 585(7825), 2020.
  * [21] R. A. Horn and C. R. Johnson.  Matrix Analysis.  Cambridge University Press, 2nd edition, 2012.
  * [22] J. D. Hunter.  Matplotlib: A 2d graphics environment.  Computing In Science & Engineering, 2007.
  * [23] M. Karami, D. Schuurmans, J. Sohl-Dickstein, L. Dinh, and D. Duckworth.  Invertible convolutional flow.  In Advances in Neural Information Processing Systems, 2019.
  * [24] S. M. Khansari-Zadeh and A. Billard.  Learning stable nonlinear dynamical systems with gaussian mixture models.  IEEE Transactions on Robotics, 2011.
  * [25] P. Kidger, J. Morrill, J. Foster, and T. Lyons.  Neural controlled differential equations for irregular time series.  In Advances in Neural Information Processing Systems, 2020.
  * [26] D. P. Kingma, T. Salimans, R. Jozefowicz, X. Chen, I. Sutskever, and M. Welling.  Improved variational inference with inverse autoregressive flow.  In Advances in Neural Information Processing Systems, 2016.
  * [27] J. Lee.  Introduction to smooth manifolds. 2nd revised ed.  Springer, 2012.
  * [28] S. Lefschetz and R. Alverson.  Stability by liapunov’s direct method with applications.  Physics Today, 1962.
  * [29] A. Lou, D. Lim, I. Katsman, L. Huang, Q. Jiang, S. N. Lim, and C. M. De Sa.  Neural manifold ordinary differential equations.  In Advances in Neural Information Processing Systems, 2020.
  * [30] S. Massaroli, M. Poli, J. Park, A. Yamashita, and H. Asama.  Dissecting neural odes.  In Advances in Neural Information Processing Systems, 2020.
  * [31] A. Norcliffe, C. Bodnar, B. Day, J. Moss, and P. Liò.  Neural ode processes.  In ICLR, 2021.
  * [32] A. Norcliffe, C. Bodnar, B. Day, N. Simidjievski, and P. Lió.  On second order behaviour in augmented neural odes.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, 2020.
  * [33] G. Papamakarios, E. T. Nalisnick, D. J. Rezende, S. Mohamed, and B. Lakshminarayanan.  Normalizing flows for probabilistic modeling and inference.  Journal of Machine Learning Research, 2021.
  * [34] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In Advances in Neural Information Processing Systems 32. 2019.
  * [35] L. Perko.  Differential Equations and Dynamical Systems.  Springer-Verlag, 1991.
  * [36] L. S. Pontryagin, E. Mishchenko, V. Boltyanskii, and R. Gamkrelidze.  The mathematical theory of optimal processes.  1962\.
  * [37] M. Raissi, P. Perdikaris, and G. Karniadakis.  Numerical gaussian processes for time-dependent and nonlinear partial differential equations.  ArXiv, 2018.
  * [38] D. Rezende and S. Mohamed.  Variational inference with normalizing flows.  In Proceedings of the 32nd International Conference on Machine Learning, pages 1530–1538, 2015.
  * [39] H. Robertson.  The solution of a set of reaction rate equations.  1966\.
  * [40] Y. Rubanova, R. T. Q. Chen, and D. Duvenaud.  Latent odes for irregularly-sampled time series.  Advances in Neural Information Processing Systems, 2019.
  * [41] L. Ruthotto, S. J. Osher, W. Li, L. Nurbekyan, and S. W. Fung.  A machine learning framework for solving high-dimensional mean field game and mean field control problems.  Proceedings of the National Academy of Sciences, 2020.
  * [42] S. Sarkka.  Bayesian Filtering and Smoothing.  Cambridge University Press, 2013.
  * [43] T. Teshima, I. Ishikawa, K. Tojo, K. Oono, M. Ikeda, and M. Sugiyama.  Coupling-based invertible neural networks are universal diffeomorphism approximators.  In Advances in Neural Information Processing Systems, 2020.
  * [44] T. Zhang, Z. Yao, A. Gholami, J. E. Gonzalez, K. Keutzer, M. W. Mahoney, and G. Biros.  Anodev2: A coupled neural ode framework.  In Advances in Neural Information Processing Systems, 2019.
  * [45] J. Zhuang, N. Dvornek, S. Tatikonda, and J. Duncan.  Mali: A memory efficient and reverse accurate integrator for neural odes.  In ICLR, 2021.

##  Appendix A Additional Results: Training times

We present the training times for directly learning ODEs with our method,
using a linear base ODE. These include the training times on the 3D Lotka-
Volterra, and the robot imitation datasets, outlined in Sections 5.1.1 and
5.1.2 of the main paper. We run training for 1000 iterations, where in each
iteration the batch includes the entire training set. We see that, by
leveraging the closed-form solution of linear ODEs, our method is able to also
drastically speed up training. Additionally the parallel nature of passing
through the invertible neural network allows more consistent training times
across datasets.

| 3D Lotka-Volterra | Imitation S | Imitation cube pick | Imitation C
---|---|---|---|---
| Per iter (s) | Total (s) | Per iter (s) | Total (s) | Per iter (s) | Total (s) | Per iter (s) | Total (s)
Ours | 0.030± 0.002 | 29.57 | 0.029± 0.002 | 29.35 | 0.031± 0.007 | 31.34 | 0.030± 0.004 | 29.78
Euler | 0.131±0.003 | 130.94 | 1.740± 0.017 | 1740.43 | 1.713± 0.016 | 1712.53 | 1.704± 0.010 | 1703.85
Midpoint | 0.235±0.006 | 235.07 | 3.228± 0.027 | 3227.51 | 3.177± 0.028 | 3177.37 | 3.207± 0.025 | 3206.92
RK4 | 0.469±0.005 | 468.52 | 6.671± 0.048 | 6671.44 | 6.388± 0.046 | 6388.10 | 6.441± 0.057 | 6440.61
Dopri5 | 0.408±0.037 | 408.36 | 1.413± 0.034 | 1413.34 | 1.246± 0.023 | 1245.65 | 1.247± 0.022 | 1246.67

Table 3: The training times in seconds with standard deviations, for 1000
iterations. By leveraging the closed-form solution of linear ODEs, training
time with our method is consistently orders of magnitude faster than by using
a differentiable numerical integrator.

##  Appendix B Additional Results: Ablation Study

We study the effects of the number of layers in the invertible neural network
and number of parameters in the sub-network, which are the main hyper-
parameters of the invertible neural networks used. To this end, we conduct
ablation studies of the speed and performance of our method on the real-world
datasets outlined in section 5.1.2 of the paper. Our basic model uses an
invertible neural network with $5$ layers and sub-networks in the invertible
network had $1500$ hidden dimension size. We alter the number of layers to be:
$2$, $3$, $4$, $5$, $6$, $7$, $8$, and hidden dimensions of the sub-networks
within the invertible network to be: $500$, $1000$, $1500$, $2000$, $2500$.
The results are presented below:

|  | Imitation S | Imitation cube pick | Imitation C
---|---|---|---|---
No. Layers | Sub-Net Hid. Dim. Size | Int. time (ms) | MSE | Int. time (ms) | MSE | Int. time (ms) | MSE
2 | 1500 | 3.551± 0.585 | 122.40 | 2.993±0.103 | 41.51 | 2.914±0.049 | 20.69
3 | 1500 | 4.589± 1.193 | 130.49 | 4.026±0.129 | 15.00 | 5.150±1.407 | 26.20
4 | 1500 | 5.418± 0.746 | 24.54 | 6.294±0.939 | 26.18 | 5.374±0.382 | 10.33
5 | 1500 | 6.461± 0.686 | 4.40 | 7.401±1.531 | 26.56 | 7.463±1.929 | 13.16
6 | 1500 | 7.529± 0.698 | 8.17 | 8.993±2.212 | 17.16 | 8.994±2.781 | 18.27
7 | 1500 | 9.426± 1.664 | 4.91 | 9.858±1.828 | 20.39 | 9.669±1.240 | 25.76
8 | 1500 | 10.636± 2.541 | 5.62 | 10.732±2.111 | 14.56 | 10.958±2.518 | 6.57
5 | 500 | 7.159± 1.475 | 5.62 | 8.018±2.109 | 19.37 | 7.315±1.483 | 6.22
5 | 1000 | 6.972± 1.247 | 6.04 | 6.376±0.203 | 11.02 | 7.297±1.311 | 6.37
5 | 1500 | 7.031± 1.289 | 4.40 | 7.321±1.236 | 26.56 | 6.901±0.802 | 13.16
5 | 2000 | 7.787± 1.363 | 6.23 | 7.443±1.457 | 14.05 | 7.776±2.514 | 9.48
5 | 2500 | 7.208± 1.628 | 10.92 | 6.521±0.082 | 11.66 | 7.611±1.687 | 5.59

Table 4: Ablation study results of different configurations for the invertible
neural network model.

We see that as we increase the number of invertible network layers, the
integration times increase, while the hidden dimension size of the sub-
networks within the invertible network does not visibly affect the integration
times. Overall, the generalisation performance improves as the number of
invertible layers are used, up to some number of layers. Beyond this number of
layers, adding layers does not vary performance significantly. Additionally,
the hidden dimension sizes, for the values tested do not greatly vary the
generalisation performance.

##  Appendix C Proofs

Proofs for Propositions 4.1 and 4.2 can be found in [27] as Propositions 8.19
and 9.6.

######  Theorem 4.1.

Suppose two ODEs $\mathbf{x}^{\prime}(t)=g(\mathbf{x}(t))$,
$\mathbf{y}^{\prime}(t)=f(\mathbf{y}(t))$ are related via
$\mathbf{y}(t)=F(\mathbf{x}(t))$, where $F$ is a diffeomorphism. If the former
ODE is asymptotically stable with $n_{e}$ equilibrium points
$\mathbf{x}^{*}_{1},\ldots,\mathbf{x}^{*}_{n_{e}}$, then the latter ODE is
also asymptotically stable, with equilibrium points
$F(\mathbf{x}^{*}_{1}),\ldots,F(\mathbf{x}^{*}_{n_{e}})$.

###### Proof.

First we show $F(\mathbf{x}^{*}_{1}),\ldots,F(\mathbf{x}^{*}_{n_{e}})$ are
equilibrium points of ODE $\mathbf{y}^{\prime}(t)=f(\mathbf{y}(t))$. By
$\mathbf{y}(t)=F(\mathbf{x}(t))$, we can write the time derivatives
$\mathbf{y}^{\prime}$ at $F(\mathbf{x})$ as

$\mathbf{y}^{\prime}(t)=f(F(\mathbf{x}(t)))=\frac{\mathrm{d}F(\mathbf{x}(t))}{\mathrm{d}t}=J_{F}(\mathbf{x}(t))g(\mathbf{x}(t)),$
(9)

where $J_{F(\mathbf{x}(t))}$ is the Jacobian of $F$. $F$ is a diffeomorphism
and hence invertible over its domain. By the inverse function theorem [12],
the Jacobian $J_{F}(\mathbf{x}(t))$ is invertible, and furthermore, by the
invertible matrix theorem [21], it has a null-space containing only the zero
vector. Therefore,
$\mathbf{y}^{\prime}(t)=f(F(\mathbf{x}(t)))=J_{F}(\mathbf{x}(t))g(\mathbf{x}(t))=0$
if and only if $g(\mathbf{x}(t))=0$. As $g(\mathbf{x}^{*}(t))=0$ for
$\mathbf{x}^{*}\in\\{\mathbf{x}_{1}^{*}\ldots\mathbf{x}_{n_{e}}^{*}\\}$, then
we also have $f(F(\mathbf{x}^{*}(t)))=0$, hence
$\mathbf{y}^{*}\in\\{F(\mathbf{x}^{*}_{1}),\ldots,F(\mathbf{x}^{*}_{n_{e}})\\}$
gives equilibrium points for $\mathbf{y}^{\prime}(t)=f(\mathbf{y}(t))$.

We now show asymptotically stability of
$\mathbf{y}^{\prime}(t)=f(\mathbf{y}(t))$, by the existence of a _Lyapunov
function_ [28], $V_{\mathbf{y}}:\mathbb{R}^{n}\rightarrow\mathbb{R}$, where
$n$ is the dimension of $\mathbf{y}$, such that
$\frac{\partial{V}_{\mathbf{y}}(\mathbf{y})}{\partial t}<0$ for all
$\mathbf{y}\in\mathbb{R}^{n}\setminus\\{F(\mathbf{x}^{*}_{1}),\ldots,F(\mathbf{x}^{*}_{n_{e}})\\}$,
and $\frac{\partial{V}_{\mathbf{y}}(\mathbf{y^{*}})}{\partial t}=0$ for
$\mathbf{y}^{*}\in\\{F(\mathbf{x}^{*}_{1}),\ldots,F(\mathbf{x}^{*}_{n_{e}})\\}$.
We assume the candidate function to be
$V_{\mathbf{y}}=V_{\mathbf{x}}(F^{-1}(\mathbf{y}))$, where $V_{\mathbf{x}}$ is
a valid Lyapunov function of the asymptotically stable
$\mathbf{x}^{\prime}(t)=g(\mathbf{x}(t))$, with
$\frac{\partial{V}_{\mathbf{x}}(\mathbf{x})}{\partial t}<0$ for
$\mathbf{x}\in\mathbb{R}^{b}\setminus\\{\mathbf{x}^{*}_{1},\ldots,\mathbf{x}^{*}_{n_{e}}\\}$
and $\frac{\partial{V}_{\mathbf{x}}(\mathbf{x}^{*})}{\partial t}=0$ for
$\mathbf{x}^{*}\in\\{\mathbf{x}^{*}_{1},\ldots,\mathbf{x}^{*}_{n_{e}}\\}$.
Consider the time derivative of the candidate function:

$\displaystyle\frac{\partial V_{\mathbf{y}}(\mathbf{y})}{\partial t}$
$\displaystyle=\frac{\partial
V_{\mathbf{y}}}{\partial\mathbf{y}}\frac{\partial\mathbf{y}}{\partial
t}=\frac{\partial V_{\mathbf{y}}}{\partial\mathbf{y}}f(\mathbf{y})$ (10)
$\displaystyle=\Big{(}\frac{\partial{V_{\mathbf{x}}}}{\partial\mathbf{x}}\frac{\partial
F^{-1}}{\partial\mathbf{y}}\frac{\partial
F}{\partial\mathbf{x}}g(\mathbf{x})\Big{)}_{\mathbf{x}=F^{-1}(\mathbf{y})}$
(11)
$\displaystyle=\Big{(}\frac{\partial{V_{\mathbf{x}}}}{\partial\mathbf{x}}J_{F}(\mathbf{x})^{-1}J_{F}(\mathbf{x})g(\mathbf{x})\Big{)}_{\mathbf{x}=F^{-1}(\mathbf{y})}$
By the inverse function theorem [12], (12)
$\displaystyle=\Big{(}\frac{\partial{V_{\mathbf{x}}}}{\partial\mathbf{x}}g(\mathbf{x})\Big{)}_{\mathbf{x}=F^{-1}(\mathbf{y})}=\Big{(}\frac{\partial
V_{\mathbf{x}}(\mathbf{x})}{\partial
t}\Big{)}_{\mathbf{x}=F^{-1}(\mathbf{y})}.$ (13)

Therefore, our candidate $V_{\mathbf{y}}$ is a valid Lyapunov function for
$\mathbf{y}^{\prime}(t)=f(\mathbf{y}(t))$. Thus, the system
$\mathbf{y}^{\prime}(t)=f(\mathbf{y}(t))$ is asymptotically stable. ∎

##  Appendix D Additional Implementation details

We run all of our experiments on a machine with an Intel i7-3770k 3.50GHz
processor, 32GB RAM and an NVIDIA GTX1080 GPU, with 8GB vRAM. For all of our
experiments, we use the optimiser ADAM with step-size $10^{-4}$, except for
the experiments in the Latent ODE, which where we use the standard set-up from
the Latent-ODE repository [4]. The dynamics models of compared ODEs have the
architecture: Input->dense(Input dimensions,
150)->tanh()->dense(150,150)->tanh()->dense(150,150)->tanh()->dense(150,150)->tanh()->dense(150,150)->tanh()->dense(150,output
dimensions)->output. Except for the Latent-ODE comparisons where settings from
the original repository [4] is used, and for the stiff system, we train for
500 iterations with step-size $10^{-4}$, and then train with step-size
$10^{-6}$ for 4500 iterations. For all of the experiments, except latent ODE
experiments where we follow the original set-up, we train for $5000$
iterations in total.

For all the experiments where we directly learn a dynamical system, we use an
invertible neural network with 5 invertible layers, and sub-networks with one
hidden layer of 1500 units. For non-linear base ODEs parameterised with a
simple neural networks, we use the architecture: Input->dense(Input
dimensions,30)->tanh()->dense(30,30)->tanh()->dense(30,30)->tanh()->dense(30,Output
dimensions)->Outputs. Additionally, all learned dynamics, both with ours and
compared methods, excepted when adhering to the original Latent-ODE set-up,
were augmented with the same number of additional zeros as original state
dimensions, for example 3 dimensional systems were augmented to 6 dimensions.

The Lotka-Volterra system used has the dynamics:

$\displaystyle x^{\prime}(t)$ $\displaystyle=x(t)(0.75-0.75y(t))$ (14)
$\displaystyle y^{\prime}(t)$ $\displaystyle=y(t)(-0.75+0.75x(t)-0.75z(t))$
(15) $\displaystyle z^{\prime}(t)$ $\displaystyle=z(t)(-0.75+0.75y(t))$ (16)

for $t\in[0,7]$ with initial conditions
$\\{(5,5,1),(2,6,6),(3,1,4),(7,1,2),(6,2,4),\\\
(3,3,1),(2,2,2),(4,4,3),(3,3,4),(1,1,5)\\}$.

The Lorenz system used has the dynamics:

$\displaystyle x^{\prime}(t)$ $\displaystyle=10(y(t)-x(t))$ (17)
$\displaystyle y^{\prime}(t)$ $\displaystyle=x(t)(28-y(t))-x(t)$ (18)
$\displaystyle z^{\prime}(t)$ $\displaystyle=x(t)y(t)-\frac{8}{3}z(t)$ (19)

for $t\in[0,2]$ with the initial conditions $(0.15,0.15,0.15)$.

The Robertson’s system used has the dynamics:

$\displaystyle x^{\prime}(t)$ $\displaystyle=-0.04x(t)+3\times 10^{4}y(t)z(t)$
(20) $\displaystyle y^{\prime}(t)$ $\displaystyle=0.04x(t)-3\times
10^{4}y(t)^{2}-10^{4}y(t)z(t)$ (21) $\displaystyle z^{\prime}(t)$
$\displaystyle=3\times 10^{4}y(t)^{2}$ (22)

for $t\in[0,120]$ with the initial conditions $(1,0,0)$.

In the latent-ODE problem setup, an observable time-series is assumed to have
latent variables which follow some ODE dynamics, and uses an Encoder -> ODE ->
Decoder architecture where an ODE is used to model the hidden latent
dimensions between the Encoder and Decoder. Note that a valid ODE is not
guaranteed in the space of observable data, but only in the latent dimensions.
Our set-up follows the repository given by [4], with the training settings for
the Encoder and Decoder architecture as below:

Periodic 100: We train the entire model for 500 epochs with Adamax optimiser
and an initial learning rate of $10^{-2}$. We sub-sample 5% of the original
time points and the size of the latent state is 10. The noise weight is set as
0.01 and the total number of time points is 100. For the Neural ODE
architectures, there is one layer in the recognition ODE and one layer in the
generative ODE, and 100 unit per layers. For the GRU unit there exists 100
units per layer for the GRU update network. All the above settings are exactly
the same as the configuration given in repository [4].

Periodic 1000: Settings are the same as _Periodic 100_ , except that the total
number of time points is set as 1000 to predict for finer time steps.

Human Activity: The model is trained for 200 epochs, with a dimensionality of
15 in the latent state. There are 4 layers in the recognition ODE and 2 layers
in the generative ODE, and 500 units per layer. The GRU unit has exists 50
units per layer. These settings are exactly the same as the configuration
given in the original repository [4].

ECG: Settings are the same as the classification task of _Human Activity_ ,
except that we use the ECG Heartbeat data available at [1].

##  Appendix E Additional Figures

We provide figures for learning an additional Lorenz system for $t\in[0,5]$,
with trajectory at initial condition $(-3.1,-1.15,8.15)$. We see that our
method, with a base ODE parameterised by a neural network, can generate
trajectories that closely match the ground truth:

Figure 6: A learned Lorenz system with the generated trajectory, at 10x data
resolution, and ground truth. Figure 7: A learned Lorenz system with the
generated trajectory, at 10x data resolution, and ground truth, rolled out in
time

We provide the change in coordinates over time, for the trajectory shown in
figure 5(a) in the paper:

Figure 8: The corresponding plot showing the coordinates over the time
interval of a learned Lorenz system over $t\in[0,2]$, at 10x data resolution,
which corresponds to the 3d figure shown as fig 5(a) in the main paper.

The following figures show how the diffeomorphism finds a mapping between a
stiff ODE with variations in very different scales across dimensions, to a
more manageable base ODE where the scales are more similar. The left plot is
in log-scale, while the right plot shows the same ODEs in linear scale:

Figure 9: The data is generated from a stiff ODE with large differences in
scales (indicated in crosses) over the different dimensions. We see that the
corresponding base ODE (indicated in solid lines) has changes with much more
similar scales.

We provide additional plots of trajectories, at different start points, from a
learned Lotka-Volterra system. The ground truth data is in blue, while
generated trajectory, of $10$x data resolution, is in red.

Figure 10: We see that trajectories from the learned Lotka-Volterra system, in
red, closely matches the ground truth, in blue.

We provide an additional figure for trajectories generated at unseen starting
points after being trained on the “imitation C” training data. The four
generated trajectories are in red, while the ground truths are in blue. Our
generated trajectories match the ground truth, and accurately capture the
motion of drawing a “C” character.

Figure 11: Robot motion trajectories in red, that imitate drawing a “C”
character. The ground truth is given in blue.

##  Appendix F Licenses for Packages

Common scientific packages used in our code include: (i) Numpy [20] (BSD
license), for general linear algebra and miscellaneous math operations (ii)
Matplotlib [22] (BSD compatible custom license), for plotting figures.

More specialised packages used include (i) FrEIA [3] (MIT license), for
invertible neural networks; (ii) TorchDiffEq [2] (MIT license), for
differentiable numerical integrators; (iii) Pytorch [34] (BSD license), for
optimisation and automatic differentiation; (iv) Latent-ODE [4] (MIT license),
for latent ODE implementation.


††thanks: These authors contribute equally to this work.††thanks: These
authors contribute equally to this work.

# Nonlinearity-induced dynamical self-organized twisted-bilayer lattices in
Bose-Einstein condensates

Rui Tian  Ministry of Education Key Laboratory for Nonequilibrium Synthesis
and Modulation of Condensed Matter,Shaanxi Province Key Laboratory of Quantum
Information and Quantum Optoelectronic Devices, School of Physics, Xi’an
Jiaotong University, Xi’an 710049, China  Yue Zhang  Ministry of Education Key
Laboratory for Nonequilibrium Synthesis and Modulation of Condensed
Matter,Shaanxi Province Key Laboratory of Quantum Information and Quantum
Optoelectronic Devices, School of Physics, Xi’an Jiaotong University, Xi’an
710049, China  Tianhao Wu  Ministry of Education Key Laboratory for
Nonequilibrium Synthesis and Modulation of Condensed Matter,Shaanxi Province
Key Laboratory of Quantum Information and Quantum Optoelectronic Devices,
School of Physics, Xi’an Jiaotong University, Xi’an 710049, China  Min Liu
Ministry of Education Key Laboratory for Nonequilibrium Synthesis and
Modulation of Condensed Matter,Shaanxi Province Key Laboratory of Quantum
Information and Quantum Optoelectronic Devices, School of Physics, Xi’an
Jiaotong University, Xi’an 710049, China  Yong-Chang Zhang  Ministry of
Education Key Laboratory for Nonequilibrium Synthesis and Modulation of
Condensed Matter,Shaanxi Province Key Laboratory of Quantum Information and
Quantum Optoelectronic Devices, School of Physics, Xi’an Jiaotong University,
Xi’an 710049, China  Shuai Li<EMAIL_ADDRESS>Ministry of
Education Key Laboratory for Nonequilibrium Synthesis and Modulation of
Condensed Matter,Shaanxi Province Key Laboratory of Quantum Information and
Quantum Optoelectronic Devices, School of Physics, Xi’an Jiaotong University,
Xi’an 710049, China  Bo Liu<EMAIL_ADDRESS>Ministry of Education Key
Laboratory for Nonequilibrium Synthesis and Modulation of Condensed
Matter,Shaanxi Province Key Laboratory of Quantum Information and Quantum
Optoelectronic Devices, School of Physics, Xi’an Jiaotong University, Xi’an
710049, China

###### Abstract

Creating crystal bilayers twisted with respect to each other would lead to
large periodic supercell structures, which can support a wide range of novel
electron correlated phenomena, where the full understanding is still under
debate. Here, we propose a new scheme to realize a nonlinearity-induced
dynamical self-organized twisted-bilayer lattice in an atomic Bose-Einstein
condensate (BEC). The key idea here is to utilize the nonlinear effect from
the intrinsic atomic interactions to couple different layers and induce a
dynamical self-organized supercell structure, dramatically distinct from the
conventional wisdom to achieve the static twisted-bilayer lattices. To
illustrate that, we study the dynamics of a two-component BEC and show that
the nonlinear interaction effect naturally emerged in the Gross-Pitaevskii
equation of interacting bosonic ultracold atoms can dynamically induce both
periodic (commensurable) and aperiodic (incommensurable) moiré structures. One
of the interesting moiré phenomena, i.e., the flat-band physics, is shown
through investigating the dynamics of the wave packet of BEC. Our proposal can
be implemented using available state-of-the-art experimental techniques and
reveal a profound connection between the nonlinearity and twistronics in cold
atom quantum simulators.

##  I Introduction

Superimposed crystal bilayers twisted with respect to each other provide a new
knob to unveil a wide range of novel electron correlated phenomena Dos Santos
et al. (2007); Sboychakov et al. (2015); Dai et al. (2016); Liu et al. (2020);
Cao et al. (2020); Chen et al. (2019); Shen et al. (2020); Regan et al.
(2020); Tang et al. (2020). In particular, moiré patterns associated with
small twist angles result in various novel physical properties Lu et al.
(2013); Woods et al. (2014); Yan et al. (2013); Alden et al. (2013); Kim et
al. (2016); Yin et al. (2016); Huder et al. (2018); Yankowitz et al. (2012);
Cao et al. (2018a); Yankowitz et al. (2019); Lu et al. (2019); Sharpe et al.
(2019); Seyler et al. (2019); Tran et al. (2019); Jin et al. (2019); Wang et
al. (2020a), including the flat band Bistritzer and MacDonald (2011); Cao et
al. (2018b), fractional Chern insulators Chen et al. (2020); Abouelkomsan et
al. (2020), quantum Hall effect Dean et al. (2013); Serlin et al. (2020);
Bultinck et al. (2020); Wu and Das Sarma (2020), moiré ferroelectricity Zheng
et al. (2020, 2020), exotic magnetism Sharpe et al. (2019); Seo et al. (2019);
Liu et al. (2019); Tschirhart et al. (2021), moiré excitons Seyler et al.
(2019); Alexeev et al. (2019); Andersen et al. (2021), non-Abelian gauge
potentials San-Jose et al. (2012); Ahn et al. (2019); Park et al. (2019);
Nuckolls et al. (2020), and unconventional superconductivity Cao et al.
(2018b, a); Yankowitz et al. (2019); Lu et al. (2019). In parallel to these
developments in solid state materials, examples of twisted-bilayer are also
studied in other artificial quantum systems, such as in optics Huang et al.
(2016); Wang et al. (2020b); Fu et al. (2020); Mao et al. (2021); Arkhipova et
al. (2023), acoustics Gardezi et al. (2021); Zheng et al. (2022); Jiang et al.
(2022); Duan et al. (2023) and cold atoms Soltan-Panahi et al. (2011);
Tarruell et al. (2012); Jo et al. (2012); Gall et al. (2021); Meng et al.
(2023). In photonic lattices, the moiré pattern and its induced phenomena,
such as the localization of light and flat band physics, have been studied in
the photorefractive medium Huang et al. (2016); Wang et al. (2020b); Fu et al.
(2020) and nanostructured photonic moiré superlattices Mao et al. (2021).

Ultracold atoms in optical lattices are another one of the most promising
platform to investigate the remarkable physical properties of twisted
artificial structures Soltan-Panahi et al. (2011); Tarruell et al. (2012); Jo
et al. (2012); Gall et al. (2021); Meng et al. (2023). It thus has stimulated
tremendous amount of efforts, including theoretically predicting moiré physics
in ultracold fermions Salamon et al. (2020a, b); Luo and Zhang (2021); Lee and
Pixley (2022) and bosons González-Tudela and Cirac (2019). In particular,
experimental breakthrough observation of the spatial moiré pattern and
momentum diffraction of a Bose-Einstein condensate (BEC) in twisted bilayer
optical lattices has recently been achieved Meng et al. (2023). However, the
conventional wisdom in past studies focuses on achieving the static twisted-
bilayer lattices Meng et al. (2023); Salamon et al. (2020a, b); Luo and Zhang
(2021); Lee and Pixley (2022); González-Tudela and Cirac (2019). How to
realize such twisted lattice structures in a dynamical way has so far eluded
efforts.

Figure 1:  Schematic illustration of the system. A two-component BEC is
considered here. Component 1 atoms are trapped in a tight harmonic trap and
component 2 atoms are loaded in a square optical lattice confined by a very
shallow harmonic trapping potential. After the system reaches equilibrium, the
tight harmonic trap for component 1 atoms is tuned to be very shallow and a
square optical lattice is also added. Distinguished from previous schemes via
adding external single-particle type interlayer coupling Meng et al. (2023);
Luo and Zhang (2021); González-Tudela and Cirac (2019), such as adding
microwave fields between two layers, our scheme only relies on the nonlinear
effect from the intrinsic atomic interactions.

In this work, we propose a new scheme to dynamically achieve twisted-bilayer
lattices that builds on the concept of nonlinearity. Our scheme relies on the
nonlinear interaction effect, which naturally emerges in the Gross-Pitaevskii
equations of interacting bosonic ultracold atoms. We find that such nonlinear
interaction effect amazingly results in a dynamical self-organized twisted-
bilayer lattice. This new scheme frees out the usual requirement of adding
external single-particle type coupling between different layers and only
relies on the nonlinear effect from the intrinsic atomic interactions,
dramatically distinguished from previous schemes Meng et al. (2023); Luo and
Zhang (2021); González-Tudela and Cirac (2019). We show that both dynamical
self-organized periodic (commensurable) and aperiodic (incommensurable)
twisted-bilayer lattice structures can be realized through our proposed
scheme. In particular, for the aperiodic case, we find the localization of the
Gaussian BEC wave packet, resulting from the emergence of the flat bands in
the nonlinearity-induced self-organized twisted-bilayer lattices. It indicates
that our proposed nonlinear interaction effect can dynamically produce the
tunable moiré lattices in BEC.

##  II Effective model

The system considered here is schematically illustrated in Fig.1. A two-
component BEC, for instance 87Rb atoms, is used, where a spin-dependent
harmonic potential is applied Arora et al. (2011); Bause et al. (2020);
Pegahan et al. (2019). Atoms in component 1 are trapped in a tight harmonic
trap and the atoms in component 2 are loaded in a square optical lattice
confined by a very shallow harmonic trapping potential. After the system
reaches equilibrium, the tight harmonic trap for component 1 atoms is tuned to
be very shallow and a square optical lattice is added. At the same time, the
inter-species interaction is turned on. Such a system can be described by the
following model Hamiltonian

$\displaystyle\mathbf{H}$ $\displaystyle=$
$\displaystyle\int\text{d}x\text{d}y\sum_{\alpha=1,2}\psi_{\alpha}^{\dagger}\left[\frac{\mathbf{k}^{2}}{2m}+V_{\text{trap},\alpha}(x,y)+V_{\alpha}(x,y)\right]\psi_{\alpha}$
(1) $\displaystyle+$
$\displaystyle\int\text{d}x\text{d}y(g_{11}\hat{n}_{1}^{2}+g_{22}\hat{n}_{2}^{2}+2g_{12}\hat{n}_{1}\hat{n}_{2}),$

where $\psi_{\alpha}$ $(\alpha=1,2)$ denote the Bose field operators for two-
component bosons, respectively. The corresponding density operators are
defined as $\hat{n}_{\alpha}=\psi^{\dagger}_{\alpha}\psi_{\alpha}$.
$V_{\text{trap},\alpha}$ denotes the harmonic trap. The spin-dependent optical
lattices are
$V_{1}(x,y)={V}_{1}[\sin^{2}(\cos\frac{\theta}{2}k_{L}x-\sin\frac{\theta}{2}k_{L}y)+\sin^{2}(\sin\frac{\theta}{2}k_{L}x+\cos\frac{\theta}{2}k_{L}y)]$
and
$V_{2}(x,y)={V}_{2}[\sin^{2}(\cos\frac{\theta}{2}k_{L}x+\sin\frac{\theta}{2}k_{L}y)+\sin^{2}(-\sin\frac{\theta}{2}k_{L}x+\cos\frac{\theta}{2}k_{L}y)]$,
respectively, with $\theta$ being the twist angle. $k_{L}=2\pi/\lambda$ is the
wavevector of the laser field with the corresponding lattice constant defined
as $a_{L}=\pi/k_{L}$, and $\lambda$ is the wavelength of the laser.
${V}_{\alpha}$ captures the corresponding lattice depth. $g_{11}$ and $g_{22}$
characterize the intra-species interaction strength for component 1 and
component 2 atoms, respectively. The inter-species interaction strength is
captured by $g_{12}$.

From Eq. (1), the following dimensionless energy functional can be constructed
under the Gross-Pitaevskii (GP) mean-field theory

$\displaystyle\varepsilon$ $\displaystyle=$
$\displaystyle\int\text{d}x^{\prime}\text{d}y^{\prime}\sum_{\alpha=1,2}\psi_{\alpha}^{*}\left[-\frac{\nabla^{2}}{2}+V^{\prime}_{\text{trap},\alpha}(x^{\prime},y^{\prime})+V^{\prime}_{\alpha}(x^{\prime},y^{\prime})\right]$
$\displaystyle\times$
$\displaystyle\psi_{\alpha}+\int\text{d}x^{\prime}\text{d}y^{\prime}(g^{\prime}_{11}|\psi_{1}|^{4}+g^{\prime}_{22}|\psi_{2}|^{4}+2g^{\prime}_{12}|\psi_{1}|^{2}|\psi_{2}|^{2}),$

where dimensionless spatial coordinates are defined as $x^{\prime}=x/\lambda$
and $y^{\prime}=y/\lambda$, respectively.
$V^{\prime}_{\text{trap},\alpha}=V_{\text{trap},\alpha}/(\hbar^{2}/m\lambda^{2})$
represents the harmonic trapping potential.
$V^{\prime}_{\alpha}=V_{\alpha}/(\hbar^{2}/m\lambda^{2})$ labels the
dimensionless lattice depth.
$g^{\prime}_{11(12,22)}=g_{11(12,22)}/(\hbar^{2}/m)$ are the corresponding
dimensionless interaction strength. The wavefunction is normalized as
$\iint|\psi_{\alpha}|^{2}\text{d}x^{\prime}\text{d}y^{\prime}=N_{\alpha}$ with
$N_{\alpha}$ being the particle number of the two-component atoms,
respectively. Based on Eq. (2), the dynamics of the system can be well
described by the following coupled GP equations

$\displaystyle i\frac{\partial\psi_{1}}{\partial t^{\prime}}=$
$\displaystyle\Big{[}-\frac{1}{2}\nabla^{2}+V^{\prime}_{\text{trap,1}}(x^{\prime},y^{\prime})+V^{\prime}_{1}(x^{\prime},y^{\prime})+g^{\prime}_{11}|\psi_{1}|^{2}$
$\displaystyle+g^{\prime}_{12}|\psi_{2}|^{2}\Big{]}\psi_{1}$ (3a)
$\displaystyle i\frac{\partial\psi_{2}}{\partial t^{\prime}}=$
$\displaystyle\Big{[}-\frac{1}{2}\nabla^{2}+V^{\prime}_{\text{trap,2}}(x^{\prime},y^{\prime})+V^{\prime}_{2}(x^{\prime},y^{\prime})+g^{\prime}_{12}|\psi_{1}|^{2}$
$\displaystyle+g^{\prime}_{22}|\psi_{2}|^{2}\Big{]}\psi_{2},$ (3b)

where $t^{\prime}=t/(m\lambda^{2}/\hbar)$ denotes the dimensionless evolution
time.

Figure 2: Time evolution of the two-component BEC in the nonlinearity-induced
dynamical self-organized twisted-bilayer lattices. The left column of (b) and
(c) show the distinct dynamics of the component 1 BEC with different twist
angles $\theta$, starting from a Gaussian wave packet as the initial state
shown in the left column of (a). In the presence of nonlinear term $w$
(defined in the main text), even without external single-particle type
coupling between two layers, one of the interesting moiré physics, i.e., both
the delocalized and localized behavior of the Gaussion wave packet of
component 1 atoms, can be achieved when $\theta$ being chosen as a Pythagorean
angle, $\theta=\arctan 3/4$ in (b) and a non-Pythagorean angle $\theta=\pi/6$
in (c), respectively. The right column shows the time evolution of the
component 2 BEC. Here we choose $g^{\prime}_{11}=1.85\times 10^{-4}$,
$g^{\prime}_{22}=2\times 10^{-2}$, $g^{\prime}_{12}=2.1\times 10^{-3}$,
$V^{\prime}_{1}=36$, $V^{\prime}_{2}=333$, $N_{1}=1\times 10^{3}$ and
$N_{2}=5\times 10^{5}$. In (b) and (c), the evolution time is chosen as
$t^{\prime}=10$.

##  III Moiré localization induced by the nonlinear interaction effect

In the following, we will demonstrate that distinguished from previous studies
Meng et al. (2023); Luo and Zhang (2021); González-Tudela and Cirac (2019),
even without external single-particle type interlayer coupling, the existence
of nonlinear interaction term in GP equation, i.e.,
$w=g^{\prime}_{11}|\psi_{1}|^{2}+g^{\prime}_{12}|\psi_{2}|^{2}$, can
dynamically result in the twisted-bilayer lattices and lead to one of the
interesting moiré phenomena. To show that, we monitor the dynamics of BEC for
different twist angles $\theta$. Here the initial state for component 1 atoms
is the ground state of the harmonic trap, which is chosen as a Gaussian wave
packet (left part in Fig. 2(a)). While for component 2 atoms, the initial
state is chosen as the ground state (right part in Fig.2(a)) determined by the
effective Hamiltonian in the right hand-side of Eq. (3b) with
$g^{\prime}_{12}=0$.

Interestingly, as shown in Fig. 2, when tuning the twist angle $\theta$, in
the presence of nonlinear term $w$, even without external single-particle type
coupling between two layers, one of the interesting moiré physics, i.e., both
the localized and delocalized behavior of component 1 atoms, can be achieved.
For instance, as shown in Fig. 2(b) and (c), when $\theta$ is chosen as a
Pythagorean angle, for instance, $\theta=\arctan 3/4$ (Fig. 2(b)), the
Gaussion wave packet for component 1 atoms displays the delocalized behavior.
While considering $\theta$ being the non-Pythagorean twisted angle, for
instance, $\theta=\pi/6$ (Fig. 2(c)), the wave packet for component 1 atoms
turns out to be localized. Note that here we consider the case with large
number imbalance between two components, i.e., $N_{1}\ll N_{2}$ and lager
intra-species interaction, i.e., $g^{\prime}_{22}>g^{\prime}_{12}$. The effect
of inter-species interaction on the dynamics of component 2 atoms is highly
suppressed and the density profile of component 2 atoms almost does not change
during the evolution, as shown in the right column of Fig. 2.

Figure 3:  Mean square displacement of the component 1 BEC wave packet as the
function of the evolution time. When the twist angle $\theta$ being a
Pythagorean angle, i.e., $\theta=\arctan 3/4$, $\sigma(t)$ tends to linearly
increase with respect to the evolution time, i.e., $\sigma(t)\propto t$,
indicating its delocalized behavior. While $\theta$ being a non-Pythagorean
angle, i.e., $\theta=\pi/6$, $\sigma(t)$ reaches a stable value in a fairly
short time, characterizing its localized behavior. Other parameters are chosen
as the same in Fig. 2.

To quantitatively distinguish between the delocalized and localized behavior
of the wave packet of component 1 atoms, we analyze the following mean square
displacement (MSD) Hiramoto and Abe (1988); Zhong et al. (2001)

$\sigma(t)=\left[\frac{\iint\text{d}x\text{d}y(x^{2}+y^{2})|\psi_{1}(t)|^{2}}{\iint\text{d}x\text{d}y|\psi_{1}(t)|^{2}}\right]^{\frac{1}{2}},$
(4)

which can characterize the width of the wave packet of component 1 atoms. The
different dynamical behaviors of the wave packet can be quantified as
$\sigma(t)\propto t^{\kappa}$ with distinct dynamical index $\kappa$. As shown
in Fig. 3, when the twist angle $\theta$ is chosen as a Pythagorean angle, for
instance, $\theta=\arctan 3/4$, $\sigma(t)$ tends to linearly increase with
respect to the evolution time, i.e., $\sigma(t)\propto t$. The dynamical index
is 1, indicating its delocalized behavior. While $\theta$ being a non-
Pythagorean angle, for instance, $\theta=\pi/6$, $\sigma(t)$ will reach a
stable value in a fairly short time, indicating the dynamical index $\kappa$
being 0. It thus characterizes the localized behavior of the wave packet.
Therefore, through tuning the twist angle $\theta$, one of the moiré physics,
i.e., both the localized and delocalized behavior of the wave packet of
component 1 atoms, can be achieved.

##  IV Nonlinearity-induced dynamical self-organized twisted-bilayer lattices

In the following, we will explain the reason for the appearance of the moiré
physics in our proposed system. Let us focus on Eq. (3a) in the GP equation
describing the two-component atomic BEC. In the presence of the nonlinear term
$w$, component 1 atoms suffer from a dynamical self-organized effective
potential $V_{eff}\equiv
V^{\prime}_{\text{trap,1}}+V^{\prime}_{1}+w=V^{\prime}_{\text{trap,1}}+V^{\prime}_{1}+g^{\prime}_{11}|\psi_{1}|^{2}+g^{\prime}_{12}|\psi_{2}|^{2}$.
As shown in Fig. 4(a), when the twist angle $\theta$ is chosen as a
Pythagorean angle, for instance, $\theta=\arctan 3/4$, the effective potential
$V_{eff}$ possess a spatially periodic structure. While considering $\theta$
being a non-Pythagorean twist angle, for instance, $\theta=\pi/6$, as shown in
Fig. 4(b), $V_{eff}$ behaves as a twisted-bilayer lattice with a spatial
aperiodic structure. Therefore, distinguished from previous studies Meng et
al. (2023); Luo and Zhang (2021); González-Tudela and Cirac (2019), even
without external single-particle type interlayer coupling, the nonlinear
effect resulting from the intrinsic atomic interactions provides a new way to
dynamically induce the twisted-bilayer lattices. Through tuning the twist
angle $\theta$, we find that both periodic (commensurable) and aperiodic
(incommensurable) moiré structures can be created in such a dynamical twisted-
bilayer lattice.

Figure 4: Nonlinearity-induced self-organized twisted-bilayer lattice
$V_{eff}$ for component 1 BEC. Both periodic (commensurable) and aperiodic
(incommensurable) structures can be created in such a nonlinearity-induced
twisted-bilayer lattice through adjusting the twist angle $\theta$ from the
Pythagorean angle to the non-Pythagorean angle , where in (a) $\theta=\arctan
3/4$ and in (b) $\theta=\pi/6$, respectively. Here, the evolution time is
chosen as $t^{\prime}=10$. Other parameters are chosen as the same in Fig. 2.

To further understand the distinct dynamical behaviors of the wave packet of
component 1 atoms, we rely on the analysis of the appearance of the flat band
physics in our proposed nonlinearity-induced self-organized twisted-bilayer
lattices. To do that, we calculate the band structure of $V_{eff}$ for both
periodic (commensurable) and aperiodic (incommensurable) moiré structures. In
particular, when considering the aperiodic moiré lattice, the band structure
can be calculated through approximating the non-Pythagorean twist angle by a
Pythagorean one. For instance, here we use $\theta=\arctan 120/209$ to
approximate $\theta=\pi/6$. The band structure of $V_{eff}$ can thus be
obtained by the plane-wave expansion method through introducing the Bloch
basis
$\psi_{n\vec{k}}=\sum_{\vec{G}}u_{n\vec{k},\vec{G}}|\vec{k}+\vec{G}\rangle$
with the Bloch vector $\vec{k}$ and the reciprocal lattice vector $\vec{G}$.
Here $n$ labels the band index. Then, the energy spectra can be obtained
through solving the eigen-problem via the following relation

$\displaystyle\frac{(\vec{k}+\vec{G})^{2}\lambda^{2}}{2}u_{n\vec{k},\vec{G}}+\sum_{\vec{G^{\prime}}}\langle\vec{k}+\vec{G}|V_{eff}|\vec{k}+\vec{G^{\prime}}\rangle
u_{n\vec{k},\vec{G^{\prime}}}$ (5)
$\displaystyle=E^{\prime}_{n\vec{k}}u_{n\vec{k},\vec{G}},$

where $E^{\prime}_{n\vec{k}}=E_{n\vec{k}}/(\hbar^{2}/m\lambda^{2})$ is the
dimensionless energy spectra.

Figure 5:  The band structure of the nonlinearity-induced self-organized
twisted-bilayer lattice. (a) and (b) show the corresponding band structure
(the lowest five bands) of the moiré lattices $V_{eff}$ demonstrated in Fig. 4
(a) and (b), respectively. It is shown that the flat bands emerge when
changing the twist angle $\theta$ from the Pythagorean angle (in (a)) to the
non-Pythagorean angle (in (b)), accompanying with the lattice structure
changing from the periodic case to the aperiodic case. Other parameters are
chosen as the same in Fig. 4.

As shown in Fig. 5, the band structures of the corresponding twisted-bilayer
lattices $V_{eff}$ demonstrated in Fig. 4 are obtained. We find that when
changing the twist angle from a Pythagorean angle to a non-Pythagorean angle,
more dispersive lower bands become extremely flat, accompanying with the
lattice structure changing from the periodic (commensurable) case to the
aperiodic (incommensurable) one. In particular, the appearance of the lowest
flat band will lead to the localization of the wave packet of component 1
atoms, since such a BEC wave packet will mainly concentrate at the lowest band
after a long enough time evolution. Furthermore, it is known that the flat
band can support quasi-nondiffracting localized modes. Therefore, the wave
packet of component 1 atoms loaded in such a nonlinearity induced self-
organized aperiodic (incommensurable) twisted-bilayer lattice will tend to be
localized.

##  V Conclusion

In summary, we have demonstrated a new approach to dynamically achieve
twisted-bilayer lattices in an atomic BEC. To emphasize a remarkable
difference, our scheme does not require external single-particle type
interlayer coupling. Instead, we only require the nonlinear effect arising
from the intrinsic atomic interactions. We further show that both periodic and
aperiodic moire lattices, i.e., nonlinearity-induced dynamical self-organized
twisted-bilayer lattices, can be designed through tuning the twist angle. One
of the interesting phenomena in twisted-bilayer lattices, i.e., the flat-band
physics, is also shown through investigating the dynamics of the wave packet
of BEC. Such nonlinearity-induced approach is rather generic to other quantum
systems, such as in photonics and electrical circuits, than restricted to the
setup considered in this work, since the nonlinear effects, such as Kerr
nonlinearity, is common in these systems. Our results reveal a profound
connection between the nonlinearity and twistronics.

##  VI Acknowledgment

We acknowledge helpful discussions with L. Chen. This work is supported by the
National Key R$\&$D Program of China (2021YFA1401700), NSFC (Grants No.
12074305, 12147137, 11774282), the Fundamental Research Funds for the Central
Universities (Grant No. xtr052023002), the Shaanxi Fundamental Science
Research Project for Mathematics and Physics (Grant No. 23JSZ003) and Xiaomi
Young Scholar Program. We also thank the HPC platform of Xi’an Jiaotong
University, where our numerical calculations was performed.

## References

  * Dos Santos et al. (2007) J. L. Dos Santos, N. Peres, and A. C. Neto, Phys. Rev. Lett. 99, 256802 (2007).
  * Sboychakov et al. (2015) A. Sboychakov, A. Rakhmanov, A. Rozhkov, and F. Nori, Phys. Rev. B 92, 075402 (2015).
  * Dai et al. (2016) S. Dai, Y. Xiang, and D. J. Srolovitz, Nano Lett. 16, 5923 (2016).
  * Liu et al. (2020) X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, Y. Ronen, H. Yoo, D. Haei Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, et al., Nature 583, 221 (2020).
  * Cao et al. (2020) Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 583, 215 (2020).
  * Chen et al. (2019) G. Chen, A. L. Sharpe, P. Gallagher, I. T. Rosen, E. J. Fox, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, et al., Nature 572, 215 (2019).
  * Shen et al. (2020) C. Shen, Y. Chu, Q. Wu, N. Li, S. Wang, Y. Zhao, J. Tang, J. Liu, J. Tian, K. Watanabe, et al., Nat. Phys. 16, 520 (2020).
  * Regan et al. (2020) E. C. Regan, D. Wang, C. Jin, M. I. Bakti Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, et al., Nature 579, 359 (2020).
  * Tang et al. (2020) Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, et al., Nature 579, 353 (2020).
  * Lu et al. (2013) C.-C. Lu, Y.-C. Lin, Z. Liu, C.-H. Yeh, K. Suenaga, and P.-W. Chiu, ACS Nano 7, 2587 (2013).
  * Woods et al. (2014) C. Woods, L. Britnell, A. Eckmann, R. Ma, J. Lu, H. Guo, X. Lin, G. Yu, Y. Cao, R. V. Gorbachev, et al., Nat. Phys. 10, 451 (2014).
  * Yan et al. (2013) W. Yan, W.-Y. He, Z.-D. Chu, M. Liu, L. Meng, R.-F. Dou, Y. Zhang, Z. Liu, J.-C. Nie, and L. He, Nat. Commun. 4, 2159 (2013).
  * Alden et al. (2013) J. S. Alden, A. W. Tsen, P. Y. Huang, R. Hovden, L. Brown, J. Park, D. A. Muller, and P. L. McEuen, Proc. Natl. Acad. Sci. U. S. A. 110, 11256 (2013).
  * Kim et al. (2016) K. Kim, M. Yankowitz, B. Fallahazad, S. Kang, H. C. Movva, S. Huang, S. Larentis, C. M. Corbet, T. Taniguchi, K. Watanabe, et al., Nano Lett. 16, 1989 (2016).
  * Yin et al. (2016) L.-J. Yin, H. Jiang, J.-B. Qiao, and L. He, Nat. Commun. 7, 11760 (2016).
  * Huder et al. (2018) L. Huder, A. Artaud, T. Le Quang, G. T. De Laissardière, A. G. Jansen, G. Lapertot, C. Chapelier, and V. T. Renard, Phys. Rev. Lett. 120, 156405 (2018).
  * Yankowitz et al. (2012) M. Yankowitz, J. Xue, D. Cormode, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero, P. Jacquod, and B. J. LeRoy, Nat. Phys. 8, 382 (2012).
  * Cao et al. (2018a) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018a).
  * Yankowitz et al. (2019) M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019).
  * Lu et al. (2019) X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, et al., Nature 574, 653 (2019).
  * Sharpe et al. (2019) A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019).
  * Seyler et al. (2019) K. L. Seyler, P. Rivera, H. Yu, N. P. Wilson, E. L. Ray, D. G. Mandrus, J. Yan, W. Yao, and X. Xu, Nature 567, 66 (2019).
  * Tran et al. (2019) K. Tran, G. Moody, F. Wu, X. Lu, J. Choi, K. Kim, A. Rai, D. A. Sanchez, J. Quan, A. Singh, et al., Nature 567, 71 (2019).
  * Jin et al. (2019) C. Jin, E. C. Regan, A. Yan, M. Iqbal Bakti Utama, D. Wang, S. Zhao, Y. Qin, S. Yang, Z. Zheng, S. Shi, et al., Nature 567, 76 (2019).
  * Wang et al. (2020a) L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, et al., Nat. Mater. 19, 861 (2020a).
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Proc. Natl. Acad. Sci. U. S. A. 108, 12233 (2011).
  * Cao et al. (2018b) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, et al., Nature 556, 80 (2018b).
  * Chen et al. (2020) G. Chen, A. L. Sharpe, E. J. Fox, Y.-H. Zhang, S. Wang, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, et al., Nature 579, 56 (2020).
  * Abouelkomsan et al. (2020) A. Abouelkomsan, Z. Liu, and E. J. Bergholtz, Phys. Rev. Lett. 124, 106803 (2020).
  * Dean et al. (2013) C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, et al., Nature 497, 598 (2013).
  * Serlin et al. (2020) M. Serlin, C. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. Young, Science 367, 900 (2020).
  * Bultinck et al. (2020) N. Bultinck, S. Chatterjee, and M. P. Zaletel, Phys. Rev. Lett. 124, 166601 (2020).
  * Wu and Das Sarma (2020) F. Wu and S. Das Sarma, Phys. Rev. Lett. 124, 046403 (2020).
  * Zheng et al. (2020) Z. Zheng, Q. Ma, Z. Bi, S. de La Barrera, M.-H. Liu, N. Mao, Y. Zhang, N. Kiper, K. Watanabe, T. Taniguchi, et al., Nature 588, 71 (2020).
  * Seo et al. (2019) K. Seo, V. N. Kotov, and B. Uchoa, Phys. Rev. Lett. 122, 246402 (2019).
  * Liu et al. (2019) J. Liu, Z. Ma, J. Gao, and X. Dai, Phys. Rev. X 9, 031021 (2019).
  * Tschirhart et al. (2021) C. Tschirhart, M. Serlin, H. Polshyn, A. Shragai, Z. Xia, J. Zhu, Y. Zhang, K. Watanabe, T. Taniguchi, M. Huber, et al., Science 372, 1323 (2021).
  * Alexeev et al. (2019) E. M. Alexeev, D. A. Ruiz-Tijerina, M. Danovich, M. J. Hamer, D. J. Terry, P. K. Nayak, S. Ahn, S. Pak, J. Lee, J. I. Sohn, et al., Nature 567, 81 (2019).
  * Andersen et al. (2021) T. I. Andersen, G. Scuri, A. Sushko, K. De Greve, J. Sung, Y. Zhou, D. S. Wild, R. J. Gelly, H. Heo, D. Bérubé, et al., Nat. Mater. 20, 480 (2021).
  * San-Jose et al. (2012) P. San-Jose, J. González, and F. Guinea, Phys. Rev. Lett. 108, 216802 (2012).
  * Ahn et al. (2019) J. Ahn, S. Park, and B.-J. Yang, Phys. Rev. X 9, 021013 (2019).
  * Park et al. (2019) M. J. Park, Y. Kim, G. Y. Cho, and S. Lee, Phys. Rev. Lett. 123, 216803 (2019).
  * Nuckolls et al. (2020) K. P. Nuckolls, M. Oh, D. Wong, B. Lian, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 588, 610 (2020).
  * Huang et al. (2016) C. Huang, F. Ye, X. Chen, Y. V. Kartashov, V. V. Konotop, and L. Torner, Sci. Rep. 6, 32546 (2016).
  * Wang et al. (2020b) P. Wang, Y. Zheng, X. Chen, C. Huang, Y. V. Kartashov, L. Torner, V. V. Konotop, and F. Ye, Nature 577, 42 (2020b).
  * Fu et al. (2020) Q. Fu, P. Wang, C. Huang, Y. V. Kartashov, L. Torner, V. V. Konotop, and F. Ye, Nat. Photonics 14, 663 (2020).
  * Mao et al. (2021) X.-R. Mao, Z.-K. Shao, H.-Y. Luan, S.-L. Wang, and R.-M. Ma, Nat. Nanotechnol. 16, 1099 (2021).
  * Arkhipova et al. (2023) A. Arkhipova, Y. V. Kartashov, S. K. Ivanov, S. Zhuravitskii, N. Skryabin, I. Dyakonov, A. Kalinkin, S. Kulik, V. O. Kompanets, S. Chekalin, et al., Phys. Rev. Lett. 130, 083801 (2023).
  * Gardezi et al. (2021) S. M. Gardezi, H. Pirie, S. Carr, W. Dorrell, and J. E. Hoffman, 2D Mater. 8, 031002 (2021).
  * Zheng et al. (2022) S. Zheng, J. Zhang, G. Duan, Z. Jiang, X. Man, D. Yu, and B. Xia, Phys. Rev. B 105, 184104 (2022).
  * Jiang et al. (2022) Z. Jiang, J. Liu, S. Zheng, G. Duan, and B. Xia, Appl. Phys. Lett. 121 (2022).
  * Duan et al. (2023) G. Duan, S. Zheng, J. Zhang, Z. Jiang, X. Man, D. Yu, and B. Xia, Appl. Phys. Lett. 123 (2023).
  * Soltan-Panahi et al. (2011) P. Soltan-Panahi, J. Struck, P. Hauke, A. Bick, W. Plenkers, G. Meineke, C. Becker, P. Windpassinger, M. Lewenstein, and K. Sengstock, Nat. Phys. 7, 434 (2011).
  * Tarruell et al. (2012) L. Tarruell, D. Greif, T. Uehlinger, G. Jotzu, and T. Esslinger, Nature 483, 302 (2012).
  * Jo et al. (2012) G.-B. Jo, J. Guzman, C. K. Thomas, P. Hosur, A. Vishwanath, and D. M. Stamper-Kurn, Phys. Rev. Lett. 108, 045305 (2012).
  * Gall et al. (2021) M. Gall, N. Wurz, J. Samland, C. F. Chan, and M. Köhl, Nature 589, 40 (2021).
  * Meng et al. (2023) Z. Meng, L. Wang, W. Han, F. Liu, K. Wen, C. Gao, P. Wang, C. Chin, and J. Zhang, Nature 615, 231 (2023).
  * Salamon et al. (2020a) T. Salamon, R. W. Chhajlany, A. Dauphin, M. Lewenstein, and D. Rakshit, Phys. Rev. B 102, 235126 (2020a).
  * Salamon et al. (2020b) T. Salamon, A. Celi, R. W. Chhajlany, I. Frérot, M. Lewenstein, L. Tarruell, and D. Rakshit, Phys. Rev. Lett. 125, 030504 (2020b).
  * Luo and Zhang (2021) X.-W. Luo and C. Zhang, Phys. Rev. Lett. 126, 103201 (2021).
  * Lee and Pixley (2022) J. Lee and J. Pixley, SciPost Phys. 13, 033 (2022).
  * González-Tudela and Cirac (2019) A. González-Tudela and J. I. Cirac, Phys. Rev. A 100, 053604 (2019).
  * Arora et al. (2011) B. Arora, M. S. Safronova, and C. W. Clark, Phys. Rev. A 84, 043401 (2011).
  * Bause et al. (2020) R. Bause, M. Li, A. Schindewolf, X.-Y. Chen, M. Duda, S. Kotochigova, I. Bloch, and X.-Y. Luo, Phys. Rev. Lett. 125, 023201 (2020).
  * Pegahan et al. (2019) S. Pegahan, J. Kangara, I. Arakelyan, and J. E. Thomas, Phys. Rev. A 99, 063620 (2019).
  * Hiramoto and Abe (1988) H. Hiramoto and S. Abe, Journal of the Physical Society of Japan 57, 1365 (1988).
  * Zhong et al. (2001) J. Zhong, R. B. Diener, D. A. Steck, W. H. Oskay, M. G. Raizen, E. W. Plummer, Z. Zhang, and Q. Niu, Phys. Rev. Lett. 86, 2485 (2001).