# Defect - related Anomalous Mobility of Small Polarons in Oxides: the Case of
Congruent Lithium Niobate

Anton Pfannstiel1, Mirco Imlau1, Marco Bazzan2 and Laura Vittadello1 1
Institute of Physics, Department of Mathematics/Informatics/Physics,
University of Osnabrück, Barbarastraße 7, D-49076 Osnabrück, Germany  2
Università di Padova, Dipartimento di Fisica e Astronomia, Via Marzolo 8,
35131, Padova, Italy<EMAIL_ADDRESS>
###### Abstract

Polarons play a major role in the description of optical, electrical and
dielectrical properties of several ferroelectric oxides. The motion of those
particles occur by elementary hops among the material lattice sites. In order
to compute macroscopic transport parameters such as charge mobility, normal
(i.e. Fickian) diffusion laws are generally assumed. In this paper we show
that when defect states able to trap the polarons for long times are
considered, significant deviations from the normal diffusion behaviour arise.
As an example of this behavior, we consider here the case of lithium niobate
(LN). This can be considered as a prototypical system, having a rich landscape
of interacting polaron types and for which a significant wealth of information
is available in literature. Our analysis considers the case of a
stoichiometric, defect-free lithium niobate containing a certain concentration
of small electron polarons hopping on regular Nb sites, and compares it to the
material in congruent composition, which is generally found in real-life
applications and which is characterized by a large concentration of antisite
$\mathrm{Nb_{Li}}$ defects. While in the first case the charge carriers are
free polarons hopping on a regular Nb sublattice, in the second case a
fraction of polarons is trapped on antisite defects. Thus, in the congruent
material, a range of different hopping possibilities arises, depending on the
type of starting and destination sites. We develop a formalism encompassing
all these microscopic processes in the framework of a switching diffusion
model which can be well approximated by a mobile-immobile transport model
providing explicit expressions for the polaron mobility. Finally, starting
from the Marcus-Holstein’s model for the polaron hopping frequency we verify
by means of a Monte Carlo approach the diffusion/mobility of the different
polarons species showing that, while free polarons obeys the laws for normal
diffusion as expected, bound polarons follow an anomalous diffusion behaviour
and that in the case of the congruent crystal where mixed free and bound
polaron transport is involved, our expressions indeed provide a satisfactory
description.

  * April 2024

##  1 Introduction

In the last decades it became apparent that the charge transport in a number
of polar materials must be understood in terms of polaron motion. These are
quasi-particles made up of an electrical charge that, by interaction with the
crystalline environment, is able to distort the local lattice creating a local
potential well: as a net result, the particle becomes self-localized [1].
Polarons play a major role when it comes to the interpretation of the optical,
electrical or dielectrical properties [2] of important technological materials
such as $\mathrm{LiNbO_{3}}$ [3, 4], $\mathrm{KNbO_{3}}$ [5, 6],
$\mathrm{LiTaO_{3}}$ [7, 8], $\mathrm{BaTiO_{3}}$ [9, 10],
$\mathrm{PbTiO_{3}}$ [11], $\mathrm{r-TiO_{2}}$ [12, 13], $\mathrm{CeO_{2}}$
[14, 15] and solid solutions of the above like
$\mathrm{LiNb_{1-x}Ta_{x}O_{3}}$ [16] or $\mathrm{KTa_{1-x}Nb_{x}O_{3}}$ [17].
In these polar materials, the electron-phonon coupling is so strong that the
exceeding charge is localized on a single lattice site, i.e. only _small_
strong - coupling polarons are formed. These polarons can move in the lattice
via random hopping mechanisms in response to the lattice thermal motion, as
described by the Marcus-Holstein theory [18]. In standard formulations used to
compute the drift mobility of the polaronic carriers, it is assumed that the
latter hop on a regular lattice of spacing $d$ among by nearest-neighboring
sites. Thus the mobility of the polarons follows the familiar law $\mu\propto
ed^{2}w/(k_{\mathrm{B}}T)$ [19, 20].

However, the Marcus- Holstein model forecasts an exponential distance-
dependent hopping frequency: depending on the characteristic hopping length,
the hopping process may involve also far away hopping sites at a distance
bigger than $d$. Even more daunting is the fact that real polar materials are
often characterized by the presence in the lattice of a considerable amount of
point defects. Due to extra Coulomb interaction, they can act as privileged
localization centers for the polarons. When sitting on those sites, the
hopping frequency at a given temperature may differ of several orders of
magnitude with respect to free polarons moving on regular sites. Thus
different types of polarons may exist in a given material and their interplay
can result in a rich and complex behaviour, like the formation of polarons
bound to defects (bound polarons), polaron complexes made up of a regular and
a bound polaron (bipolarons) [21, 22] or a hole polaron and an electron
polaron (self-trapped excitons) [23]. In recent works [24, 4, 25, 26] it has
been shown that the relative weight of the different microscopic hopping
processes is strongly dependent on the temperature and on the defect
concentration.

The goal of the present work is to provide a novel formulation for the
calculation of the polaron mobility taking into account the complications
induced by the presence of defects and of the different types of polarons
arising from them. As a case study, we will consider $\mathrm{LiNbO_{3}}$,
(LN). This is a technologically important material widely used in a great
variety of applications ranging from electro-optic modulators, waveguides
[27], nonlinear optical devices, photorefractive holography and, more
recently, for novel applications such as the manipulation of small organic
[28] and inorganic objects [29], as substrate for photorefractive tweezers
[30] and as a platform for microfluidic chips [31, 32, 33, 34]. The structure
of stoichiometric LN is constituted by a quasi-cubic Nb sublattice
intercalated by Li and O atoms. In an ideal defect-free material
(stoichiometric LN, sLN) polarons can only form at the regular Nb sites
$\mathrm{Nb}_{\mathrm{Nb}}^{5+}$ ions and are indicated as _F_ polarons.
However, LN is normally produced at the congruent composition (cLN) which is
characterized by a high concentration (typically, up to 1 mol.%) of
substitutional “antisite” defects $\mathrm{Nb_{Li}^{5+}}$ [35]. These defects
constitute the preferential sites for the formation of polarons indicated as
bound polarons (_B_ polarons). While further polaron types and localization
sites exist, they contribute to a lesser extent to the transport of electrons
and are thus omitted in this work. In LN the relevant microscopic parameters
necessary to compute the Marcus-Holstein hopping frequency are known from
previous studies [3, 25, 4], making it an excellent test case. Our approach
can anyway be adapted to different materials.

##  2 Theory

Polarons move through the crystal via hopping transport. According to the
Marcus-Holstein theory, each transition is attributed a hopping frequency of
the following type: [25, 18]:

$w_{if}(r,T)=w_{if}^{0}(T)\exp\left(-\frac{r_{if}}{a_{if}}\right)$ (1a)
$w_{if}^{0}(T)=\frac{1}{2}\left(\frac{\pi}{k_{\mathrm{B}}T\lambda_{if}}\right)^{\frac{1}{2}}\frac{I_{if}^{2}}{\hbar}\exp\left(-\frac{U_{if}}{k_{\mathrm{B}}T}\right)$
(1b)
$U_{if}=\frac{(2E_{i}-\varepsilon_{i}+\varepsilon_{f}-e(\textbf{r}_{i}-\textbf{r}_{f})\cdot\textbf{F})^{2}}{4(E_{i}+E_{f})}$
(1c)

where $k_{\mathrm{B}}$ is the Boltzmann constant, $T$ the absolute
temperature, $\hbar$ the reduced Planck constant and $e$ the elementary
charge. The lattice reorganization energy $\lambda_{if}=E_{i}+E_{f}$ is spent
to rearrange the lattice upon hopping, $I_{if}$ is the hopping integral pre-
factor, describing the electronic wave-function overlap between neighboring
sites, $a_{if}$ is the orbital parameter describing the localization strength
of the electronic wave-function, $U_{if}$ is the energy barrier for the
hopping process, $\varepsilon_{i}$ is the pre-localization energy of the
electron at zero deformation, $E_{i}$ is the polaron stabilization energy,
$\mathbf{F}$ is an applied electrical field and $\textbf{r}_{i}$,
$\textbf{r}_{f}$ the position vector of the initial or final site
respectively.

###  2.1 Free polaron hopping

As a first case, we consider undoped, defect free stoichiometric LN. In this
case excess charge carriers introduced in the lattice e.g. by thermal
reduction treatments, form free polarons hopping randomly within the ordered
structure of the $\mathrm{Nb_{Nb}}$ sublattice. In these conditions, all the
sites are equivalent and the distribution of the waiting times between the
jumps satisfies the hypothesis of the Central Limit Theorem, which lead to the
normal diffusion behavior [36]. Macroscopically, this results in the standard
Fick’s law. In a 1-dimensional case:

$\frac{\partial P(z,t)}{\partial
t}=D_{\mathrm{FF}}\frac{\partial^{2}P(z,t)}{\partial z^{2}}$ (2)

where $D_{\mathrm{FF}}$ is the diffusion coefficient for free polarons and
$P(z,t)$ is the Particle Distribution Function (PDF). For free polarons, the
diffusion coefficient can be calculated in the following form [37]:

$D_{\mathrm{FF}}=\frac{1}{6}\sum_{j=1}^{N}r_{j}^{2}w_{\mathrm{FF}}(r_{j},T)=\frac{w_{\mathrm{FF}}^{0}(T)}{6}\sum_{j=1}^{N}r_{j}^{2}\exp{\left(-\frac{r_{j}}{a_{\mathrm{FF}}}\right)}$
(3)

where the summation runs over all the sites of the lattice, $r_{j}$ is the
distance from the origin to the $j$-th neighbour, and $w_{\mathrm{FF}}$ is the
hopping frequency of equation 1a. Even though the exponential decay of the
Holstein’s frequency guarantees that the summation converges quite fast, in
order to get accurate results it is necessary to consider beyond-nearest
neghbouts contributions. For practical purposes, the sum 3 it can be truncated
up to a certain coordination sphere around the starting site [4]. The solution
to eq. 3 at time $t$ for an initial delta-like distribution placed at the
origin $P(z,0)=\delta(z)$ is the Gaussian distribution:

$P(z,t)=\frac{1}{\sqrt{4\pi
D_{\mathrm{FF}}t}}\exp{\left(-\frac{z^{2}}{4D_{\mathrm{FF}}t}\right)}$ (4)

As it is well known, this distribution has zero mean and variance increasing
linearly with time of the form

$\sigma^{2}_{\rm FF}=2D_{\mathrm{FF}}t$ (5)

If an electric field is present, the average of the PDF remains no longer at
zero, but instead increases linearly with time with a given mobility $\mu$:

$<z>=\mu_{\mathrm{FF}}\mathrm{F}t.$ (6)

with $\mathrm{F}=|\bf{F}|$ the electric field strength, which here is assumed
oriented along $z$.

The knowledge of the diffusion coefficient is sufficient to describe
completely the mobility via the Einstein’s relationship
$\mu_{\mathrm{FF}}=eD_{\mathrm{FF}}/(k_{\mathrm{B}}T)$ as long as the
stochastic behaviour dominates the electric field one [38]. This condition is
satisfied when $e\bf{r\cdot F}$ $\ll k_{\mathrm{B}}T$.

###  2.2 Bound polaron hopping

Congruent lithium niobate is the most commonly encountered composition of this
material, since it is the energetically most stable and the one that is most
easily grown. It contains about $[\mathrm{Nb_{Li}}]=19.09\cdot
10^{25}\,\mathrm{m^{-3}}$ antisite defects [35]. Those antisites are
localization centers for B polaron formation. As it was reported in several
recent papers [4, 25, 26], for sufficiently low temperatures and congruent
composition, the main channel for polaron transport is direct hopping among
the antisites which can be viewed collectively as the sites of a disordered
lattice. Due to this increased disorder level and to the exponential
dependence of the hopping frequency, the Central Limit Theorem is no longer
applicable to the distribution of the waiting times and the normal diffusion
theory breaks down. The theoretical framework to deal with this problem is the
Continuous Time Random Walk (CTRW) model [36] in which a particle, after
hopping to a given site, waits a certain time before moving again. For this
situation, the Fick’s second law (Eq. 2) is replaced by a more general
equation governing the diffusion. In 1-D, it can be expressed in the following
form [36]:

$\frac{\partial P(z,t)}{\partial
t}={{}_{0}}D_{t}^{1-\alpha}D_{\mathrm{BB}}\frac{\partial^{2}P(z,t)}{\partial
z^{2}}$ (7)

where P(z,t) is the PDF, ${{}_{0}}D_{t}^{1-\alpha}$ is the fractional
derivative operator defined in Ref. [36], $D_{\rm BB}$ is the generalized
diffusion coefficient for bound polarons, having the dimension $[D_{\rm
BB}]=\frac{m^{2}}{s^{\alpha}}$ and $\alpha$ is the anomalous parameter
classifying the type of diffusion. For $\alpha=$1, Eq. (7) becomes the
standard Fick’s equation so that this case describes _normal_ diffusion while
all other cases are termed _anomalous_. Thus we may expect that, in contrast
to the free polaron case occurring in defect-free materials, bound polarons
obey eq. 7 if the proper conditions are met. In the case of bound polarons
diffusing in LN, as it will be shown below, it is found $\alpha<1$, which
corresponds to the so-called sub-diffusive case. The most general fundamental
solution of equation (7) for the sub-diffusive case is given in terms of the
Fox H-function [39] which can be expressed via the following series expansion:

$P(z,t)=\frac{1}{\sqrt{4D_{\rm
BB}t^{\alpha}}}\sum_{n=0}^{\infty}\frac{(-1)^{n}}{n!\Gamma(1-\alpha[n+1]/2)}\left(\frac{z^{2}}{D_{\rm
BB}t^{\alpha}}\right)^{n/2}$ (8)

The 1D variance of the distribution of equation (8) has the form of:

$\sigma^{2}_{\rm BB}(t)=\frac{2D_{\rm BB}t^{\alpha}}{\Gamma(1+\alpha)}$ (9)

and it is no longer proportional to the first power of time.

If the stochastic behaviour dominates the electric field one, it can be shown
that the _anomalous_ mobility can still be computed from the generalized
diffusion coefficient via the Einstein’s equation:

$\mu_{\alpha}=\frac{eD_{\rm BB}}{k_{B}T}$ (10)

The first moment of the distribution is now described by the relationship:

$\left\langle z(t)\right\rangle_{F}=\frac{F\mu_{\rm
BB}t^{\alpha}}{\Gamma(1+\alpha)}$ (11)

###  2.3 Free and bound mixed hopping

The most general situation, which is also the one most often encountered in
experiments at room temperature with congruent LN is the one where both bound
and free polarons contribute to the transport. In this case two more
elementary hopping processes besides FF and BB hopping are involved in the
transport, i.e. BF and FB. The idea here is to relate the two transport modes
detailed in sections 2.1 and 2.2 via some coupled transport models.

####  2.3.1 Switching Diffusion model.

Figure 1: Schematic representation of switching time between F (green dots)
and B (blue dots) polarons used on the MIM and SW model. In the scheme
depicted here $\tau_{\mathrm{FB}}=\tau_{1}+\tau_{2}+\tau_{3}$ and
$\tau_{\mathrm{BF}}=\tau_{4}+\tau_{5}$.

In the switching diffusion model, the population of F and B polarons is
coupled through some switching rates $k_{ij}$, whose physical meaning is
sketched in Fig. 1. They represent the inverse of the average time between
sequential changes of species, irrespective of the path that the polaron takes
between the species change: $\tau_{ij}=k_{ij}^{-1}$. In the case studied here,
the model described in Ref. [40] is used as reference and extended to our case
to account for the fact that F polarons obey the normal diffusion case, while
B follow an anomalous behavior. The coupled diffusion equations then take the
form:

$\frac{\partial P_{\rm F}(x,t)}{\partial t}=D_{\rm FF}\frac{\partial^{2}P_{\rm
F}(x,t)}{\partial x^{2}}-k_{\rm FB}P_{\rm F}(x,t)+k_{\rm BF}P_{\rm B}(x,t)$
(12a) $\frac{\partial P_{\rm B}(x,t)}{\partial t}={}_{0}D_{t}^{1-\alpha}D_{\rm
BB}\frac{\partial^{2}P_{\rm B}(x,t)}{\partial x^{2}}+k_{\rm FB}P_{\rm
F}(x,t)-k_{\rm BF}P_{\rm B}(x,t)$ (12b)

with $P_{\rm F}(x,t)$ and $P_{\rm B}(x,t)$ the PDFs of the free and bound
polarons and $D_{\rm FF}$ and $D_{\rm BB}$ their respective diffusion
coefficient. The Riemann-Liouville fractional derivative operator
${}_{0}D_{t}^{1-\alpha}$ in equation 12b takes into account the anomaly of the
bound polaron transport [36, 41].

Application of the Laplace transform
$\mathscr{L}\\{P_{i}(x,t)\\}=\hat{P_{i}}(x,s)$ simplifies the expressions
according to:

$s\hat{P_{\rm F}}(x,s)-P_{\rm F}(x,t=0)=D_{\rm
FF}\frac{\partial^{2}\hat{P_{\rm F}}(x,s)}{\partial x^{2}}-k_{\rm
FB}\hat{P_{\rm F}}(x,s)+k_{\rm BF}\hat{P_{\rm B}}(x,s)$ (13a) $s\hat{P_{\rm
B}}(x,s)-P_{\rm B}(x,t=0)=s^{1-\alpha}D_{\rm BB}\frac{\partial^{2}\hat{P_{\rm
B}}(x,s)}{\partial x^{2}}+k_{\rm FB}\hat{P_{\rm F}}(x,s)-k_{\rm BF}\hat{P_{\rm
B}}(x,s)$ (13b)

where we applied
$\mathscr{L}\\{{}_{0}D_{t}^{1-\alpha}f(x,t)\\}=s^{1-\alpha}\hat{f}(x,s)$ [36].
This expression is analytically solvable in the Laplace space and a real space
solution can be derived by numerical inverse Laplace transform. The process is
performed separately for the two cases of polarons starting as F or as B
polarons respectively. In the former case the boundary conditions are chosen
as $P_{\rm F}(x,t=0)=\delta(x)$ and $P_{\rm B}(x,t=0)=0$ and reciprocally
$P_{\rm F}(x,t=0)=0$ and $P_{\rm B}(x,t=0)=\delta(x)$ for the latter, with
$\delta(x)$ the Dirac delta function. Additionally, for all cases
$\lim_{b\to\infty}P_{\rm F}(\pm b,t)=\lim_{b\to\infty}P_{\rm B}(\pm b,t)=0$ is
chosen as an additional boundary condition. It is then possible to have access
to both PDFs from which the variance of the distribution can be computed as a
function of time as well as the mobility thanks to the Einstein’s equation.

####  2.3.2 Mobile - Immobile model.

The previous treatment allow for an exact numerical solution of the problem,
but does not have an analytical solution. Since the hopping frequency of bound
polarons for all the processes considered here is several orders of magnitude
smaller than the one of free polarons, we expect that a good approximation of
the previous situation may be provided by a Mobile-Immobile diffusion model
(MIM)[42] which on the contrary is analytically solvable. In the MIM model,
free polarons are considered the mobile specie, normally diffusing with
$D_{\mathrm{FF}}>0$ and bound polaron are considered immobile with
$D_{\mathrm{BB}}\approx 0$. The mobile species immobilize with a rate $k_{\rm
FB}$ and become mobile again with a rate $k_{\rm BF}$. In this case the
diffusion equations take the form:

$\frac{\partial P_{\rm F}(x,t)}{\partial t}=D_{\rm FF}\frac{\partial^{2}P_{\rm
F}(x,t)}{\partial x^{2}}-k_{\mathrm{FB}}P_{\rm F}(x,t)+k_{\mathrm{BF}}P_{\rm
B}(x,t)$ (14a)

$\frac{\partial P_{\rm B}(x,t)}{\partial t}=k_{\mathrm{FB}}P_{\rm
F}(x,t)-k_{\mathrm{BF}}P_{\rm B}(x,t)$ (14b)

In the MIM model, the variances of the distribution can be computed
analytically as [42]:

$\sigma^{2}=\frac{2D_{\rm
FF}}{1+k_{\mathrm{FB}}/k_{\mathrm{BF}}}\left[t+\frac{k_{\mathrm{FB}}/k_{\mathrm{BF}}^{2}}{1+k_{\mathrm{FB}}/k_{\mathrm{BF}}}\left(1-e^{-(k_{\mathrm{BF}}+k_{\mathrm{FB}})t}\right)\right]$
(15a) $\sigma^{2}=\frac{2D_{\rm
FF}}{1+k_{\mathrm{FB}}/k_{\mathrm{BF}}}\left[t-\frac{k_{\mathrm{BF}}^{-1}}{1+k_{\mathrm{FB}}/k_{\mathrm{BF}}}\left(1-e^{-(k_{\mathrm{BF}}+k_{\mathrm{FB}})t}\right)\right]$
(15b)

for free and bound starting site, respectively.

The displacement of the polaron distribution under the effect of a bias can be
readily computed from the Einstein’s relation providing:

$<z>=\frac{eF\sigma^{2}}{2k_{B}T}$ (16)

where $\sigma$ stems for the variance in one of the two cases of initial F- or
B- polaron population. Clearly, in this case we shall define a time-dependent
_instantaneous_ mobility $\mu(t)=\frac{d\langle z\rangle}{dt}$ as the average
displacement is no longer proportional to time.

##  3 Monte Carlo simulations of polaron hopping transport in LN

Monte Carlo (MC) simulations have proven themselves a valuable tool in
modelling polaron transport in LN [25, 4]. A dedicated Monte Carlo algorithm
is developed to study the material transport properties using the same
methodology detailed in the cited papers. The code generates a stoichiometric
$\mathrm{LiNbO_{3}}$ structure with periodic boundary conditions and
subsequently a certain number of $\mathrm{Nb_{Li}}$ antisites is randomly
placed in the structure in the appropriate lattice sites, according to the
crystal composition.

The hopping frequencies for a given initial site $i$ towards a destination $f$
included in a suitable volume around the starting site are computed by eq. 1a
using the values reported in Tab. 1 [25, 4, 3]. The final destination is
chosen via the Gillespie algorithm. The time needed to perform the hop is then
recorded and the process is repeated until a given time has elapsed. A new
polaron is then created and the whole simulation iterated until a satisfactory
statistics has been achieved. For each run, the program records the final
polaron positions with respect to its original site, the number of different
sites encountered during its walk and the time needed to reach the final site.
All the information is then analyzed to compute the variance and the mean
displacement of the final polaron distribution as a function of time. In the
simulation we consider the possible presence of a static electric field F
directed along the +z direction to study polaron mobility, while for the study
of diffusion processes the field is set at zero.

Parameters | Value | Reference
---|---|---
$E_{\mathrm{F}}$ | 0.545 eV | [3]
$E_{\mathrm{B}}$ | 0.750 eV | [25]
$\varepsilon_{\mathrm{F}}$ | 0.000 eV | [3]
$\varepsilon_{\mathrm{B}}$ | 1.690 eV | [25]
$a_{\mathrm{FF}}=a_{\mathrm{FB}}=a_{\mathrm{BB}}=a_{\mathrm{BF}}$ | 1.60 Å | [25]
$I_{\mathrm{FF}}=I_{\mathrm{FB}}=I_{\mathrm{BB}}=I_{\mathrm{BF}}$ | 0.30 Å | [25]
Table 1: Parameters used to compute the hopping frequency in Eq. 1a

##  4 Results

###  4.1 Free polaron hopping

Figure 2: (a) Variance of the final electron position as a function of the
time on an ordered $\mathrm{Nb_{Nb}}$ structure at $\mathrm{T=300\,K}$ (b)
Normalized probability distribution at two different times of the population
of F polarons diffusing on the same lattice as in (a) at $\mathrm{T=300\,K}$.
Both histograms are fitted with a Gaussian distribution (black curves).

As a first case, we consider a system of free polarons hopping in the
$\mathrm{Nb_{Nb}}$ sublattice. Fig. 2a) show the variance of the polaron
distribution as a function of the time along the $x$, $y$, and $z$
crystallographic directions in a stoichiometric LN crystal at room
temperature. As expected, the variance increases linearly with time. This is
the fingerprint that the polaron is normally diffusing, as described in
section 2.1. The linear fit of the variance gives a diffusion coefficient of
$D_{\rm FF,x}$= (1.89 $\pm$ 0.04)$\cdot 10^{9}$ Å$\mathrm{{}^{2}/s}$, $D_{\rm
FF,y}$= (1.87 $\pm$ 0.04)$\cdot 10^{9}$ Å$\mathrm{{}^{2}/s}$, $D_{\rm FF,z}$=
(1.89 $\pm$ 0.03)$\cdot 10^{9}$ Å$\mathrm{{}^{2}/s}$. This values are in
agreement with the one calculated via equation 3, being equal to $D_{\rm
FF}=1.86\cdot 10^{9}$ Å$\mathrm{{}^{2}/s}$. This is obtained by summing the
contribution of the neighbours till $r_{j}\approx 80$ Å corresponding to the
8th coordination shell, as already found in [26]. Note that limiting the
summation in eq. 3 to the first neighbors would have seriously underestimated
this result.

The temperature dependence of the diffusion coefficient is ruled by the term
1b which can be factored out from the summation in eq. 3. Thus for free
polarons the diffusion is thermally activated with a characteristic energy
$U_{FF}=0.27$ eV given by eq. 1c.

Figure 3: (a) Diffusion coefficient in function of the electric field along x,
y and z. Horizontal black line corresponds to the value computed with F=0. (b)
Mean value of the distribution of the final position in function of the time
for the case $\mathrm{F=1\cdot 10^{7}}$ V/m.

Let us now consider the effect of an electric field on the free polaron PDF.
In fig. 3 (b) it is shown the motion of the polaron distribution as a function
of time for a given applied field. As expected the first moment of the polaron
distribution is linearly increasing with time, which once again confirms that
free polarons obey normal diffusion laws. The slope corresponds to a free
polaron mobility at room temperature in stoichiometric LN equal to
$\mu=(7.30\pm 0.3)\cdot 10^{10}$ Å2/(Vs). This result is in agreement with the
one expected from the Einstein’s relation as long as the applied field is
below a value of $\mathrm{F=1\cdot 10^{7}}$ V/m, see fig. 3 (a). For higher
field values the variance of the distribution starts to increase because the
condition $e\bf{r\cdot F}$ $\ll k_{\mathrm{B}}T$ is no longer satisfied [38].
It is interesting to note that experimentally measured internal fields in cLN
at room temperature are generally below this value [43, 4], so that Einstein’s
relation can indeed be assumed valid for free polaron transport in realistic
conditions.

###  4.2 Bound polarons hopping

Figure 4: (a) Variance of the final electron position as a function of the
time for a case of a pure BB jump mechanism at $\mathrm{T=300\,K}$ (b)
Probability distribution function of the polaron position at a fixed time at
$\mathrm{T=300\,K}$. The histogram is described analytically by the Fox
function.

Fig. 4a) shows the variance of the distribution of the final polaron position
along $x$, $y$, and $z$ as a function of time in a case of a pure BB jump
mechanism for which free polarons are forbidden in the simulation. The
relationship between these two quantities is not linear anymore, confirming
that bound polarons are anomalously diffusing on the lattice [36]. The
simulation results are in agreement with Eq. 9 with $\alpha=0.72\pm 0.01$ and
$D_{\rm BB}=(7.1\pm 0.4)\cdot 10^{3}\,\AA^{2}/\mathrm{s}^{\alpha}$. The value
of $\alpha$ reveals that the bound polarons are sub-diffusing. The actual
value of this parameter is linked to the antisite concentration in the crystal
and on the $a_{BB}$ value but not on temperature. Also in this case, the
thermal dependence of the diffusion coefficient can be factored out in the
summation of eq. 1b. Thus, for pure bound polaron transport, the diffusion
coefficient is thermally activated with energy $U_{\rm BB}$ (see eq. 1c). Fig.
4b) shows the simulated PDF for the bound polarons population at a fixed time.
The black line is the the Fox’ function reported in equation 8 with the same
$\alpha$ and $D_{\rm BB}$ as reported above, which is in good agreement with
simulation results, as expected [36].

Figure 5: (a) Diffusion coefficient in function of the electric field along x,
y and z. Horizontal line corresponds to the value computed with F=0. (b) Mean
value of the distribution of the final position in function of the time for
the case $\mathrm{F=5\cdot 10^{6}}$ V/m.

The effect of an electric field on a polaronic system composed of only bound
polarons is reported in Fig. 5 (b). The PDF average is again changing with
time under the effect of the field, but this time with a sub-linear behaviour
as expected for anomalous diffusion. From the fit with equation 11, the
anomalous mobility and the anomalous coefficient are equal to
$\mu_{\alpha}=(2.6\pm 0.4)\cdot 10^{5}\AA^{2}/\mathrm{V}\mathrm{s}^{\alpha}$
and $\alpha=0.71\pm 0.04$, in agreement with the Einstein’s relation. Also in
this case, when the electric field is too strong, the above-mentioned
behaviors are no longer observed. Fig. 5a) shows the diffusion coefficient
computed in function of the electric field at room temperature. The horizontal
line shows the values of the diffusion coefficient computed in the situation
$\mathrm{F=0}$ for comparison. Again, above $\mathrm{F\approx 1\cdot
10^{7}\,V/m}$ the variance of the PDF begins to be field-dependent, with a
steep increase.

###  4.3 Congruent LN

In this section, the real-life case of a cLN crystal is considered. As
detailed in the Introduction, this material is characterized by a high
concentration of intrinsic defects, namely Li vacancies ($\mathrm{V_{Li}}$)
and Nb antisites ($\mathrm{Nb_{Li}}$). In this case polarons do not only hop
among alike sites, but conversion processes where a free polaron is captured
by an antisite defect becoming a bound one and vice-versa are present, mixing
up the transport processes described in previous sections. The relative weight
of these processes depends not only on the sample composition, but also on
temperature [4] as the activation energies for the various processes are
different, as per eq. 1c. Two situations may be considered depending whether
the initial delta-like polaron distribution is assumed to be in the free
polaron state or in the bound one. This is very much dependent on the
experimental conditions and on the time - scale of the experiment one is
looking at. In the following both situations are simulated and discussed
focusing the attention to the case of a congruent LN crystal with an antisite
concentration $[\mathrm{Nb_{Li}}]=1.9\cdot 10^{20}\,\mathrm{cm}^{-3}$ at room
temperature.

Figure 6: Variance of the polaron PDF as a function of the time in a congruent
LN crystal at $\mathrm{T=300\,K}$ (blue dots). Theoretical variance according
to the switching diffusion model (SW, green dots). Theoretical variance using
the MIM model, eq. 15a (yellow lines).

Fig. 6 shows the variance of the polaron PDF along the z direction in function
of the time for the case that either F or B is the starting polaron of the
hopping. The small blue dots are the results obtained from the MC simulations.
If polarons begin their journey as F-polarons (upper curve), three distinct
transport phases are visible. Initially, the variance increases almost
linearly in the double-logarithmic representation, indicating a power law
behaviour. This initial transport phase continues until around 10-7 s when a
static phase follows. The transport is then frozen for almost two decades of
time, from around $10^{-7}$ s until $10^{-5}$ s when the diffusion starts
again. In this third phase, the variance evidences again a power-law
behaviour. If polarons start as B-polaron (lower blue dots), the PDF basically
does not change until $10^{-6}$ s because in such a short time frame the
possibility that polarons may escape from the antisites is negligible. From
here on, the variance curve shows a linear dependence in the double log plot.
This linearity continues until around $10^{-3}$ s where a slight slope change
is visible. From here, the evolution of the variance reaches asymptotically
the same one observed for a starting F-polaron distribution. The physical
interpretation of these results is quite straightforward. In the case of an
initial F-polaron distribution, the particles start to quickly move on the
$\mathrm{Nb_{Nb}}$ sublattice until they are gradually trapped by
$\mathrm{Nb_{Li}}$ antisites that slow down the diffusion. After this initial
stage, polaron motion may occur either by rare $B\rightarrow B$ hopping events
or by conversion to free polarons followed by a sequence of hops:
$B\rightarrow F\rightarrow F\cdots\rightarrow B$ .

The green large dots in fig. 6 are the numerical solutions for the variance
obtained from eq. 12b assuming the parameters discussed in the previous
sections ($D_{\rm FF}=1.89\cdot 10^{9}$ Å2/s, $D_{\rm BB}=7.3\cdot 10^{3}$
Å${}^{2}/s^{\alpha}$, $\alpha=0.71$) for pure F and B transport. The switching
rates ($k_{\rm FB}=3.76\cdot 10^{9}\,\mathrm{s}^{-1}$ , $k_{\rm BF}=2.14\cdot
10^{3}\,\mathrm{s}^{-1}$) are computed by comparing the MC simulation results.
The characteristic times $\tau_{FB(BF)}=k_{FB(BF)}^{-1}$ mark the positions
along the time axis of the knees of the variance curve. Finally, the yellow
lines in fig. 6 are calculated using equations 15a, 15b using the same
parameters used for the switching diffusion model. As it can be seen, the MIM-
model is a good approximation of the observed trends.

Figure 7: a) Average value of polaron position in function of the time when
the starting site is a free polaron (upper blue dots) or bound polaron (lower
blue dots) computed for $F=5\cdot 10^{6}$V/m. Data are superimposed to
equation 16 with the parameters as described in the text.

Considering now the case of an applied electric field, Fig. 7 shows the
average displacement of the polaron PDF along the polar axis as a function of
the time under the effect of an applied bias, here assumed $\mathrm{F=5\cdot
10^{6}\,V/m}$. Again, two cases are considered for an initial F-polaron
population (upper curve) or for a B-polaron population (lower curve). After an
initial swift motion with high mobility equal to $\mu_{FF}$, the polarons
remain still under the effect of the field unless sufficient time is elapsed
so that polarons can de-trap from antisites and move as free polarons or
directly hop on another antisite. As observed in precedent paragraphs this
value is sufficiently low so that Einstein’s relationship should hold, which
is confirmed by the fact that the variance of the polaron PDF under bias is
unchanged with respect to the one considered in fig. 6. The simulation results
are therefore superposed to eq. 16 (yellow curve in Fig. 7). The parameters
are $\mathrm{\tau_{FB}=(2.1\pm 0.6)\cdot 10^{-8}\,s}$ and
$\mathrm{\tau_{BF}=(3\pm 0.7)\cdot 10^{-4}\,s}$, indicated by the vertical
black lines in the figure and in agreement with previous results.

##  5 Discussion

| $t\ll\tau_{FB}$ | $t\gg\tau_{BF}$
---|---|---
Stoichiometric LN | $\mu_{FF}$ | $\mu_{FF}$
Congruent LN, F starting site | $\mu_{FF}$ | $\mu_{eff}$
Congruent LN, B starting site | $\approx 0$ | $\sim\mu_{eff}$
Table 2: Limiting values of the polaron mobility for different experimental
cases at $T=300K$. The parameters are: $\mathrm{\tau_{FB}=(2.1\pm 0.6)\cdot
10^{-8}\,s}$; $\mathrm{\tau_{BF}=(3\pm 0.7)\cdot 10^{-4}\,s}$;
$\mu_{FF}=7.3\cdot 10^{4}\,\mathrm{cm}^{2}\mathrm{V}^{-1}\mathrm{s}^{-1}$;
$\mu_{eff}=5.1\,\mathrm{cm}^{2}\mathrm{V}^{-1}\mathrm{s}^{-1}$.

In the case of defect-free lithium niobate the polaron diffusion and mobility
are normal, as expected from theoretical models and verified by Monte-Carlo
simulations. The polaron mobility is thermally activated with an energy
$U_{FF}=0.27$ eV and has a value of $\mu_{\rm FF}=7.3\cdot
10^{4}\mathrm{cm}^{2}/Vs$ at room temperature, corresponding to a diffusion
coefficient $D_{\rm FF}=1.89\cdot 10^{9}\AA^{2}/s$. When the more realistic
situation of a non-stoichiometric lithium niobate sample is considered, the
phenomenology becomes significantly reacher. The system can be well described
by a switching diffusion model, embodied by eqs. 12b which can in principle be
solved numerically, but for which the rate constants and the anomalous
coefficients parameters are difficult to compute analytically. Our MC
simulations provide for these quantities $\alpha=0.72\pm 0.01$ and
$\mathrm{D_{BB}}=(7.1\pm 0.4)\cdot 10^{3}\,\AA^{2}/\mathrm{s}^{\alpha}$.
However we showed that for a LN sample in standard conditions, i.e at room
temperature and congruent composition (i.e. $\mathrm{Nb_{Li}}=1.9\cdot
10^{20}\mathrm{cm}^{-3}$) the influence of direct $B\rightarrow B$ hops is
generally small, therefore the system can be described with a good accuracy in
terms of a Mobile - Immobile model for which the polaron diffusion and
mobility are described by eqs. 15a, 15b and 16. These equations show that the
polarons behave differently depending on the experimental situation and on the
time scale considered, as detailed below. The main results obtained in the MIM
approximation are reported in table 2.

###  5.1 F starting site

In experiments where a polaron population is suddenly created e.g. by means of
a short laser pulse such as in fs- or TAS spectroscopy, it can be assumed that
just after the pulse the majority of the polarons are free. For short times
$t\ll\tau_{FB}=\mathrm{k_{FB}}^{-1}$ it is easy to check that eq. 15a can be
approximated as $\sigma^{2}\sim 2D_{FF}t$, so that the system behaves as a
defect-free material, with the same mobility.

Figure 8: Average number of jumps per unit time $\overline{N}_{ij}$, resolved
for the different transition types in congruent LN as at T=300K as function of
simulation time. a) For the case of polarons starting as F and b) for polarons
starting as B. Vertical dashed lines represent $\tau_{\rm FB}$ and $\tau_{\rm
BF}$.

For $\tau_{\mathrm{FB}}<t<\tau_{\mathrm{BF}}$ F- polarons are almost
completely transformed in B- polarons which diffuse much slower than F
polarons, so that in this stage the charge diffusion is abruptly decelerated.
In order to investigate in greater detail the type of conduction in this time
frame, we can analyze the relative abundance of the different hop types as
retrieved from the MC simulations. In fig. 8 (a) the amount of hops per unit
time for the four elementary hopping processes is shown as a function of time.
As it can be seen after the initial stage, the number of FF hops drops down to
about $10^{4}\,\mathrm{hops}/s$ in the intermediate regime, while
$\overline{N}_{FB}=\overline{N}_{BF}$, indicating that the F and B polarons
populations are in equilibrium. F polarons do not contribute anymore to the
transport, polarons are simply shifting quickly between antisites and regular
sites without transport. In addition to those processes, we observe a constant
amount of BB hops which is the responsible for the diffusion. Thus,
experiments working in this time frame would see the charge carriers diffuse
very slowly and anomalously as B polarons with the same parameters described
in sec. 4.2. If this contribution is neglected, we end up with the MIM model
description. In this time range eq. 15a gives $\sigma^{2}\sim 2D_{\rm
FF}\tau_{FB}=2L_{\rm D}^{2}$ where $L_{\rm D}$ is the diffusion length that
can be measured in holographic experiments (see e.g. [44]) and this case equal
to $L_{\rm D}\sim 7\AA$. The effect of an applied bias in this time window
results in a change of the average polaron position reached when the PDF
”freezes” because of the antisites. The plateau value visible in Fig. 7 can be
computed within the MIM model as:

$\langle z\rangle_{P}=\frac{eL_{\rm D}^{2}}{k_{B}T}F$ (17)

Finally, for long times $t\gg\tau_{BF}$, detrapping processes are sufficiently
probable so that the chance to have a few polarons contributing to the
transport with FF hops is non negligible. Although their number is very
limited, the difference in mobility is so large that the FF contribution
becomes visible and eventually the dominating one, so that at long times
diffusion appears normal again, with an offset deriving from the previous
stages. In the MIM approximation, we obtain from eq. 15a:

$\sigma^{2}\sim 2D_{\rm eff}t+2L_{\rm D}^{2}$ (18)

with $D_{\mathrm{eff}}\approx D_{\rm FF}\frac{\tau_{FB}}{\tau_{BF}}=(1.25\pm
0.1)\cdot 10^{5}$ Å2/s i.e. about four orders of magnitude smaller than
$D_{\rm FF}$. In the regime of validity of the Einstein’s relation, this
corresponds also to the difference in the long-time limit of the polaron
mobility between a stoichiometric and a congruent LN crystal at room
temperature. Indeed, for sufficiently long times, the average displacement
under the effect of a field is:

$\langle z\rangle\approx\frac{eD_{\rm eff}}{kT}Ft+\langle
z\rangle_{P}=\mu_{\rm eff}Ft+\langle z\rangle_{P}\sim\mu_{\rm eff}Ft$ (19)

with an effective polaron mobility $\mu_{\rm
eff}=5.1\,\mathrm{cm}^{2}\mathrm{V}^{-1}\mathrm{s}^{-1}$.

###  5.2 B starting site

Another commonly encountered experimental case is the one of reduced LN
samples of congruent composition. If the reduction degree is not too high,
charge transport can be attributed to small bound polarons [45] which are much
more abundant in the lattice than free polarons. Bipolarons are present as
well [3], but they do not appear to contribute significantly to charge
transport and in the following will be disregarded.

At short times ($t\ll\tau_{BF}$), the processes that happen more frequently
are BB hops, typically those involving $\mathrm{Nb_{Li}}$ antisites at close
distance one to the other. In this regime $D\propto t^{\alpha}$, i.e. the
diffusion is anomalous with $\alpha=0.71$, as visible in Fig. 4. While eqs.
12b may be used for an accurate calculation, the MIM model, which neglects the
bound polaron contribution, underestimates the diffusion coefficient in this
time frame. For $t\gg\tau_{BF}$ the role of free polarons is more and more
impactful and the situation ends up to the same one considered the previous
section, as described by eq. 18, with an effective diffusion coefficient
$D_{\rm eff}$. As before, the polaron mobility (as shown in Fig. 7) follows
the behaviour expected from the diffusion coefficient as long as the
Einstein’s relation holds true.

##  6 Conclusions

The motion of polarons, described by a random hopping process among the sites
of a polar crystals is the underlying physical process establishing the
electrical properties of the material, which are of paramount importance for
several applications. In this work the role of defects on polaron diffusion
and mobility has been explored theoretically for the specific case of lithium
niobate. We have shown that the presence of defects requires a significant
revision of the theoretical framework needed to describe the polaron motion
and obtained realistic estimates of the physical quantities necessary to
describe quantitatively the process. In particular, the presence of intrinsic
defects, as those encountered in the commonly used LN with congruent
composition, introduces a time-dependent diffusion and mobility. Depending on
the initial sample condition and on the time scale of the experiment one is
interested in, several regimes have been evidenced, with the defects affecting
the transport behavior to increasing extents. Our theoretical models,
supported by MC simulations based on the Marcus-Holstein hopping frequency,
are able to describe quantitatively the situation under consideration,
providing both formal expressions and numerical parameters for the case of
congruent LN at room temperature. In particular, the Mobile - Immobile
description provides reliable and simple analytical expressions for diffusion
and mobility in a time range from $10^{-8}$ s up to the seconds and beyond.

These expressions should replace the ones used in general formulations for
diffusion and mobility whenever non-stoichiometric LN is considered. Our
formalism can be readily extended to other temperatures, provided that the
rate constants entering eqs. 12b are known. Furthermore, our approach may be
extended to include additional defects and or dopants with the capability to
affect the polaron hopping process, such as Fe traps in the case of Fe-doped
lithium niobate or $\mathrm{Ta_{Nb}}$ substitutional defects as in the mixed
lithium niobate - tantalate system.

## References

## References

  * [1] Cesare Franchini, Michele Reticcioli, Martin Setvin, and Ulrike Diebold.  Polarons in materials.  Nature Reviews Materials, 6(7):560–586, 2021.
  * [2] Mario Maglione.  Free charge localization and effective dielectric permittivity in oxides.  Journal of Advanced Dielectrics, 06(02):1630006, 2016.
  * [3] O. F. Schirmer, M. Imlau, C. Merschjann, and B. Schoke.  Electron small polarons and bipolarons in LiNbO3.  Journal of Physics: Condensed Matter, 21(12):123201, 2009.
  * [4] Laura Vittadello, Marco Bazzan, Anush Danielyan, Edvard Kokanyan, Laurent Guilbert, and Michel Aillerie.  A polaron approach to photorefractivity in Fe:LiNbO3.  Journal of Physics Communications, 2(12):125003, 2018.
  * [5] L. Grigorjeva, D.K. Millers, V. Pankratov, R.T. Williams, R.I. Eglitis, E.A. Kotomin, and G. Borstel.  Experimental and theoretical studies of polaron optical properties in knbo3 perovskite.  Solid State Communications, 129(11):691–696, 2004.
  * [6] S. Torbrügge, M. Imlau, B. Schoke, C. Merschjann, O. F. Schirmer, S. Vernay, A. Gross, V. Wesemann, and D. Rytz.  Optically generated small electron and hole polarons in nominally undoped and Fe-doped KNbO3 investigated by transient absorption spectroscopy.  Physical Review B, 78(12), sep 2008.
  * [7] L. A. Kappers, K. L. Sweeney, L. E. Halliburton, and J. H. W. Liaw.  Oxygen vacancies in lithium tantalate.  Physical Review B, 31(10):6792–6794, 1985.
  * [8] Youwen Liu, Kenji Kitamura, Shunji Takekawa, Masaru Nakamura, Yasunori Furukawa, and Hideki Hatano.  Two-color photorefractive properties in near-stoichiometric lithium tantalate crystals.  Journal of Applied Physics, 95(12):7637–7644, 2004.
  * [9] Naoki Tsunoda, Yu Kumagai, and Fumiyasu Oba.  Stabilization of small polarons in BaTiO3 by local distortions.  Physical Review Materials, 3(11):114602, 2019.
  * [10] E. Possenriede, H. Kröse, T. Varnhorst, R. Scharfschwerdt, and O. F. Schirmer.  Shallow acceptor and electron conduction states in batio3.  Ferroelectrics, 151(1):199–204, 1994.
  * [11] Elaheh Ghorbani, Lorenzo Villa, Paul Erhart, Andreas Klein, and Karsten Albe.  Self-consistent calculations of charge self-trapping energies: A comparative study of polaron formation and migration inPbTiO3.  Physical Review Materials, 6(7):074410, 2022.
  * [12] V. N. Bogomolov and N. Mirlin.  Optical absorption by polarons in rutile (TiO2) single crystals.  In Volume 27, Number 1 May 1, pages 443–454. De Gruyter, January 2022.
  * [13] Alex J. Tanner, Bo Wen, Yu Zhang, Li-Min Liu, Helen H. Fielding, Annabella Selloni, and Geoff Thornton.  Photoexcitation of bulk polarons in rutile ${Ti}{O}_{2}$.  Physical Review B, 103(12):L121402, March 2021.
  * [14] H.L. Tuller and A.S. Nowick.  Small polaron electron transport in reduced ceo2 single crystals.  Journal of Physics and Chemistry of Solids, 38(8):859–867, 1977\.
  * [15] José J. Plata, Antonio M. Márquez, and Javier Fdez. Sanz.  Electron Mobility via Polaron Hopping in Bulk Ceria: A First-Principles Study.  The Journal of Physical Chemistry C, 117(28):14502–14509, 2013\.
  * [16] Niklas Dömer, Julian Koelmann, Mira Hesselink, Tobias Hehemann, Anton Pfannstiel, Felix Sauerwein, Laura Vittadello, Steffen Ganschow, and Mirco Imlau.  Long-lived, pulse-induced absorption in LiNb1-xTaxO3 solid solutions: the case of three intrinsic defect sites for electron localization with strong coupling, March 2024.  type: article.
  * [17] Huimin Liu, Richard C. Powell, and L. A. Boatner.  Effect of niobium doping on the properties of picosecond laser-induced transient gratings in KT1-xNbxO3.  Physical Review B, 44(6):2461–2469, 1991.
  * [18] T. Holstein.  Studies of polaron motion.  Annals of Physics, 8(3):343–389, 1959.
  * [19] I. G. Austin and N. F. Mott.  Polarons in crystalline and non-crystalline materials.  Advances in Physics, 50(7):757–812, 2001.
  * [20] Polarons in Advanced Materials.  Springer Netherlands, 2007.
  * [21] J Koppitz, O. F Schirmer, and A. I Kuznetsov.  Thermal dissociation of bipolarons in reduced undoped LiNbO3.  Europhysics Letters (EPL), 4(9):1055–1059, 1987.
  * [22] Falko Schmidt, Agnieszka L. Kozub, Timur Biktagirov, Christof Eigner, Christine Silberhorn, Arno Schindlmayr, Wolf Gero Schmidt, and Uwe Gerstmann.  Free and defect-bound (bi)polarons in LiNbO3 : Atomic structure and spectroscopic signatures fromab initiocalculations.  Physical Review Research, 2(4), October 2020.
  * [23] S Messerschmidt, A Krampf, F Freytag, M Imlau, L Vittadello, M Bazzan, and G Corradi.  The role of self-trapped excitons in polaronic recombination processes in lithium niobate.  Journal of Physics: Condensed Matter, 31(6):065701, 2018.
  * [24] I. Mhaouech and L. Guilbert.  Temperature dependence of small polaron population decays in iron-doped lithium niobate by monte carlo simulations.  Solid State Sciences, 60:28–36, 2016.
  * [25] Laurent Guilbert, Laura Vittadello, Marco Bazzan, Imed Mhaouech, Simon Messerschmidt, and Mirco Imlau.  The elusive role of NbLi bound polaron energy in hopping charge transport in Fe:LiNbO3.  Journal of Physics: Condensed Matter, 30(12):125701, 2018.
  * [26] Laura Vittadello, Laurent Guilbert, Stanislav Fedorenko, and Marco Bazzan.  Polaron trapping and migration in iron-doped lithium niobate.  Crystals, 11(3):302, March 2021.
  * [27] Laura Vittadello, Annamaria Zaltron, Nicola Argiolas, Marco Bazzan, Nicola Rossetto, and Raffaella Signorini.  Photorefractive direct laser writing.  Journal of Physics D: Applied Physics, 49(12):125103, 2016-02.
  * [28] Alfonso Blázquez-Castro, Juan C. Stockert, Begoña López-Arias, Angeles Juarranz, Fernando Agulló-López, Angel García-Cabañes, and Mercedes Carrascosa.  Tumour cell death induced by the bulk photovoltaic effect of LiNbO3:Fe under visible light irradiation.  Photochemical & Photobiological Sciences, 10(6):956–963, 2011\.
  * [29] L. Lucchetti, K. Kushnir, A. Zaltron, and F. Simoni.  Liquid crystal cells based on photovoltaic substrates.  Journal of the European Optical Society: Rapid Publications, 11, 2016.
  * [30] Javier Villarroel, Héctor Burgos, Ángel García-Cabañes, Mercedes Carrascosa, Alfonso Blázquez-Castro, and Fernando Agulló-López.  Photovoltaic versus optical tweezers.  Optics Express, 19(24):24320, 2011.
  * [31] Zezheng Wu, Hongqiang Jiang, Lingling Zhang, Kezhen Yi, Heng Cui, Fubing Wang, Wei Liu, Xingzhong Zhao, Fuling Zhou, and Shishang Guo.  The acoustofluidic focusing and separation of rare tumor cells using transparent lithium niobate transducers.  Lab on a Chip, 19(23):3922–3930, 2019.
  * [32] Leonardo Zanini, Annamaria Zaltron, Enrico Turato, Riccardo Zamboni, and Cinzia Sada.  Opto-microfluidic integration of the bradford protein assay in lithium niobate lab-on-a-chip.  Sensors, 22(3):1144, 2022.
  * [33] Sebastian Cremaschini, Alberto Cattelan, Davide Ferraro, Daniele Filippi, Filippo Marinello, Alessio Meggiolaro, Matteo Pierno, Cinzia Sada, Annamaria Zaltron, Paolo Umari, et al.  Trifurcated splitting of water droplets on engineered lithium niobate surfaces.  ACS Applied Materials & Interfaces, 16(3):4271–4282, 2024.
  * [34] R. Zamboni, A. Zaltron, M. Chauvet, and C. Sada.  Real-time precise microfluidic droplets label-sequencing combined in a velocity detection sensor.  Scientific Reports, 11(1), 2021-09.
  * [35] Tatyana Volk.  Lithium Niobate.  Springer-Verlag Berlin Heidelberg, 2008.
  * [36] Ralf Metzler and Joseph Klafter.  The random walk's guide to anomalous diffusion: a fractional dynamics approach.  Physics Reports, 339(1):1–77, 2000.
  * [37] J Metzger.  Effective-medium theory for further-neighbour hopping transport.  Journal of Physics C: Solid State Physics, 17(29):5153–5161, 1984\.
  * [38] A. V. Nenashev, F. Jansson, S. D. Baranovskii, R. Österbacka, A. V. Dvurechenskii, and F. Gebhard.  Effect of electric field on diffusion in disordered materials. II. two- and three-dimensional hopping transport.  Physical Review B, 81(11):115204, March 2010.
  * [39] Elliott W. Montroll and George H. Weiss.  Random walks on lattices. $\mathrm{II}$.  Journal of Mathematical Physics, 6(2):167–181, 1965.
  * [40] Denis S. Grebenkov.  Time-averaged mean square displacement for switching diffusion.  Physical Review E, 99(3):032133, 2019.
  * [41] Ralf Metzler, Ashish Rajyaguru, and Brian Berkowitz.  Modelling anomalous diffusion in semi-infinite disordered systems and porous media.  New Journal of Physics, 24(12):123004, 2022-12.
  * [42] Timo J. Doerries, Aleksei V. Chechkin, and Ralf Metzler.  Apparent anomalous diffusion and non-gaussian distributions in a simple mobile–immobile transport model with poissonian switching.  Journal of The Royal Society Interface, 19(192), jul 2022.
  * [43] K. Peithmann, A. Wiebrock, and K. Buse.  Photorefractive properties of highly-doped lithium niobate crystals in the visible and near-infrared.  Applied Physics B: Lasers and Optics, 68(5):777–784, 1999.
  * [44] Jaime Frejlich.  Photorefractive Materials.  Wiley-Interscience, 2006.
  * [45] Ajay Dhar and Abhai Mansingh.  Polaronic hopping conduction in reduced lithium niobate single crystals.  Philosophical Magazine B, 61(6):1033–1042, 1990.


# On the effectiveness of Randomized Signatures as Reservoir for Learning
Rough Dynamics

1st Enea Monzio Compagnoni  Department of Mathematics
and Computer Science
University of Basel
Basel, Switzerland
<EMAIL_ADDRESS>2nd Anna Scampicchio  Institute for Dynamic
Systems and Control
ETH Zürich
Zürich, Switzerland
<EMAIL_ADDRESS>3rd Luca Biggio  Department of Computer Science
ETH Zürich
Zürich, Switzerland
<EMAIL_ADDRESS>4th Antonio Orvieto  Department of Computer Science
ETH Zürich
Zürich, Switzerland
<EMAIL_ADDRESS>5th Thomas Hofmann  Department of Computer
Science
ETH Zürich
Zürich, Switzerland
<EMAIL_ADDRESS>6th Josef Teichmann  Department of Mathematics
ETH Zürich
Zürich, Switzerland
<EMAIL_ADDRESS>Enea Monzio Compagnoni1, Anna Scampicchio2, Luca
Biggio3, Antonio Orvieto3,
Thomas Hofmann3 and Josef Teichmann4  1Department of Mathematics and Computer
Science, University of Basel, Basel, Switzerland,
Email<EMAIL_ADDRESS>2Institute for Dynamic Systems and
Control, ETH Zürich, Zürich, Switzerland. Email<EMAIL_ADDRESS>3Department of Computer Science, ETH Zürich, Zürich, Switzerland
Emails: $\\{$luca.biggio, antonio.orvieto<EMAIL_ADDRESS>4Department of Mathematics, ETH Zürich, Zürich, Switzerland. Email:
<EMAIL_ADDRESS>
###### Abstract

Many finance, physics, and engineering phenomena are modeled by continuous-
time dynamical systems driven by highly irregular (stochastic) inputs. A
powerful tool to perform time series analysis in this context is rooted in
rough path theory and leverages the so-called Signature Transform. This
algorithm enjoys strong theoretical guarantees but is hard to scale to high-
dimensional data. In this paper, we study a recently derived random projection
variant called Randomized Signature, obtained using the Johnson-Lindenstrauss
Lemma. We provide an in-depth experimental evaluation of the effectiveness of
the Randomized Signature approach, in an attempt to showcase the advantages of
this reservoir to the community. Specifically, we find that this method is
preferable to the truncated Signature approach and alternative deep learning
techniques in terms of model complexity, training time, accuracy, robustness,
and data hungriness.

###### Index Terms:

stochastic differential equations, reservoir computing, signature transform,
randomized signatures

##  I Introduction

We consider dynamical systems that are described by the stochastic
differential equation

$dY_{t}=f(Y_{t})dX_{t},\qquad Y_{0}=y_{0}\in\mathbb{R}^{m},$ (1)

where $f(\cdot):\mathbb{R}^{m}\rightarrow\mathbb{R}^{d}$ is an unknown smooth
map, and $X:[0,T]\rightarrow\mathbb{R}^{d}$ is a piece-wise smooth stochastic
process, also known as control, forcing the system evolution111For instance,
letting $d=2$, one can have $X_{t}=[t,\quad W_{t}]^{\top}$, where $W_{t}$ is a
$1-$dimensional Wiener process.. The problem under investigation consists in
predicting the solution $\bar{Y}_{t}$ of (1) given a new, unseen control
$\bar{X}_{t}$, using an algorithm trained on a set of control-output
trajectories. Two main challenges are the fact that data is observed on a
possibly unevenly spaced time grid, and that controls are often highly
irregular (we will mainly focus on rough paths, which are formally defined,
e.g., in [30], Chapter 3). This setting is of particular interest, e.g., in
high-frequency trading, as tick-level prices are not observed regularly in
time [14], and where fundamental quantities such as the volatility of prices
are rough [17].

A first strategy consists in performing system identification, i.e., learning
$f(\cdot)$ from data. The theory is well established for the particular case
of linear systems, both for parametric [27] and non-parametric approaches
leveraging the theory of Reproducing Kernel Hilbert Spaces [39, 1], and is
mainly deployed for discrete-time systems; its extension to the continuous-
time case is investigated, e.g., in [16] and references therein. Nonlinear
system identification is an active area of research, and the methods currently
deployed revolve around kernel techniques [38], sparse regression [4], random
features [40, 41], and deep neural networks (see, e.g., [33] for discrete-time
models, and [12] for continuous-time ones).

The advantage of estimating $f$ is that it allows for retrieving the output
trajectory also for different initial conditions and different time intervals.
However, solving the identification problem of a general continuous-time
system is typically hard if observations are made on an uneven time grid and
$f$ is nonlinear. Moreover, if the control signal is highly irregular, as
typical in finance and physical models [28, 25], integrating the differential
equation with the estimated $f$ is far from trivial. Thus, an alternative
viewpoint consists in focusing on estimating the solution $Y_{t}$ directly. To
this aim, deep neural networks have been successfully deployed [11, 15]: see
also the works on neural controlled differential equations [24, 35].
Nevertheless, their outstanding performance comes at the price of over-
parametrization, data hungriness, and expensive training cost [36, 32].
Furthermore, the resulting models learn representations of the input data that
are highly specialized to the training task. In addition, the remarkable
performance of these methods is often the result of a substantial engineering
effort and is not supported by theoretical results.

Another approach consists in reservoir computing [43], in which learning is
divided into two phases: first, data go through an untrained reservoir that
extracts a set of task-independent features; second, a simple and efficient-
to-train linear map (the readout map) projects such features into the desired
output. An example is Echo State Networks [23]. The critical point is that the
design of the reservoir determines the expressiveness of the features, and
several alternatives can be found in the literature (see [18] and references
therein).

A powerful reservoir is offered by the Signature Transform, often simply
referred to as Signature, stemming from rough path theory [20, 13]. The
Signature of a path is an infinite-dimensional tensor. Intuitively, it
consists of enhancing the path with additional curves corresponding to
iterated integrals of the curve with itself. A strong mathematical result
supports the choice of the Signature as a reservoir: it can be shown [26] that
the solution of a rough differential equation can be approximated arbitrarily
well by a linear map of the Signature of the controls. On the other hand, it
is often the case that this reservoir is very high-dimensional, and hence
particularly expensive to calculate and computationally intractable for use in
downstream tasks. Furthermore, the high dimensionality of the Signature poses
additional challenges for modern gradient-based optimizers, as convergence
rates suffer from a linear dependence in the model dimension [3].
Inspired by the remarkable theoretical properties of the Signature reservoir
and motivated to fix its practical pitfalls, the so-called Randomized
Signatures was introduced in [7, 8]. The Randomized Signature is obtained by
numerically integrating a set of random linear stochastic differential
equations driven by the control path. Importantly, based on a non-trivial
Johnson-Lindenstraus argument, [8] showed that calculating the Randomized
Signature of a path this way is equivalent to projecting its Signature using a
random linear operator. These random features provably retain the expressive
power of Signature, yet dramatically reduce the computational burden. In fact,
differently from [34], calculating the Randomized Signature does not require
computing the (truncated) Signature of the path: the projection can be
obtained directly in the compressed space. However, the lack of an in-depth
experimental study comparing its performance to Signature, Reservoir
Computing, and Deep Learning limits its popularity to the theoretical
community.

The contribution of the present paper is twofold: first, we extend the
theoretical analysis of Randomized Signature, using results from Malliavin
Calculus to prove that the Randomized Signature has the power of representing
the behavior of any dynamical system of interest; second, we provide a rich
set of experiments showing that this approach achieves performances comparable
with, if not superior to, competitive Deep Learning, Reservoir Computing, and
Signature-based models. In particular, we find that Randomized Signature
requires less trainable parameters, i.e. has lower model complexity, which in
turn implies a reduced training time and memory usage. In terms of
performance, our models are more accurate out-of-sample, more robust, and less
data-hungry – especially in high dimensions.

#### Notation

The canonical basis for $\mathbb{R}^{d}$ will be denoted as
$\\{e_{1},\dots,e_{d}\\}$. The symbol $\otimes$ represents a tensor product:
e.g., $e_{i}\otimes e_{j}$ is the $d\times d$ matrix of all zeros except for
the term at the $i-$th row and $j-$th columns, which equals 1. In general,
$\left(\mathbb{R}^{d}\right)^{\otimes l}$ is the space of tensors of shape
$(d,\ldots,d)$ given by $\mathbb{R}^{d}\otimes\cdots\otimes\mathbb{R}^{d}$ for
$l$ times. The tensor algebra on $\mathbb{R}^{d}$, and its truncated version
of order $M\geq 0$, are written as
$\mathcal{T}\left(\mathbb{R}^{d}\right):=\prod_{l=0}^{\infty}\left(\mathbb{R}^{d}\right)^{\otimes
l}$ and
$\mathcal{T}^{M}\left(\mathbb{R}^{d}\right):=\prod_{l=0}^{M}\left(\mathbb{R}^{d}\right)^{\otimes
l}$, respectively. Given two vector fields $V_{1}$ and $V_{2}$ mapping
$\mathbb{R}^{k}$ into itself, and denoting with $DV_{i}(z)$ the Fréchet
derivative of $V_{i}$ evaluated at $z\in\mathbb{R}^{k}$, the Lie bracket is
defined as $[V_{1},V_{2}](z)=DV_{1}(z)V_{2}(z)-DV_{2}(z)V_{1}(z)$.

##  II Background

Referring to the stochastic differential equation (1), we let the control
$X=\left(X^{1},\cdots,X^{d}\right):[0,T]\to\mathbb{R}^{d}$ be a continuous and
piece-wise smooth path – in particular, we will mainly regard $X$ as a rough
path222For the rigorous definition, we refer the reader to, e.g., [29, 28],
due to space limitations. For ease of visualization, one can think of $X$ as a
$d-$dimensional fractional Brownian motion with Hurst coefficient $H>1/4$
(Theorem D.3.2, [2]). on $\mathbb{R}^{d}$. We start by defining its Signature,
which is a tensor of iterated integrals of $X$ with itself.

######  Definition II.1 (Signature).

For any $t\in[0,T]$, the Signature of $X:[0,T]\to\mathbb{R}^{d}$ on $[0,t]$ is
the countable collection
$\mathbf{S}_{t}:=\left(1,S_{t}^{1},S_{t}^{2},\ldots\right)\in\mathcal{T}\left(\mathbb{R}^{d}\right)$
where, for each $l\geq 1$, the entries $S_{t}^{l}$ are defined as

$\small S_{t}^{l}:=\sum_{\begin{subarray}{c}\left(i_{1},\ldots,i_{l}\right)\\\
\in\\{1,\ldots,d\\}^{l}\end{subarray}}\left(\int_{0\leq s_{1}\leq\cdots\leq
s_{l}\leq t}dX^{i_{1}}_{s_{1}}\ldots
dX^{i_{l}}_{s_{l}}\right)e_{i_{1}}\otimes\dots\otimes e_{i_{l}}.$

Given that this object is infinite-dimensional, to actually compute the
reservoir, one can use only a finite amount of terms $S_{t}^{l}$. Therefore,
we consider the following object:

######  Definition II.2 (Truncated Signature).

The Truncated Signature of $X$ of order $M\geq 0$ is defined as

$\small\mathbf{S}^{M}_{t}:=\left(1,S_{t}^{1},\ldots,S_{t}^{M}\right)\in\mathcal{T}^{M}\left(\mathbb{R}^{d}\right).$
(2)

To give an intuition on how to compute the (truncated) Signature, we provide
the following

######  Example II.3.

Let $X:[0,T]\to\mathbb{R}$; then $S^{1}_{t}=\int_{0}^{t}dX_{s}$, which is
exactly $X_{t}-X_{0}$. To get $S^{2}_{t}$, we instead have to compute the
following iterated integral:
$S^{2}_{t}=\int_{0}^{t}\left(\int_{0}^{v}dX_{s}\right)dX_{v}$. Iterated
integrals of higher order $S^{j}_{t}$ are computed in a similar way, by
iteratively integrating the path $j$ times. As a practical example, let
$X_{t}=t$. Then it is easy to see that $S^{j}_{t}=\frac{t^{j}}{j!}$. Now let
$Y_{t}$ be an analytic function of time for which we have
$Y_{t}=\sum_{j=0}^{\infty}Y^{(j)}_{0}\frac{t^{j}}{j!}$. Taylor’s theorem
combined with the previous computation implies that $Y$ can be approximated as
a linear map of the Truncated Signature of $t$. Finally, note that $S^{j}_{t}$
gets smaller and smaller in magnitude as $j$ increases, which can be proven in
general. This suggests that the Truncated Signature can be safely used to
approximate $Y$. $\square$

The following result given in [8] proves that the solutions of differential
equations of the type given in (1) can be expanded in terms of the iterated
integrals stored in the Signature:

######  Theorem II.4 (Theorem 2.3, [8]).

Let $V_{i}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m},i=1,\ldots,d$ be vector
fields regular enough such that
$dY_{t}=\sum_{i=1}^{d}V^{i}\left(Y_{t}\right)dX^{i}_{t},Y_{0}=y\in\mathbb{R}^{m}$,
admits a unique solution $Y_{t}:[0,T]\rightarrow\mathbb{R}^{m}$. Then, for any
smooth test function $F:\mathbb{R}^{m}\rightarrow\mathbb{R}$ and for every
$M\geq 0$ there is a time-homogeneous linear operator
$L:\mathcal{T}^{M}\left(\mathbb{R}^{d}\right)\to\mathbb{R}$, which depends
only on $\left(V_{1},\ldots,V_{d},F,M,y\right)$, such that

$F\left(Y_{t}\right)=L\left(\mathbf{S}^{M}_{t}\right)+\mathcal{O}\left(t^{M+1}\right),\quad
t\in[0,T].$ (3)

Interpreting the linear operator $L$ as a readout layer, such a result
strongly motivates the use of (truncated) Signature as a valuable reservoir
under rough dynamics (note that smooth controls constitute a particular case).
The drawback of using $\mathbf{S}^{M}_{t}$ is that it has an
$\mathcal{O}(d^{M})$ computational complexity333Indeed, consider $d=2$:
${S}^{2}_{t}$ is a $2\times 2$ matrix with elements
$\int_{0}^{t}\left(\int_{0}^{v}dX^{1}_{s}\right)dX^{1}_{v}$,
$\int_{0}^{t}\left(\int_{0}^{v}dX^{1}_{s}\right)dX^{2}_{v}$,
$\int_{0}^{t}\left(\int_{0}^{v}dX^{2}_{s}\right)dX^{1}_{v}$ and
$\int_{0}^{t}\left(\int_{0}^{v}dX^{2}_{s}\right)dX^{2}_{v}$. For $M=3$, the
object to compute is instead a $2\times 2\times 2$ tensor, containing all
integrals of the type
$\int_{0}^{t}\left(\int_{0}^{w}\left(\int_{0}^{v}dX^{i_{1}}_{s}\right)dX^{i_{2}}_{v}\right)dX^{i_{3}}_{w}$
for all $i_{1},i_{2},i_{3}\in\\{1,2\\}$. Hence, the complexity — as well as
the dimensionality of the features — scales exponentially in $M$., which
becomes intractable for high-dimensional systems and/or for large values of
$M$ aimed at obtaining a finer representation of the solution. To cope with
this issue, instead of calculating the Signature, one can extract a new
quantity, called Randomized Signature, which is easier to compute and inherits
the expressiveness and inductive bias of the Signature.

######  Definition II.5 (Randomized Signature).

Given $k\in\mathbb{N}$ and random matrices $A_{1},\ldots,A_{d}$ in
$\mathbb{R}^{k\times k}$, random shifts $b_{1},\ldots,b_{d}$ in
$\mathbb{R}^{k\times 1}$, random starting point $z$ in $\mathbb{R}^{k}$, and
any fixed activation function $\sigma$, the Randomized Signature of $X$ in
$t\in[0,T]$ is the solution of the differential equation

$dZ_{t}=\sum_{i=1}^{d}\sigma\left(A_{i}Z_{t}+b_{i}\right)dX^{i}_{t},\quad
Z_{0}=z\in\mathbb{R}^{k}.$ (4)

The Randomized Signature is constructed in [8] as a random projection of the
Truncated Signature according to an argument based on the Johnson-
Lindenstrauss Lemma [45]. We refer the reader to [8] for all the details on
the theoretical derivation. The key message is the following:

######  Theorem II.6 (Informal).

For any number of features $k$ big enough, the Randomized Signature of $X$
defined in (4) can be linearly mapped to the solution of any differential
equation controlled by it, up to a small error vanishing at $k\to\infty$.

This result leads to the practical recipe for extracting the Randomized
Signature, summarized in Algorithm 1.

Algorithm 1 Generate Randomized Signature

0: $X\in\mathbb{R}^{d}$ sampled at $0=t_{0}<\dots<t_{N}=T$, dimension $k$ of
the Randomized Signature, and activation function $\sigma$.

Initialize: $Z_{0}\in\mathbb{R}^{k},A_{i}\in\mathbb{R}^{k\times
k},b_{i}\in\mathbb{R}^{k}$ to have i.i.d. standard normal entries.

for $n=1,\cdots,N$ do

$Z_{t_{n}}=Z_{t_{n-1}}+\sum_{i=1}^{d}\sigma\left(A_{i}Z_{t_{n-1}}+b_{i}\right)\left(X^{i}_{t_{n}}-X^{i}_{t_{n-1}}\right)$

end for

The computational complexity for calculating $Z$ is $\mathcal{O}(k^{2}d)$, and
its dimensionality is $\mathcal{O}(k)$. In Section V-B we show experimentally
that, in order to match the approximation capabilities of the Truncated
Signature of order $M$, the number $k$ of required Randomized Signatures is
fairly small – in particular, it is not exponential in $M$. This confirms that
working with Randomized Signatures is often less computationally demanding and
results in lower-dimensional – yet expressive – features.

##  III Randomized Signature as Reservoir: the procedure

Combining Theorems II.4 and II.6, one can perform linear (ridge) regression to
find the sought readout map. This is computed using observed (sampled)
control-output trajectories, and can then be used to predict the solution of
(1) given a new control sequence. The complete procedure for retrieving the
output sequence $\bar{Y}$ given a new control $\bar{X}$ is summarized in
Algorithm 2.

Algorithm 2 Simulate solution of (1)

0: Time grid $\mathcal{D}=\left\\{0=t_{0},\cdots,t_{N}=T\right\\}$;
$N_{train}$ input-output trajectories indexed by $m$,
$\\{(X_{t}(m),Y_{t}(m))\\}_{t\in\mathcal{D}}$, with common initial condition
$y_{0}\in\mathbb{R}^{m}$; new control $\\{\bar{X}_{t}\\}_{t\in\mathcal{D}}$;
order of Randomized Signature $k$; regularization parameter $\lambda$.

for $m=1,...,N_{train}$ do

compute the Randomized Signature $\\{Z_{t}(m)\\}_{t\in\mathcal{D}}$ via
Algorithm 1.

end for

Define $\mathbf{Y}\in\mathbb{R}^{(N+1)*N_{train}\times m}$ and
$\mathbf{Z}\in\mathbb{R}^{(N+1)*N_{train}\times k}$ such that
$\mathbf{Y}={\scriptsize\begin{bmatrix}Y_{t_{0}}(1)^{\top}\\\ \vdots\\\
Y_{t_{N}}(1)^{\top}\\\ \vdots\\\
Y_{t_{N}}(N_{train})^{\top}\end{bmatrix}},\quad\mathbf{Z}={\scriptsize\begin{bmatrix}Z_{t_{0}}(1)^{\top}\\\
\vdots\\\ Z_{t_{N}}(1)^{\top}\\\ \vdots\\\
Z_{t_{N}}(N_{train})^{\top}\end{bmatrix}}$

Solve  $\hat{\beta}=\arg\min_{\beta\in\mathbb{R}^{k\times
m}}\|\mathbf{Y}-\mathbf{Z}\beta\|^{2}+\lambda\|\beta\|^{2}.$

Compute the Randomized Signature of $\bar{X}$,
$\\{\bar{Z}_{t}\\}_{t\in\mathcal{D}}$, and store it in
$\mathbf{\bar{Z}}=[\bar{Z}_{t_{0}},\cdots,\bar{Z}_{t_{N}}]^{\top}$.

Retrieve
$\mathbf{\bar{Y}}=[\bar{Y}_{t_{0}},\cdots,\bar{Y}_{t_{N}}]=\mathbf{\bar{Z}}\hat{\beta}$.

Note that, while the choice of the activation function $\sigma$ does not
affect the theoretical results [8, 7], selecting it carefully positively
impacts expressiveness. Inspired by seminal works on the stability of deep
linear networks [19] and by the connection to neural ordinary differential
equations [6], it turns out that a good choice for $\sigma$ is a linear
function with $\frac{1}{d\times\sqrt{k}}$ as slope444The dynamics or
Randomized Signature is intrinsically exponential. This initialization
guarantees that the growth does not depend on the number of controls $d$ nor
on the number of features $k$.. The performance is further affected by the
Randomized Signature order $k$, and by the regularization parameter $\lambda$:
we select the first via cross-validation, and typically set the latter to the
value 0.001. Further investigation into these choices will be carried out in
future work.

##  IV Theoretical contribution

We now provide a novel insight into the expressive power of Randomized
Signature built as in Algorithm 1. Instead of relying on the theory of rough
paths or on approximation estimates, as done in the backbone of Theorem II.6,
we consider tools from Malliavin Calculus [31, 5].

######  Theorem IV.1.

Let us assume that $k\geq 2$, that $X$ is a $d$-dimensional Brownian motion,
and that the random matrices $A_{i}\in\mathbb{R}^{k\times k}$ and shifts
$b_{i}\in\mathbb{R}^{k}$ are independent and identically distributed following
a law absolutely continuous with respect to the Lebesgue measure on the space
of matrices. If the activation function $\sigma$ is real and analytic, then
the Randomized Signature $Z_{t}$ at time $t$ has a density with respect to
Lebesgue measure on $\mathbb{R}^{k}$ for almost all initial values
$Z_{0}=z\in\mathbb{R}^{k}$.

Proof. Considering the vector fields $V_{i}(z)=\sigma(A_{i}z+b_{i})$ for
$i=1,\cdots,d$, it holds that the Lie bracket $[V_{i},V_{j}](z)$ is
independent with respect to $V_{i}(z)$ and $V_{j}(z)$ almost surely: in fact,
independent random samples only meet with probability zero into zero sets of
non-constant analytic functions [9]. By an inductive argument, it follows that
the vector fields $z\to\sigma(A_{i}z+b_{i})$ satisfy Hörmander condition, i.e.
the Lie algebra generated by $\\{V_{i}(z)\\}_{i=1}^{d}$ spans
$\mathbb{R}^{k}$. The conclusion follows by applying, e.g., Theorem 7.4 in
[37]. $\blacksquare$
Theorem IV.1 shows that the process $Z_{t}$ will move in all directions with
positive probability: in other words, the obtained coordinate curves $t\to
Z_{t}^{i}$ form $k$ curves which are almost surely linearly independent in
time. As a further consequence, if the control $X_{t}$ is a Brownian motion,
for any partition $\mathcal{D}=\left\\{t_{0},\cdots,t_{N}\right\\}$ of $[0,T]$
of size $N+1$, also the sampling
$\left(Z^{i}_{t_{0}},\cdots,Z^{i}_{t_{N}}\right)$ are almost surely linearly
independent among each other if $k\geq N+1$. Therefore, for an appropriate
choice of $k\geq N+1$, Randomized Signature allows representing the behavior
of any target dynamical system on time grids $\mathcal{D}$. We highlight that
while this result is only proven when the control $X$ is a Brownian motion,
Theorem II.6 holds for any, possibly time-varying, rough path (control) X. In
the next Section, we show this experimentally.

##  V Numerical Experiments

We test the effectiveness of the Randomized Signature as a reservoir computer
in multiple challenging scenarios. We start by demonstrating the robustness of
our approach (Section V-A), as we show that the predictions of Algorithm 2
over multiple random initializations are consistent up to a negligible
deviation. Then, we display that it is an effective and efficient low-
dimensional compression of the Truncated Signature (Section V-B), and then we
show the resulting advantage in terms of data hungriness and computational
time (Section V-C). Next, Section V-D compares the performance of Randomized
Signatures against state-of-the-art techniques for simulation and system
identification methods in the presence of a control that is so irregular that
it does not even allow a formal definition of Signature. In Section V-E we
deepen such a comparison on the real-world scenario of an electrochemical
battery, where measurements are affected by noise. Next, we use the enzyme-
substrate model [22] to show the generalization property of the Randomized
Signature on out-of-distribution trajectories. Finally, we show on a scalar
Langevin equation with double-well potential that our proposed approach can
effectively deal with irregularly sampled time grids, which is a main
criticality in most of the state-of-the-art methods for trajectory prediction.

###  V-A Robustness over different random initializations

In this experiment, we show that the outputs of Algorithm 2 are stable across
different realizations of $A_{i}$, $b_{i}$, and $Z_{0}$. We consider the
Fractional Ornstein-Uhlenbeck process

$dY_{t}=\Theta\left(\mu-Y_{t}\right)dt+\Sigma dB^{(H)}_{t},\quad
Y_{0}=y_{0}\in\mathbb{R}^{m},$ (5)

where $B_{t}^{(H)}$ is an $m$-dimensional fractional Brownian motion of Hurst
parameter $H\in(0,1)$, $\mu\in\mathbb{R}^{m}$, and $\Theta$, $\Sigma$ are both
$m\times m$ positive semi-definite matrices. Relating this model to (1), we
have that $X_{t}=[t,\>(B_{t}^{(H)})^{\top}]^{\top}\in\mathbb{R}^{d}$ with
$d=m+1$, and $f(Y_{t})=[\Theta(\mu-
Y_{t}),\>\Sigma]\in\mathbb{R}^{m\times(m+1)}$. In this experiment, we take
$m=1$, $y_{0}=1$ and $\left(\mu=2,\Theta=1,\Sigma=2\right)$; we let $H=0.2$,
and the partition $\mathcal{D}$ of $[0,1]$ is made of $N=101$ equally spaced
times. For $10$ different random seeds, we draw different instances of
$A_{i}$, $b_{i}$, and $Z_{0}$, generate the reservoir $Z_{t}$ with $k=100$ and
apply Algorithm 2 to map $N_{\text{Train}}=100$ train samples of $Z$ into the
respective solution $Y_{t}$, to which we add white noise with variance $0.01$.
Figure 1 shows the average prediction ($\pm 3\times$ standard deviation) on a
test sample across the above-mentioned $10$ random seeds. Because the signal-
to-noise ratio is $\approx 50$, this shows that the model is robust to
different realizations of the Randomized Signature.

Figure 1: Experiment of Section V-A: average prediction with $\pm 3\times$
standard deviation bounds.

We conclude by highlighting that $k=100$ is a relatively low value with
respect to those used in the other experiments: we selected it to make the
error bars clearly visible.

######  Remark V.1.

We observed the same behavior consistently in all the other proposed
experiments, so we omit the Monte Carlo study in the next sections.

###  V-B Randomized Signature as compression of the Truncated one

Consider a $10$-dimensional control
$X_{t}=\left[t,(W_{t})^{\top}\right]^{\top}$ where $W_{t}$ is a
$9$-dimensional Brownian motion with independent components, and fix the order
of truncation of the Signature to $M=6$. Divide the time interval $[0,1]$
uniformly into $0=t_{0}<\cdots<t_{N}=1$ with $N=100$, and for each element in
the grid we compute both the Truncated Signature and the Randomized Signature
of order $k$, with $k$ taking values in $\\{1,\cdots,200\\}$. Reshaping the
two objects into matrices with dimension
$\left(N\times\left(\left(d^{(M+1)}-1\right)/(d-1)-1\right)\right)$ and
$N\times k$, respectively, we perform linear regression to find
$\mathbf{\beta}\in\mathbb{R}^{k\times\left(\left(d^{(M+1)}-1\right)/(d-1)-1\right)}$
mapping the Randomized Signature into the Truncated Signature. We observed
that, in order to obtain an approximation error of $10^{-4}$, we needed the
Randomized Signature to be of dimension approximately $k=190$. Therefore,
instead of calculating $\left(\left(d^{(M+1)}-1\right)/(d-1)-1\right)=1111110$
integrals per time step, we could just perform $k^{2}d=36100$ calculations per
time step, which is 3 times cheaper.

###  V-C Effectiveness of Randomized versus Truncated Signatures

Figure 2: Experiment of Section V-C: Randomized Signature vs. Truncated
Signature model. (Left) The number of trainable parameters for Randomized
Signature is significantly smaller regardless of the number of controls.
(Middle) Truncated Signature is much slower than Randomized Signature in high
dimensions. (Right) As opposed to Randomized Signature, the performance of
Truncated Signature degrades as the number of controls increases, and even
more so when the training set is small, thus indicating data hungriness.

We now deploy Truncated and Randomized Signatures to estimate the dynamics of
the Fractional Ornstein-Uhlenbeck process given in (5). In this experiment, we
fix $y_{0}=\mathbf{1}$, ${\mu}=\mathbf{1}$, ${\Sigma}=I_{d}$,
$[\Theta]_{i,j}=i/j$, the partition $\mathcal{D}$ of $[0,1]$ to have $N=101$
equally spaced time steps, and $H=0.3$. The order of truncation for the
Signature is set to $M=3$, and we consider different experiments with
increasing values of $m$ taking range in $\\{20,\cdots,80\\}$. To have a more
complete picture, we repeat the experiment in two cases, i.e. when the number
of training trajectories is $N_{train}=20$ and $N_{train}=50$. To keep the
computational cost of extracting features equal to $\mathcal{O}(d^{3})$ in
both the models given by Truncated and Randomized Signatures, we let $k=d$. As
a result, the number of features for the two models are $\mathcal{O}(d^{3})$
and $\mathcal{O}(d)$, respectively, strongly impacting the computational time
(middle panel of Figure 2). The right panel of Figure 2 shows also that the
performance of the Truncated Signature degenerates as the underlying
optimization problem explodes in dimension, while that of the Randomized one
is stable. This result clearly highlights the data hungriness of Signature-
based models when the number of dimensions is high.

###  V-D Comparison with baseline methods

In this experiment, we consider again the Fractional Ornstein-Uhlenbeck
process presented in (5) with $m=1$, $y_{0}=1$, and the same time grid
$\mathcal{D}$ of $101$ equally spaced points on [0,1], but we take a Hurst
coefficient $H=0.1$ corresponding to a highly irregular control. We benchmark
a Randomized Signature of order $k=50$ with the following: (a) Neural
Controlled Differential Equations (NCDEs) [24], a model which parametrizes the
vector fields of a latent controlled differential equation of dimension
$n_{latent}=100$ with feedforward neural networks with $1$ hidden layer of
$n_{nodes}=70$ nodes each, followed by a linear layer mapping the latent
variable into the output; (b) Echo State Networks (ESN) [23], which evolve the
input state according to an update rule which is that of an untrained
recurrent neural network that is ultimately linearly mapped into the output.
We chose the internal state to be of size $50$ (such that we have the same
number of trainable parameters as the model based on Randomized Signature) and
the activation functions to be hyperbolic tangents. The spectral radius and
leaking rate have been selected in cross-validation and set to $0.7$ and
$0.4$, respectively; (c) Neural Network Autoregressive model with Exogenous
Input (NNARX) as presented in [42], i.e. with $n_{a}=n_{b}=12$, $n_{k}=1$, and
using a feedforward neural network with input dimension $n_{a}+n_{b}+1=25$ and
2 hidden layers each with $100$ hidden units; (d) a Long Short-Term Memory
(LSTM) neural network [21] with $2$ hidden recursive layers of dimension $35$.
In this experiment, we use $N_{train}=1000$ trajectories to train the models,
and $N_{test}=1000$ to test the results. For the NNARX and NCDE models, we
minimized the mean square error optimizing with Adam with a learning rate of
$0.01$ for $100$ epochs. Similarly, for the LSTM model, we used Adam with a
learning rate of $0.001$ for $1000$ epochs. The results, showing the superior
performance of the Randomized Signature in terms of accuracy and computational
load, are presented in Table I.

| Average $L^{2}$ relative error | Training time [s] | # parameters
---|---|---|---
RS | $\mathbf{(1.02\pm 1.67)\cdot 10^{-5}}$ | $\mathbf{1.59}$ | $\mathbf{50}$
NCDE | $(7.57\pm 7.95)\cdot 10^{-2}$ | $3296.63$ | $14471$
ESN | $(4.24\pm 3.27)\cdot 10^{-2}$ | $3.01$ | $50$
NNARX | $(2.96\pm 8.33)\cdot 10^{-5}$ | $323.73$ | $12801$
LSTM | $(4.49\pm 6.95)\cdot 10^{-4}$ | $535.21$ | $15436$

TABLE I: Results for baseline comparison presented in Section V-D

###  V-E Real-world experiment: electrochemical battery model with noisy
observations

Figure 3: Experiment of Section V-E with electrochemical battery model. (Left)
Comparison between the ground truth and our prediction on a sample trajectory.
(Middle) Predictions on Test Sample for different values of $k$. Consistently,
the fit quality improves as $k$ increases. (Right) Comparison with NNARX.

In this experiment, we learn the dynamics of the electrochemical battery model
proposed in [10], which returns the voltage $Y$ as the current $X$ is injected
into the battery. This system is of real-world relevance, as it relies on
high-dimensional nonlinear physic-based differential equations that ensure the
high fidelity of the simulated data. We use the open-source NASA Prognostic
Model Package [44] to simulate voltage trajectories given input current
control paths. On a fixed equally spaced partition $\mathcal{D}$ of $[0,500]$,
we model the input current with step functions taking values $0$ or $1$ on
random sub-intervals of $[0,500]$. We apply Algorithm 2 to map
$N_{Train}=1000$ instances of $k$-dimensional Randomized Signature of the
controls into the respective solutions to which we add white noise with
variance 0.01. We consider $k=\\{50,166,1000\\}$. We compare our results with
NNARX – which, consistently with the experiment in Section V-D, is the best-
performing benchmark on this task. Specifically, we choose the parameters as
$n_{a}=n_{b}=12$ and $n_{k}=1$, which leads to the best results, and we use a
feedforward neural network with input dimension $n_{a}+n_{b}+1=25$ and 2
hidden layers each with either $22$ (NNARX22) or $1000$ (NNARX1000) hidden
units. We minimized the mean square error optimizing with Adam with a learning
rate of $0.01$ for $100$ epochs. The results are presented in Figure 3.
Furthermore, Table II presents the out-of-sample comparison in terms of Mean
Squared Error with respect to the ground truth averaged across $N_{test}=1000$
test trajectories. Note that the NNARX22 has around $1000$ trainable
parameters just like our model with $k=1000$ and that NNARX1000, which matches
our best model in terms of MSE, has around $10^{3}$ more trainable parameters.

| $k=50$ | $k=166$ | $k=1000$ | NNARX 22 | NNARX 1000
---|---|---|---|---|---
$N_{Test}=1000$ | $3.15\cdot 10^{-4}$ | $2.66\cdot 10^{-4}$ | $\mathbf{2.59\cdot 10^{-4}}$ | $3.88\cdot 10^{-4}$ | $2.62\cdot 10^{-4}$

TABLE II: Electrochemical battery model, experiment in Section V-E: MSE Error
Comparison

###  V-F Out-of-sample generalization on enzyme-substrate model

We consider the controlled differential equation describing the reaction
between concentrations of a substrate $S_{t}$ and of an enzyme
$E_{t}=1-S_{t}$, yielding the enzyme-substrate complex $C_{t}$ according to
the Michaelis-Menten model. Additional substrate is injected through a control
$X_{t}$, and the observed quantity of interest $Y_{t}$ is the chemical product
of the reaction – for instance, the latter can be glucose obtained from
lactose-lactase reaction. The overall kinetics can be described by the model

$\begin{cases}dS_{t}&=\left(k_{-1}C_{t}-k_{1}S_{t}\left(1-C_{t}\right)\right)dt+X_{t}dt\\\
dC_{t}&=-\left(k_{-1}C_{t}-k_{1}S_{t}\left(1-C_{t}\right)\right)dt-
k_{2}C_{t}dt\\\ dY_{t}&=k_{2}C_{t}dt.\end{cases}$

Following [22], we choose $(k_{1},k_{-1},k_{2})=(30,1,10)$, set
$\left(S_{0},C_{0},Y_{0}\right)=(0,0,0)$ and consider the evolution on
$t\in[0,1]$. We fix the time grid to have $N=101$ equally spaced time steps
and the control $X_{t}$ to follow the law of $W^{2}_{t}$ where $W_{t}$ is a
$1$-dimensional Brownian Motion (to ensure positivity). We apply Algorithm 2
to map $10^{5}$ instances of $222$-dimensional Randomized Signature $Z$ of the
controls into the respective solution $Y_{t}$. On the top of Figure 4, we plot
the comparison of the true and the generated time series on a test sample. As
we can see, the model has learned to correctly map a trajectory of $X_{t}$ to
the respective system response $Y_{t}$. More surprisingly, the bottom of such
a figure shows that our model is able to predict the correct output even if we
stimulate the system with a substrate injection that follows a completely
different law with respect to those used in training, i.e.
$X_{t}=0.5\cdot\mathbb{1}_{\left\\{W^{2}_{t}>0.5\right\\}}$. This suggests
that the system was correctly identified even out-of-distribution.

Figure 4: Experiment of Section V-F. Enzyme-Substrate Reactions stimulated
with: (Top) squared Brownian motion; (Bottom) step function, which is out-of-
distribution.

###  V-G Test on irregularly sampled grid

We consider the $1$-dimensional Langevin equation with double-well potential
given by

$dY_{t}=\theta Y_{t}\left(\mu-Y^{2}_{t}\right)dt+\sigma dW_{t},\quad
Y_{0}=y_{0}\in\mathbb{R},$ (6)

where $t\in[0,1]$, $W_{t}$ is a $1$-dimensional Brownian motion, and
$\left(\mu,\theta,\sigma\right)\in\mathbb{R}\times\mathbb{R}^{+}\times\mathbb{R}^{+}$;
in this experiment, we fix $y_{0}=1$ and
$\left(\mu=2,\theta=1,\sigma=1\right)$. For each train and test sample, the
partition $\mathcal{D}$ of $[0,1]$ is made of $N$ randomly drawn times. More
precisely, $\mathcal{D}=\\{0,t_{1},\cdots,t_{N-1},1\\}$ such that
$t_{k}=1/(1-\exp(-s_{k}))$ and $\\{s_{1},\cdots,s_{N-1}\\}$ are $N-2$
independent realizations of a uniform distribution $\mathcal{U}[0,1]$ sorted
in increasing order. As a result, the probability that two samples share the
same $\mathcal{D}$ is null. We apply Algorithm 2 with $N_{\text{Train}}=10000$
train samples, and Figure 5 shows the comparison on an out-of-sample generated
and true trajectory. Finally, Table III shows the Relative $L^{2}$ Error on
$10000$ test samples as we vary the number of time steps $N$ and $k$, and we
compare it to the respective experiment in case the time grid is regularly
spaced. As we can see, even though the performance is worse than the regularly
sampled setup, this technique proves to be anyway reliable on irregularly
sampled regimes.

Figure 5: Experiment of Section V-G with Langevin equation. Simulation of an
out-of-sample trajectory using an irregularly sampled time grid.

| $(N,k)=(11,111)$ | $(N,k)=(101,222)$ | $(N,k)=(1001,332)$
---|---|---|---
Irregular | $0.082735$ | $0.016885$ | $0.010902$
Regular | $0.026759$ | $0.004465$ | $0.003004$

TABLE III: Relative $L^{2}$ Error comparison with regularly and irregularly
sampled grid (Section V-G)

##  VI Conclusions

A challenging problem emerging in a plethora of fields consists in solving a
controlled stochastic differential equation. The main difficulties that can
arise in this situation may be: (a) the law governing the differential
equation is unknown, so one needs to rely on sampled input/output
trajectories; (b) the samples are observed on an irregular time grid; (c) the
input trajectory is highly irregular, e.g., is a rough path. To cope with
them, this work investigated the power of Randomized Signature as a reservoir.
Such an approach proved to be very effective in estimating the solution of the
stochastic differential equation driven by a new control input, showing its
low data hungriness and robustness compared with state-of-the-art system
identification and deep learning-based methods. Further investigations will
aim at providing deeper theoretical results on the generalization capability
of Randomized Signature.

## References

  * [1] Nachman Aronszajn.  Theory of reproducing kernels.  Transactions of the American Mathematical Society, 68:337–404, 1950\.
  * [2] Francesca Biagini, Yaozhong Hu, Bernt Øksendal, and Tusheng Zhang.  Stochastic Calculus for Fractional Brownian Motion and Applications.  Springer, 2008.
  * [3] Léon Bottou, Frank E. Curtis, and Jorge Nocedal.  Optimization methods for large-scale machine learning.  Siam Review, 60(2):223–311, 2018.
  * [4] Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the National Academy of Sciences, 113(15):3932–3937, 2016.
  * [5] Thomas Cass and Peter K. Friz.  Densities for rough differential equations under hormander’s condition.  Annals of Mathematics, 171:2115–2141, 2007.
  * [6] Ricky T.Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud.  Neural ordinary differential equations.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 6572–6583, 2018.
  * [7] Christa Cuchiero, Lukas Gonon, Lyudmila Grigoryeva, Juan-Pablo Ortega, and Josef Teichmann.  Discrete-time signatures and randomness in reservoir computing.  IEEE Transactions on Neural Networks and Learning Systems, 2021\.
  * [8] Christa Cuchiero, Lukas Gonon, Lyudmila Grigoryeva, Juan-Pablo Ortega, and Josef Teichmann.  Expressive power of randomized signature.  In The Symbiosis of Deep Learning and Differential Equations, 2021\.
  * [9] Christa Cuchiero, Martin Larsson, and Josef Teichmann.  Deep neural networks, generic universal interpolation, and controlled ODEs.  SIAM Journal on Mathematics of Data Science, 2(3):901–919, 2020\.
  * [10] Matthew Daigle and Chetan S. Kulkarni.  Electrochemistry-based battery modeling for prognostics.  In Annual Conference of the PHM Society, volume 5, 2013.
  * [11] Hassan Ismail Fawaz, Germain Forestier, Jonathan Weber, Lhassane Idoumghar, and Pierre-Alain Muller.  Deep learning for time series classification: a review.  Data mining and knowledge discovery, 33(4):917–963, 2019.
  * [12] Marco Forgione and Dario Piga.  Continuous-time system identification with neural networks: Model structures and fitting criteria.  European Journal of Control, 59:69–81, 2021.
  * [13] Peter K. Friz and Martin Hairer.  A course on rough paths.  Springer, 2020.
  * [14] Tak-chung Fu.  A review on time series data mining.  Engineering Applications of Artificial Intelligence, 24(1):164–181, 2011.
  * [15] John Cristian Borges Gamboa.  Deep learning for time-series analysis.  arXiv preprint arXiv:1701.01887, 2017.
  * [16] Hugues Garnier and Peter C. Young.  The advantages of directly identifying continuous-time transfer function models in practical applications.  International Journal of Control, 87(7):1319–1338, 2014.
  * [17] Jim Gatheral, Thibault Jaisson, and Mathieu Rosenbaum.  Volatility is rough.  Quantitative finance, 18(6):933–949, 2018.
  * [18] Daniel J. Gauthier, Erik Bollt, Aaron Griffith, and Wendson A.S. Barbosa.  Next generation reservoir computing.  Nature Communications, 12(1), Sep 2021.
  * [19] Xavier Glorot and Yoshua Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In Proceedings of the thirteenth international conference on artificial intelligence and statistics, pages 249–256. JMLR Workshop and Conference Proceedings, 2010.
  * [20] Ben Hambly and Terry J. Lyons.  Uniqueness for the signature of a path of bounded variation and the reduced path group.  Annals of Mathematics, 171:109–167, 2010.
  * [21] Sepp Hochreiter and Jürgen Schmidhuber.  Long Short-Term Memory.  Neural Computation, 9(8):1735–1780, 11 1997.
  * [22] Brian P. Ingalls.  Mathematical modeling in systems biology: an introduction.  MIT press, 2013.
  * [23] Herbert Jaeger.  Adaptive nonlinear system identification with echo state networks.  NIPS, 06 2003.
  * [24] Patrick Kidger, James Morrill, James Foster, and Terry J. Lyons.  Neural controlled differential equations for irregular time series, 2020\.
  * [25] Peter E. Kloeden and Eckhard Platen.  Numerical Solution of Stochastic Differential Equations.  Stochastic Modelling and Applied Probability. Springer Berlin Heidelberg, 2013.
  * [26] Daniel Levin, Terry J. Lyons, and Hao Ni.  Learning from the past, predicting the statistics for the future, learning an evolving system.  arXiv preprint arXiv:1309.0260, 2013.
  * [27] Lennart Ljung.  System Identification, Theory for the User.  Prentice Hall, 1997.
  * [28] Terry J. Lyons.  Rough paths, signatures and the modelling of functions on streams.  arXiv preprint arXiv:1405.4537, 2014.
  * [29] Terry J. Lyons, Michael Caruana, and Thierry Lévy.  Differential Equations Driven by Rough Paths.  Springer, 2004.
  * [30] Terry J. Lyons and Zhongmin Qian.  System Control and Rough Paths.  Oxford University Press, 12 2002.
  * [31] P. Malliavin.  Stochastic Analysis.  Grundlehren der mathematischen Wissenschaften. Springer Berlin Heidelberg, 1998.
  * [32] Gary Marcus.  Deep learning: A critical appraisal.  arXiv preprint arXiv:1801.00631, 2018.
  * [33] Daniele Masti and Alberto Bemporad.  Learning nonlinear state-space models using deep autoencoders.  In 2018 IEEE Conference on Decision and Control (CDC), pages 3862–3867. IEEE, 2018.
  * [34] James Morrill, Adeline Fermanian, Patrick Kidger, and Terry J. Lyons.  A generalised signature method for multivariate time series feature extraction.  arXiv preprint arXiv:2006.00873, 2020.
  * [35] James Morrill, Patrick Kidger, Lingyi Yang, and Terry J. Lyons.  Neural controlled differential equations for online prediction tasks.  arXiv preprint arXiv:2106.11028, 2021.
  * [36] Behnam Neyshabur, Zhiyuan Li, Srinadh Bhojanapalli, Yann LeCun, and Nathan Srebro.  The role of over-parametrization in generalization of neural networks.  In International Conference on Learning Representations, 2018.
  * [37] David Nualart.  Malliavin Calculus and Its Applications.  Regional conference series in mathematics. Conference Board of the Mathematical Sciences, 2009.
  * [38] Gianluigi Pillonetto.  System identification using kernel-based regularization: New insights on stability and consistency issues.  Automatica, 93:321–332, 2018.
  * [39] Gianluigi Pillonetto, Tianshi Chen, Alessandro Chiuso, Giuseppe De Nicolao, and Lennart Ljung.  Regularized system identification - Learning dynamic models from data.  Communications and Control Engineering. Springer Cham, 2022.
  * [40] Ali Rahimi and Benjamin Recht.  Random features for large-scale kernel machines.  In J. Platt, D. Koller, Y. Singer, and S. Roweis, editors, Advances in Neural Information Processing Systems, volume 20. Curran Associates, Inc., 2008.
  * [41] Alessandro Rudi and Lorenzo Rosasco.  Generalization properties of learning with random features.  In NIPS, pages 3215–3225, 2017.
  * [42] Johan Schoukens and Lennart Ljung.  Nonlinear system identification: A user-oriented road map.  IEEE Control Systems Magazine, 39(6):28–99, 2019.
  * [43] Benjamin Schrauwen, David Verstraeten, and Jan Van Campenhout.  An overview of reservoir computing: theory, applications and implementations.  In Proceedings of the 15th european symposium on artificial neural networks. p. 471-482 2007, pages 471–482, 2007.
  * [44] Christopher Teubert, Matteo Corbetta, Chetan Kulkarni, and Matthew Daigle.  Prognostics models python package, August 2021.
  * [45] Santosh S. Vempala.  The Random Projection Method.  DIMACS Series. American Mathematical Society, 2005.