# Stability of Exponentially Damped Oscillations under Perturbations of the
Mori-Chain

Robin Heveling<EMAIL_ADDRESS>Department of Physics, University of
Osnabrück, D-49076 Osnabrück, Germany  Jiaozi Wang<EMAIL_ADDRESS>Department
of Physics, University of Osnabrück, D-49076 Osnabrück, Germany  Christian
Bartsch<EMAIL_ADDRESS>Department of Physics, University of Osnabrück,
D-49076 Osnabrück, Germany  Jochen Gemmer<EMAIL_ADDRESS>Department of
Physics, University of Osnabrück, D-49076 Osnabrück, Germany

###### Abstract

There is an abundance of evidence that some relaxation dynamics, e.g.,
exponential decays, are much more common in nature than others. Recently,
there have been attempts to trace this dominance back to a certain stability
of the prevalent dynamics versus generic Hamiltonian perturbations. In the
paper at hand, we tackle this stability issue from yet another angle, namely
in the framework of the recursion method. We investigate the behavior of
various relaxation dynamics with respect to alterations of the so-called
Lanczos coefficients. All considered scenarios are set up in order to comply
with the “universal operator growth hypothesis”. Our numerical experiments
suggest the existence of stability in a larger class of relaxation dynamics
consisting of exponentially damped oscillations. Further, we propose a
criterion to identify “pathological” perturbations that lead to uncommon
dynamics.

##  I Introduction

The apparent emergence of irreversibility from the underlying reversible
theory of quantum mechanics is a long-standing puzzle that lacks an entirely
satisfying answer to this day Gogolin and Eisert (2016). Over the course of
the last decades, fundamental concepts like the “eigenstate thermalization
hypothesis” Srednicki (1994); Deutsch (1991); D’Alessio  _et al._ (2016) and
“quantum typicality” Lloyd (2013); Goldstein  _et al._ (2006); Reimann (2007)
have crystallized, which give conditions under which isolated quantum systems
eventually reach an equilibrium state. However, while these mechanisms ensure
eventual equilibration, they make no statement in which manner the equilibrium
state is actually reached, i.e., they do not narrow down the eligible routes
to equilibrium. In contrast, it is evidently true that some relaxation
dynamics, e.g., exponential decays, are much more common in nature than
others, e.g., recurrence dynamics.
There are recent attempts to explain this prevalence with the idea that some
dynamics are stable versus perturbations of the Hamiltonian. These attempts
include, for example, investigations based on Hamiltonian perturbations on the
level of random matrices Knipschild and Gemmer (2018). Further, there exist
advances that consider an entire ensemble of permissible perturbations Dabelow
and Reimann (2020); Dabelow  _et al._ (2020); Dabelow and Reimann (2021). It
is then analytically shown that weak “typical” perturbations lead to an
exponential damping of the original dynamics Dabelow and Reimann (2020). This
renders exponential decays stable since only the decay constant changes,
whereas recurrence dynamics are exponentially suppressed.
With the paper at hand, we address the issue of stability of certain [classes
of] relaxation dynamics in the framework of the recursion method Viswanath and
Müller (2008); Mattis (1981). Central quantities that appear within this
framework are the so-called Lanczos coefficients, real numbers that
characterize the complexity growth of operators over the course of time. As
will be presented in Sec. II, the Lanczos coefficients can be interpreted as
hopping amplitudes in a tight-binding model. In this manner, many physical
problems, like calculating correlation functions, can practically be reduced
to a one-dimensional [finite or semi-infinite] chain, which we refer to as the
“Mori-chain” in the title of this paper. In the following, we consider
perturbations on the level of Lanczos coefficients and examine the resulting
effect on various kinds of relaxation dynamics.
The particular advantage of our approach within the recursion method framework
is that all considered scenarios can be set up to directly comply with the
universal operator growth hypothesis Parker  _et al._ (2019), which previous
approaches have been lacking Knipschild and Gemmer (2018). Said hypothesis
concerns the asymptotic growth of the Lanczos coefficients and it basically
states that the coefficients should eventually attain linear growth [with a
logarithmic correction in one dimension]. The hypothesis is backed up by
analytical as well as numerical evidence Parker  _et al._ (2019); Noh (2021);
Heveling  _et al._ (2022).
The paper at hand is structured as follows: Sec. II constitutes a preliminary
section on the recursion method and the operator growth hypothesis.
Afterwards, in Sec. III, the concrete strategy to study the stability of
certain classes of dynamics is explained in detail. In Sec. IV we present and
discuss our numerical results. We conclude in Sec. V.

##  II Preliminaries: Recursion method and operator growth hypothesis

In this section, we briefly recall the basics of the recursion method
Viswanath and Müller (2008); Mattis (1981) as well as the universal operator
growth hypothesis Parker  _et al._ (2019). We consider a system described by
some Hamiltonian $\mathcal{H}$. An observable of interest represented by a
Hermitian operator $\mathcal{O}$ gives rise to a corresponding autocorrelation
function

$C(t)=\text{Tr}[\mathcal{O}(t)\mathcal{O}]\,,$ (1)

where
$\mathcal{O}(t)=e^{\mathrm{i}\mathcal{H}t}\mathcal{O}e^{-\mathrm{i}\mathcal{H}t}$
is the time-dependent operator in the Heisenberg picture ($\hbar=1$). In the
following, it is convenient to work directly in Liouville space, i.e., the
Hilbert space of operators, and denote its elements $\mathcal{O}$ as states
$|\mathcal{O})$. Elements of the Liouville space evolve under the Liovillian
$\mathcal{L}=[\mathcal{H},\cdot\,]$, i.e.,
$|\mathcal{O}(t))=e^{\mathrm{i}\mathcal{L}t}|\mathcal{O})$, similar to wave
functions that evolve under the Hamiltonian $\mathcal{H}$. Using the
Liovillian superoperator, the autocorrelation function may be written as
$C(t)=(\mathcal{O}|e^{\mathrm{i}\mathcal{L}t}|\mathcal{O})$. The Liouville
space is equipped with an infinite-temperature inner product
$(\mathcal{O}_{1}|\mathcal{O}_{2})=\text{Tr}[\mathcal{O}_{1}^{\dagger}\mathcal{O}_{2}]$,
which induces a norm via $||\mathcal{O}||=\sqrt{(\mathcal{O}|\mathcal{O})}$.
In the following, the central object of interest is the Liouvillian
$\mathcal{L}$ represented in a particular basis $\\{|\mathcal{O}_{n})\\}$, the
so-called Krylov basis. The Krylov basis is routinely constructed as part of
the Lanczos algorithm. In this basis, which is determined by some “seed”
observable $\mathcal{O}$, the representation of the Liouvillian is
tridiagonal. To initialize the algorithm, we take the normalized state
$|\mathcal{O}_{0})=|\mathcal{O})$, i.e., $(\mathcal{O}|\mathcal{O})=1$, and
set $b_{1}=||\mathcal{L}\mathcal{O}_{0}||$ as well as
$|\mathcal{O}_{1})=\mathcal{L}|\mathcal{O}_{0})/b_{1}$. Then we iteratively
compute

$\displaystyle|\mathcal{Q}_{n})$
$\displaystyle=\mathcal{L}|\mathcal{O}_{n-1})-b_{n-1}|\mathcal{O}_{n-2})\,,$
(2) $\displaystyle b_{n}$ $\displaystyle=||\mathcal{Q}_{n}||\,,$
$\displaystyle|\mathcal{O}_{n})$ $\displaystyle=|\mathcal{Q}_{n})/b_{n}\,.$

The tridiagonal representation of the Liouvillian in the Krylov basis
$\\{|\mathcal{O}_{n})\\}$ is then given by

$L_{mn}=(\mathcal{O}_{m}|\mathcal{L}|\mathcal{O}_{n})=\begin{pmatrix}0&b_{1}\vphantom{\vdots}&0&...\\\
\vphantom{...}b_{1}&0\vphantom{\vdots}&b_{2}&\\\ 0&b_{2}&0&\ddots\\\
\vdots&&\ddots&\ddots\end{pmatrix}_{mn}\,,$ (3)

where the Lanczos coefficients $b_{n}$ are real, positive numbers output by
the algorithm. For our purposes, it is sufficient to assume a finite-
dimensional space. Then, the algorithm halts at step $n=d+1$, where $d$ is the
dimension of the Liouville space, and $L$ is a $d\times d$-matrix. Rewriting
the Heisenberg equation of motion in the Krylov basis yields

$\partial_{t}\varphi=-b_{n+1}\varphi_{n+1}+b_{n}\varphi_{n}\,,$ (4)

where we defined
$\varphi_{n}:=\mathrm{i}^{-n}(\mathcal{O}_{n}|\mathcal{O}(t))$. The initial
condition is given by $\varphi_{n}(0)=\delta_{n0}$ and both $\varphi_{n}$ and
$b_{n}$ are set to zero by convention when $n$ is “out of bounds”. The above
Eq. (4) takes the form of a discrete Schrödinger equation and can be
numerically solved by familiar means of, e.g., exact diagonalization or
iterative schemes like Runge-Kutta and Chebyshev polynomials.
As mentioned, the Lanczos coefficients $b_{n}$ can be interpreted as hopping
amplitudes in a tight-binding model. Then, the correlation function coincides
with the amplitude of the first site, i.e., $C(t)=\varphi_{0}(t)$.
For later reference, we introduce the spectral function $\Phi(\omega)$ as the
Fourier transform of the correlation function, i.e.,

$\Phi(\omega)=\int_{-\infty}^{\infty}e^{-\mathrm{i}\omega
t}\,C(t)\,\text{d}t\,.$ (5)

It can be shown that the Lanczos coefficients $b_{n}$ appear in the continued
fraction expansion of $\Phi(\omega)$.

$\Phi(\omega)=\text{Re
}\dfrac{2}{\mathrm{i}\omega+\dfrac{b_{1}^{2}}{\mathrm{i}\omega+\dfrac{b_{2}^{2}}{\mathrm{i}\omega+...}}}$
(6)

Consequently, there exists a (non-linear) one-to-one map between the Lanczos
coefficients $b_{n}$ and the autocorrelation function $C(t)$. Thus, a set of
$b_{n}$’s uniquely determines $C(t)$ and vice versa.
Lastly, we present the universal operator growth hypothesis as brought forth
in Ref. Parker  _et al._ (2019). The hypothesis concerns the asymptotic
behavior of the Lanczos coefficients $b_{n}$ and basically states that in
generic, non-integrable systems the Lanczos coefficients of local, few-body
observables grow asymptotically linear, i.e., above some $n$ the growth is
given by

$b_{n}\sim\alpha n+\gamma+o(1)\,,$ (7)

where $\alpha>0$ and $\gamma$ are real constants and $o(g_{n})$ denotes some
real sequence $f_{n}$ with $\lim_{n\rightarrow\infty}|f_{n}/g_{n}|=0$. In the
special case of a one-dimensional system, the asymptotic growth is sub-linear
due to an additional logarithmic correction Parker  _et al._ (2019).

##  III Probing for Stability

In this section, we devise our strategy to probe the stability of classes of
dynamics with respect to certain types of perturbations. This section consists
of the following parts: We present ways to find Lanczos coefficients
corresponding to a chosen [class of] correlation functions [Sec. III.1]. Next,
we present four physically motivated requirements that the Lanczos
coefficients should fulfill [Sec. III.2]. Then, we discuss the particular type
of perturbations we consider [Sec. III.3]. Finally, we introduce quantifiers
that allow to evaluate the stability of the respective dynamics [Sec. III.4].

###  III.1 Dynamics and Lanczos coefficients

The preliminary section already touched on the relation between a correlation
function $C(t)$ and the corresponding Lanczos coefficients $b_{n}$. As
mentioned, this correspondence is one-to-one, thus, a set of $b_{n}$’s
uniquely determines $C(t)$ and vice versa. However, the convergence can be
quite subtle, i.e., there may be fairly similar dynamics with vastly different
Lanczos coefficients. On the other side, similar Lanczos coefficients may lead
to quite different dynamics.
In the numerical experiments below, we proceed as follows. First, we choose a
correlation function $C^{\prime}(t)$ whose stability we want to probe. We may
specify $C^{\prime}(t)$ as a concrete analytical function, or as a member of a
class of functions, e.g., exponential decays. The exact Lanczos coefficients
corresponding to $C^{\prime}(t)$ are denoted by $b^{\prime}_{n}$. These
coefficients are, in general, unknown and obtaining them is quite a difficult
task. Our objective is to find Lanczos coefficients $b_{n}$ corresponding to a
correlation function $C(t)$ that should be practically indistinguishable from
$C^{\prime}(t)$ for all intents and purposes [$C(t)\simeq C^{\prime}(t)$].
These “approximate” coefficients $b_{n}$ may be obtained in three different
ways:
i. The exact coefficients $b^{\prime}_{n}$ may be analytically known. In this
case $b_{n}=b^{\prime}_{n}$ as well as $C(t)=C^{\prime}(t)$.
ii. The coefficients may be obtained via an educated guess. In this case, the
$b_{n}$ may be quite different from $b^{\prime}_{n}$ while still $C(t)\simeq
C^{\prime}(t)$.
iii. The coefficients may be “reverse-engineered” from the correlation
function $C^{\prime}(t)$. For details see App. A.
In the following, we will omit the technical distinction between $C(t)$ and
$C^{\prime}(t)$, since we assume that Lanczos coefficients can be found for
which the two correlation functions are practically indistinguishable.

###  III.2 Design of the Lanczos coefficients

In the following, we compare dynamics of correlation functions $C_{A}(t)$ and
$C_{B}(t)$ with respect to their stability under a certain class of
perturbations. Corresponding sets of Lanczos coefficients $b^{A}_{n}$ and
$b^{B}_{n}$ may be obtained by one of the three methods mentioned in the last
section. We demand that these Lanczos coefficients satisfy the following
[partly physically motivated] requirements sufficiently well:
i. The Lanczos coefficients $b^{X}_{n}$ should reproduce the respective
correlation function $C_{X}(t)$ to an acceptable accuracy, where $X=A,B$ (this
is sort of obvious and has already been addressed above).
ii. The Lanczos coefficients $b^{X}_{n}$ should comply with the universal
operator growth hypothesis, i.e., they should eventually attain linear growth.
iii. The resulting correlation functions $C_{A}(t)$ and $C_{B}(t)$ should
decay on more or less the same time scale. This is to ensure a fair comparison
between the two, since faster dynamics are typically less affected by
perturbations than slower ones.
iv. The Lanczos coefficients $b^{A}_{n}$ and $b^{B}_{n}$ should be similar in
magnitude. In particular, we demand that the quantity
$\sum_{n}(b^{X}_{n})^{2}$ is comparable for $X=A,B$. We show in App. B that
this sum is related to the spectral variance of the Hamiltonian. Hence, this
condition fosters the notion that the two correlation functions $C_{A}(t)$ and
$C_{B}(t)$ originate from different observables while the underlying
Hamiltonians are quite similar. In practice, a value of

$q(b^{\text{A}}_{n},b^{\text{B}}_{n})=\dfrac{\sum_{n}(b^{\text{A}}_{n})^{2}}{\sum_{n}(b^{\text{B}}_{n})^{2}}$
(8)

close to unity is desirable.

###  III.3 Design of the perturbations

The considered perturbations are designed on the level of Lanczos
coefficients. The coefficients $b_{n}$ corresponding to some correlation
function $C(t)$ will be slightly altered according to

$\tilde{b}_{n}=b_{n}+\lambda v_{n}\,,$ (9)

where the perturbed coefficients are denoted by $\tilde{b}_{n}$. Here,
$\lambda$ is the perturbation strength and $v_{n}$ is specified below. We show
in App. C that this particular form of the perturbation yields a sensible
scaling of the perturbed Hamiltonian with $\lambda$.
In general, it is an intricate problem to determine how a perturbation in the
form of $v_{n}$ corresponds to a perturbation $\mathcal{V}$ on the level of
Hermitian matrices. This is an issue that we do not attempt to tackle in this
paper.
Here, we try to make as few assumptions as possible and model the perturbation
$v_{n}$ as random numbers [with the only restriction of a minimal correlation
length, see below]. Concretely, we set

$v_{n}=\sum_{k=1}^{N_{\text{f}}}x_{k}\cos(2\pi nk/d)+y_{k}\sin(2\pi nk/d)\,,$
(10)

where the $x_{k},y_{k}$ are real, random numbers from a Gaussian distribution
with zero mean and unit variance. They are normalized as
$\sum_{k}x_{k}^{2}+y_{k}^{2}=1$.
The sum in Eq. (10) is capped at a number $N_{\text{f}}$. This corresponds to
excluding shorter wavelength or higher frequencies, which induces a minimal
correlation length in the coefficients $\tilde{b}_{n}$ [reddish noise]. No
truncation at all, i.e., $N_{\text{f}}=d$, would be equivalent to adding
uncorrelated random numbers [white noise] to the coefficients $b_{n}$. This
choice can be justified by invoking the tight-binding model interpretation in
which the Lanczos coefficients represent hopping amplitudes between
neighboring sites. Simply altering the hopping amplitudes in a random,
uncorrelated manner leads to localization effects similar to Anderson
localization Anderson (1958). Through unsystematic varying of $N_{\text{f}}$,
we find that the transition from unlocalized to to localized behavior is quite
sharp. Since the aim is to study the stability of certain relaxation dynamics,
we choose $N_{\text{f}}$ as large as possible, but small enough to avoid said
localization effects. In practice, this amounts to $N_{\text{f}}\approx d/3$.
In Sec. IV.3, we ease this restriction to identify “pathological”
perturbations.

###  III.4 Assessing stability

Once the perturbation has acted on the coefficients $b_{n}$ and yielded the
perturbed coefficients $\tilde{b}_{n}$, we can calculate the corresponding
perturbed dynamics $\tilde{C}(t)$ as laid out in Sec. II. Now, some tools are
needed in order to judge how strongly the dynamics were altered by the
perturbation. In particular, it is important to assess to what extent the
perturbed dynamics still falls into the original class of functions to which
the unperturbed dynamics belonged. For example, we may ask if $\tilde{C}(t)$
is still exponential, given that the original dynamics $C(t)$ was exponential
[the decay constants may differ]. To this end, we attempt to fit the perturbed
dynamics with a function that also described the unperturbed dynamics, e.g.,
we may try to fit $\tilde{C}(t)$ with $f(t)=Ae^{-\mu t}$. In practice, the fit
is obtained via a standard fitting routine.
Before we introduce quantifiers that measure the effect of the perturbation,
some remarks on the fitting ansatz are in order. Correlation functions are
symmetric in time, thus, all odd moments necessarily vanish. In particular,
correlation functions have zero slope at $t=0$. Further, we normalize the
correlation function to unity at $t=0$. However, for fitting we choose
functions that lack these features, i.e., the above exponential ansatz may
yield an $A$ that slightly differs from unity. Further, the slope of an
exponential at $t=0$ is never zero. This “negligence” is due to the fact that
we are more interested in an overall description of relaxation dynamics,
rather than in the details of the short-time behavior.
We assess the quality of a fit $f(t)$ by calculating “how far off” it is from
the given perturbed dynamics. Concretely, we define a measure of the “error”
or “deviation” $\epsilon$ by the expression

$\epsilon=\sqrt{\dfrac{1}{N_{\text{eq}}}\sum_{n=0}^{N_{\text{eq}}}(\tilde{C}(t_{n})-f(t_{n}))^{2}}\,.$
(11)

Here, the respective functions are evaluated at available points in time
$t_{n}=n\delta t$, where $\delta t$ is the time step used to solve the
equation of motion. The upper bound $N_{\text{eq}}$ corresponds to a time at
which the dynamics in question has seemingly equilibrated not .
To get a feeling of which numerical value of $\epsilon$ constitutes a good or
bad fit, we refer to the exemplary numerical data displayed in Figs. 3, 4, 8
and 9.
We introduce a second quantifier $\sigma$ that measures how strongly the
unperturbed dynamics are altered due to the perturbation in the first place.

$\sigma=\sqrt{\dfrac{1}{N_{\text{eq}}}\sum_{n=0}^{N_{\text{eq}}}(\tilde{C}(t_{n})-C(t_{n}))^{2}}\,.$
(12)

The construction of this quantity is similar to the construction of the
quantity $\epsilon$ that measures the fit quality.

##  IV Numerical analysis

In this section, we apply the presented strategy to certain relaxation
dynamics. In particular, we investigate and compare the stability of two
classes of dynamics. The first class in question is the class of damped
oscillations, whose damping is due to an ordinary exponential factor. This
class also includes simple exponential decays, which can be viewed as damped
oscillations with zero frequency. These dynamics are ubiquitous in nature. For
example, slow exponential dynamics may commonly arise whenever a system
interacts weakly with an environment or whenever long-wavelength Fourier
components of spatial densities of conserved quantities are considered.
Further, exponentially damped oscillations are routinely observed as short-
wavelength components.
The possible choices for the second class of dynamics are of course manifold.
Since an exhaustive analysis is impossible, we decide on a “Gaussianized”
version of the first class. This is supposed to mean that the class includes
Gaussian decays as well as oscillations damped by Gaussian factors.
The restriction of the second class of course stifles any aspiration of
generality. Thus, the following results should be viewed as mere numerical
experiments that corroborate the existence of stability in the first class.
Concretely, in Sec. IV.1, we compare exponential and Gaussian decay. In Sec.
IV.2, the comparison is between two damped oscillations, one with an
exponential damping factor and the other with a Gaussian damping factor.
Lastly, in Sec. IV.3, we identify “pathological” perturbations that destroy
the considered dynamics.
Throughout the next sections, we set the dimension of the Liouville space to
$d=10000$ and the frequency cut-off to $N_{\text{f}}=3333\approx d/3$. We
repeat the strategy $N=1000$ times by drawing random numbers $x_{k},y_{k}$ and
present statistics for the quantifiers $\epsilon$ and $\sigma$.

###  IV.1 Exponential vs. Gaussian decay

In the first round of our stability investigation, an exponential decay
competes against a Gaussian decay. Following the strategy laid out in the
previous sections, the first step is to find suitable Lanczos coefficients
that comply with the four conditions presented in Sec. III.2.
We begin with the Gaussian decay. The Lanczos coefficients that exactly
correspond to Gaussian decay of the form $C(t)=e^{-t^{2}/2}$ are analytically
known [method one in Sec. III.1)], i.e., $b_{n}=\sqrt{n}$ Joslin and Gray
(1986). Nevertheless, the second condition, the compliance with the operator
growth hypothesis, needs to be satisfied. To this end, we choose a cutoff-
point $n^{\star}$ at which the coefficients are continued in a linear fashion.
This yields the following Lanczos coefficients [“g” = “Gaussian”]

$b_{n}^{\text{g}}=\begin{cases}\sqrt{n}\,,\quad\quad\,\,\;n\leq n^{\star}\\\
\alpha n+\gamma\,,\quad n>n^{\star}\,.\end{cases}$ (13)

The parameters $\alpha$ and $\gamma$ are determined by the requirement of a
smooth transition from square-root growth to linear growth, i.e.,
$\alpha=1/(2\sqrt{n^{\star}})$ and $\gamma=\sqrt{n^{\star}}/2$. The change
from square-root to linear growth for $n>n^{\star}$ does not strongly affect
the “Gaussianity” of the correlation dynamics [compare Fig. 2] such that
condition one remains fulfilled.
Next, coefficients for an (approximate) exponential decay need to be found. As
mentioned, a correlation function can never be truly exponential, since it
necessarily features zero slope at $t=0$. Thus, we only require the
exponential decay to be present after a short Zeno-time, which is usually
exceedingly short compared to the relaxation time Breuer and Petruccione
(2006); Joos  _et al._ (2003). We achieve an approximate exponential decay via
an educated guess [method two in Sec. III.1]. Slow dynamics are characterized
by a relatively small first coefficient $b_{1}$, followed by a jump to a
larger remainder of coefficients $b_{n\geq 2}$. Thus, we make the following
ansatz for the coefficients [“e” = “exponential”]

$b_{n}^{\text{e}}=\begin{cases}a\,,\quad\quad\quad\;\;n=1\\\ \alpha
n+\gamma\,,\quad n\geq 2\,.\end{cases}$ (14)

The parameter $a$ is set to $1.2$, a similar magnitude as the first Gaussian
coefficient $b_{1}^{\text{g}}$. The other parameters $\alpha$ and $\gamma$ are
the same as in the Gaussian case. Thus, condition two, i.e., compliance with
the universal operator growth hypothesis, is fulfilled. The Lanczos
coefficients as defined in Eq. (13) and Eq. (14) are depicted in Fig. 1. Since
both sets of coefficients coincide for $n>n^{\star}$, condition four is
satisfied as the quantity $q(b^{\text{g}}_{n},b^{\text{e}}_{n})=0.999993$ is
sufficiently close to unity.

Figure 1: Lanczos coefficients corresponding to Gaussian decay (red) and
exponential decay (blue). The square-root growth transitions to linear growth
at $n=n^{\star}$.

The corresponding dynamics are depicted in Fig. 2. The ansatz for the
coefficients $b^{\text{e}}_{n}$ in Eq. (14) indeed yields a nice exponential
decay. A fit of the form $Ae^{-\lambda t}$ with $A=1.02$ and $\lambda=0.24$
captures the dynamics quite well. Thus, condition one is sufficiently
fulfilled for the exponential decay as well.
Further, we can extract from Fig. 2 that both dynamics decay on more or less
the same time scale. Therefore, we also view condition three as satisfied. If
any, the Gaussian dynamics is faster and thus less prone to perturbations.
Now that we have checked all four conditions presented in Sec. III.2, we are
ready to apply the perturbation as laid out in Sec. III.3. We set
$\lambda=0.5$. Three exemplary perturbed dynamics with respective fits for the
exponential case are depicted in Fig. 3.
It is evident that the perturbed dynamics are still nicely described by an
exponential, only the decay constant changes. For the Gaussian decay, three
exemplary perturbed dynamics with respective fits are depicted in Fig. 4. Two
displayed perturbed Gaussian curves feature oscillations, which can not
possibly be captured by a Gaussian fit ansatz. The three depicted fits are
much worse than in the exponential case. Insets show corresponding Lanczos
coefficients.

Figure 2: Unperturbed correlation functions $C(t)$ calculated from the
respective Lanczos coefficients. Dashed curves indicate an exponential fit
(cyan) and an exact Gaussian $e^{-t^{2}/2}$ (black). Figure 3: Three exemplary
correlation functions originating from the perturbed exponential Lanczos
coefficients $\tilde{b}_{n}^{\text{e}}$. Dashed, black lines indicate
exponential fits. The quantifier $\epsilon$ of all three dynamics is
relatively close to the mean value, see Fig. 5. Inset: corresponding perturbed
Lanczos coefficients.

Figure 4: Three exemplary correlation functions originating from the perturbed
Gaussian Lanczos coefficients $\tilde{b}_{n}^{\text{g}}$. Dashed, black lines
indicate Gaussian fits. The quantifier $\epsilon$ of all three dynamics is
relatively close to the mean value, see Fig. 5. Inset: corresponding perturbed
Lanczos coefficients.

The impression from these six exemplary curves is confirmed in Fig. 5, which
shows histograms of deviations $\epsilon$ of the fits from the respective
perturbed dynamics. The symbol $\Omega$ denotes the number of values
$\epsilon$ within a bin of size $5\cdot 10^{-4}$. There is a clear division
between the exponential cases (blue) and the Gaussian cases (red).
The deviations $\epsilon$ of the exponential decays are much smaller than
those of the Gaussian decays. In particular, the mean value
$\bar{\epsilon}_{\text{e}}=0.002$ (indicated by a blue, dashed line) in the
exponential case is about twenty times smaller than in the Gaussian case
$\bar{\epsilon}_{\text{g}}=0.042$. For a visualization of this comparison see
Fig. 3 and Fig. 4, whose exemplary curves feature values of $\epsilon$ close
to the respective mean values.
We have to be aware that this apparent stability of the exponential may be due
to the possibility that exponential decays are generally less affected by our
constructed perturbation than the Gaussian decays. To check this possibility
and to show that this is not the case, we consider the quantifier $\sigma$,
which measures the difference between the perturbed and the unperturbed
dynamics. The inset of Fig. 5 depicts a scatter plot of all $1000$ pairs
$(\sigma_{i},\epsilon_{i})$. If one decay would be consistently less affected
than the other, a cluster of points to the left [small values of $\sigma$]
would emerge. This is evidently not the case, as the distribution along the
horizontal axis is relatively similar for both decays. The existence of the
diagonal edge is expected, since there can be no fit that is “farther away”
from the perturbed dynamics than the unperturbed dynamics itself, which is, of
course, also an eligible candidate for fitting.
Thus, we conclude that exponential decays is indeed stable with respect to
perturbations [as in Eq. (10)]. In contrast, Gaussian decays seem to be quite
unstable. This is the first main result of the paper at hand.

Figure 5: Histogram of the fit quality measure $\epsilon$ with a bin size of
$5\cdot 10^{-4}$. Dashed lines indicate respective mean values. The stability
of the exponential decay is evident. In contrast, the Gaussian decay does not
seem to be stable, as the deviations $\epsilon$ become quite large. Inset:
scatter plot of all points $(\sigma_{i},\epsilon_{i})$. Both dynamics are
equally affected by the perturbation.

###  IV.2 Damped oscillations: exponential vs. Gaussian

In the second round, we investigate and compare two types of damped
oscillations, one with an exponential damping factor, the other with a
Gaussian damping factor. As in the last section, first, two sets of Lanczos
coefficients need to be found that satisfy all four conditions imposed in Sec.
III.2.
We again begin with the Gaussian case. The coefficients $b_{n}^{\text{gdo}}$
[“gdo” $=$ “Gaussian damped oscillation”] that correspond to an oscillation
damped by a Gaussian factor a neither analytically known, nor can we make an
educated guess. This only leaves the third method, in which we “reverse-
engineer” the coefficients from the dynamics itself. To this end, we choose a
particular correlation function $C(t)=e^{-t^{2}/8}\cos(2t)$ and proceed as
detailed in App. A. This procedure allows us to obtain about $50$ Lanczos
coefficients. After that, the coefficients are continued “by hand” in a manner
that first, respects the pattern exhibited by the coefficients so far, and
second, eventually becomes linear.
For the exponentially damped oscillation we can again make an educated guess
[method two] for the coefficients $b_{n}^{\text{edo}}$ [“edo” $=$
“exponentially damped oscillation”].
We set the first two coefficients $b_{1},b_{2}$ to some values ($b_{1}=2.0$,
$b_{2}=1.6$), followed by a jump to larger coefficients $b_{n\geq 3}$, which
then grow in a linear fashion. The slope is determined by and coincides with
the slope in the Gaussian case.

Figure 6: Lanczos coefficients corresponding to an oscillation damped by a
Gaussian (red) and an oscillation damped by an exponential (blue). Figure 7:
Unperturbed correlation functions $C(t)$ calculated from the respective
Lanczos coefficients. Dashed curves indicate a fit of an exponentially damped
oscillation (cyan) and the oscillation damped by a Gaussian
$e^{-t^{2}/8}\cos(2t)$ (black).

The unperturbed Lanczos coefficients and corresponding correlation functions
are depicted in Fig. 6 and Fig. 7, respectively. The “reverse-engineered and
continued by hand” coefficients $b_{n}^{\text{gdo}}$ indeed reproduce the
given correlation function $C(t)=e^{-t^{2}/8}\cos(2t)$. Further, the guess for
$b_{n}^{\text{edo}}$ yields a nice exponentially damped curve, i.e., the fit
ansatz with [$A=1.04$, $\mu=0.57$, $\omega=0.32$, $\phi=2.19$] nicely captures
the dynamics. Therefore, condition one is satisfied for both cases. The two
sets of coefficients were designed to eventually attain linear growth, thus
condition two is fulfilled. It is evident from Fig. 7 that both dynamics decay
on similar time scales, rendering condition three satisfied. Lastly, the
quantity $q(b_{n}^{\text{gdo}},b_{n}^{\text{edo}})$ practically equals unity.
Thus, all four conditions are met.

Figure 8: Three exemplary correlation functions originating from the perturbed
Lanczos coefficients $\tilde{b}_{n}^{\text{edo}}$. Dashed, black lines
indicate $4$-parametric fits of exponentially damped oscillations. The
quantifier $\epsilon$ of all three dynamics is relatively close to the mean
value, see Fig. 10. The perturbation not only changes the decay constant, but
also the frequency. Inset: corresponding perturbed Lanczos coefficients.
Figure 9: Three exemplary correlation functions originating from the perturbed
Lanczos coefficients $\tilde{b}_{n}^{\text{gdo}}$. Dashed, black lines
indicate $4$-parametric fits of oscillations damped by a Gaussian. The
quantifier $\epsilon$ of all three dynamics is relatively close to the mean
value, see Fig. 10. Inset: corresponding perturbed Lanczos coefficients.

Next, we switch on the perturbation and set $\lambda=0.1$. Exemplary perturbed
dynamics are displayed in Fig. 8 and Fig. 9, respectively. The perturbed
dynamics are fitted with the $4$-parametric ansatz $C(t)=Ae^{-\mu
t}\cos(\omega t-\phi)$ in the exponential and $C(t)=Ae^{-\mu t^{2}}\cos(\omega
t-\phi)$ in the Gaussian case. While the perturbed coefficients in the
exponential case still lead to correlation functions that are within the class
of exponentially damped oscillations [with different decay constants $\mu$ and
frequencies $\omega$], the same can not be said for the Gaussian case, where
the perturbed dynamics decay too slowly and can not by captured by the fit
ansatz above.

Figure 10: Histogram of the fit quality measure $\epsilon$ with a bin size of
$5\cdot 10^{-4}$. Dashed lines indicate respective mean values. The stability
of the exponentially damped oscillations is evident. In contrast, the Gaussian
counterpart does not seem to be stable, as the deviations $\epsilon$ become
quite large.
Inset: scatter plot of all points $(\sigma_{i},\epsilon_{i})$. Both dynamics
are more or less equally affected by the perturbation.

The impression from the six exemplary curves is confirmed in Fig. 10, which
displays histograms of the quantifier $\epsilon$. Again, the division into two
distinct peaks is evident. The respective mean values
$\overline{\epsilon}_{\text{gdo}}=0.026$ and
$\overline{\epsilon}_{\text{edo}}=0.006$ are indicated as dashed, vertical
lines and differ by a factor of about five. For a visualization of what these
values mean for the fit quality, see Fig. 8 and Fig. 9, whose curves feature
values of $\epsilon$ close to the respective mean values. Thus, we conclude
that the exponentially damped dynamics are quite stable and, in particular,
more stable than the oscillations damped by a Gaussian. The disparity between
the two is not as strong as for the decays in the previous section [however,
different values of $\lambda$ may not be directly comparable.]
The inset of Fig. 10 again shows a scatter plot of all points
$(\sigma_{i},\epsilon_{i})$. The blue cluster of dots is a little more
concentrated to small values of $\sigma$, indicating that the exponentially
damped oscillations are a little less altered by the perturbation than their
Gaussian counterpart. However, the marginal distribution over $\sigma$ is
still less partitioned into two peaks than the one over $\epsilon$ [which is
the data in the histogram itself]. Hence, the stability of exponentially
damped oscillations is still due to the nature of the dynamics and not due to
a smaller effect of the perturbation on the dynamics.
We conclude that exponentially damped oscillations are stable with respect to
perturbations [as in Eq. (10)]. In contrast, oscillations damped by a Gaussian
seem to be quite unstable. This is the second main result of the paper at
hand.

Figure 11: Three exemplary dynamics originating from perturbed Lanczos
coefficients $\tilde{b}^{\text{e}}_{n}$. The perturbation consists of
uncorrelated random numbers, i.e., $N_{\text{f}}=d$ in Eq. (10). The
exponential decay is completely destroyed and the correlation function does
not seem to equilibrate at all. Inset: corresponding perturbed Lanczos
coefficients.

###  IV.3 Pathological perturbations

In this section, we lift the restriction imposed on the $v_{n}$ earlier, which
was that frequencies were cut at $N_{\text{f}}\approx d/3$. Instead, we
include all frequencies in the construction of the perturbation, which amounts
to setting $N_{\text{f}}=d$. This case simply corresponds to adding
uncorrelated random numbers to the unperturbed coefficients $b_{n}$.
Otherwise, the strategy is pursued in the same manner as above.
Three exemplary dynamics for exponential decays are depicted in Fig. 11 [for
conciseness, we refrain from showing more exemplary data for the other cases].
As is evident, the exponential decay is completely absent and replaced by
quite irregular dynamics, which do not seem to reach an equilibrium
[$N_{\text{eq}}$ is set as large as possible not  in these cases]. These
curves hint at the presence of localization effects, which prevent the
“particle” [in the tight-binding picture] to leave the first site.
A histogram of the Gaussian vs. exponential data can be viewed in Fig. 12. The
mean values read $\bar{\epsilon}_{\text{g}}=0.15$ in the Gaussian case and
$\bar{\epsilon}_{\text{e}}=0.12$ in the exponential case, which is some orders
of magnitude larger than before when shorter wavelengths were excluded.
Further, the accumulation of points along the diagonal in the inset suggests
that the curves no longer resemble their original form, which is expected
judging from Fig. 11.
The histogram for the oscillating cases is displayed in Fig. 13. The mean
deviations read $\bar{\epsilon}_{\text{gdo}}=0.035$ and
$\bar{\epsilon}_{\text{edo}}=0.025$. These values are comparable to the
Gaussian case in Fig. 10. However, since the quantifier $\epsilon$ measures
more or less the deviation “per time step”, the comparison of values
$\bar{\epsilon}$ between dynamics that do and do not equilibrate can be void
of meaning. In practice, judging from all figures displaying exemplary curves,
a value $\epsilon$ of about one percent corresponds to fits that “look good to
the eye”.

Based on these observations, we specify a property of perturbations that lead
to non-generic or pathological dynamics. Namely, a perturbation is
“untypical”, if it gives rise to Lanczos coefficients that vary non-smoothly,
i.e., there is no minimal correlation length within the coefficients. Again,
what this means on the level of Hermitian matrices is difficult to assess and
beyond the scope of this work. This is the third main result of the paper at
hand.

Figure 12: Histogram of the fit quality measure $\epsilon$ with a bin size of
$5\cdot 10^{-3}$. Dashed lines indicate respective mean values. For both
oscillating cases the deviations $\epsilon$ are quite large. Thus, both
relaxation dynamics are unstable with respect to the perturbation including
all wave lengths. Inset: scatter plot of all points
$(\sigma_{i},\epsilon_{i})$. Points accumulate along the diagonal edge
indicating that the respective dynamics are heavily altered due to the
perturbation.

##  V Conclusion

Figure 13: Histogram of the fit quality measure $\epsilon$ with a bin size of
$5\cdot 10^{-4}$. Dashed lines indicate respective mean values. For both the
exponential and the Gaussian case the deviations $\epsilon$ are extremely
large. Thus, both relaxation dynamics are unstable with respect to the
perturbation including all wave lengths. Inset: scatter plot of all points
$(\sigma_{i},\epsilon_{i})$. Points accumulate along the diagonal edge
indicating that the respective dynamics are heavily altered due to the
perturbation.

In this paper, we performed numerical experiments to probe the stability of
two classes of relaxation dynamics. The first class consisted of exponentially
damped oscillations, which also includes exponential decays. The second class
was chosen as a Gaussian counterpart to the first class, i.e., including
Gaussian decays and oscillations damped by a Gaussian.
The whole strategy was formulated in the framework of the recursion method, in
particular, the perturbations were constructed as an alteration of the Lanczos
coefficients. Unperturbed coefficients $b_{n}$ and perturbation $v_{n}$ were
chosen to satisfy certain physically motivated conditions.
The first main message of the paper at hand is that the exponential class of
dynamics is relatively stable under the considered perturbations. In contrast,
the Gaussian counterpart is found to be quite unstable. These findings confirm
and extent upon previous results based on random matrices Knipschild and
Gemmer (2018), which did not comply with the operator growth hypothesis.
We want to emphasize that within this work the main focus should be put on the
stability of the former class of relaxation dynamics, as these are ubiquitous
in nature, examples are given at the beginning of Sec. IV.
The investigation of the latter class should just be taken as an exemplary
comparison. In fact, the choice of a Gaussian counterpart to the first class
is quite arbitrary. Any number of relaxation dynamics could have been
investigated instead.
The second main message of the paper at hand concerns the nature of the
perturbations themselves. Not only, but also in the context of the works on
“typicality of perturbations” Dabelow and Reimann (2020); Dabelow  _et al._
(2020); Dabelow and Reimann (2021), it could be interesting to find properties
of perturbations that lead to non-generic, “pathological” dynamics. Here, we
identified such a criterion, namely that the perturbation should yield Lanczos
coefficients whose minimal correlation length is still above some threshold
value.
As is evident from Fig. 11, including short correlations seems to lead to
unorthodox dynamics. The displayed curves do not seem to reach an equilibrium,
at least on the available time scale. Rather, localization-like effects are
introduced, which cause some part of the wave function to remain on the first
site. On the other hand, sufficiently smooth coefficients alter the relaxation
dynamics in a “controlled yet non-trivial” manner.
Perturbations in various numerical investigations based in random matrices
Knipschild and Gemmer (2018); Dabelow and Reimann (2020) as well as spin
lattice models Heveling  _et al._ (2020); Richter  _et al._ (2020) seem to
naturally possess the above specified property. Hence, a more systematic
investigation on the existence of a minimal correlation length in the
coefficients for realistic setups could be a possible prospect for future
research.

## Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft (DFG) within
the Research Unit FOR 2692 under Grant No. 397107022 (GE 1657/3-2) and No.
397067869 (STE 2243/3-2).

## References

  * Gogolin and Eisert (2016) C. Gogolin and J. Eisert, “Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantum systems,” Reports on Progress in Physics 79, 056001 (2016).
  * Srednicki (1994) M. Srednicki, “Chaos and quantum thermalization,” Phys. Rev. E 50, 888–901 (1994).
  * Deutsch (1991) J.M. Deutsch, “Quantum statistical mechanics in a closed system,” Phys. Rev. A 43, 2046–2049 (1991).
  * D’Alessio  _et al._ (2016) L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, “From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics,” Advances in Physics 65, 239–362 (2016).
  * Lloyd (2013) S. Lloyd, “Pure state quantum statistical mechanics and black holes,” (2013), arXiv:1307.0378 .
  * Goldstein  _et al._ (2006) S. Goldstein, J.L. Lebowitz, R. Tumulka, and N. Zanghì, “Canonical typicality,” Phys. Rev. Lett. 96, 050403 (2006).
  * Reimann (2007) P. Reimann, “Typicality for generalized microcanonical ensembles,” Phys. Rev. Lett. 99, 160404 (2007).
  * Knipschild and Gemmer (2018) L. Knipschild and J. Gemmer, “Stability of quantum dynamics under constant Hamiltonian perturbations,” Phys. Rev. E 98, 062103 (2018).
  * Dabelow and Reimann (2020) L. Dabelow and P. Reimann, “Relaxation theory for perturbed many-body quantum systems versus numerics and experiment,” Phys. Rev. Lett. 124, 120602 (2020).
  * Dabelow  _et al._ (2020) L. Dabelow, P. Vorndamme, and P. Reimann, “Modification of quantum many-body relaxation by perturbations exhibiting a banded matrix structure,” Phys. Rev. Research 2, 033210 (2020).
  * Dabelow and Reimann (2021) L. Dabelow and P. Reimann, “Typical relaxation of perturbed quantum many-body systems,” Journal of Statistical Mechanics: Theory and Experiment 2021, 013106 (2021).
  * Viswanath and Müller (2008) V.S. Viswanath and G. Müller,  _The Recursion Method: Applications to Many-Body Dynamics_ (Springer, New York, 2008).
  * Mattis (1981) D.C. Mattis, “How to reduce practically any problem to one dimension,” in  _Physics in One Dimension_ (Springer Berlin Heidelberg, 1981) pp. 3–10.
  * Parker  _et al._ (2019) D.E. Parker, X. Cao, A. Avdoshkin, T. Scaffidi, and E. Altman, “A Universal Operator Growth Hypothesis,” Phys. Rev. X 9, 041017 (2019).
  * Noh (2021) J.D. Noh, “Operator growth in the transverse-field ising spin chain with integrability-breaking longitudinal field,” Phys. Rev. E 104, 034112 (2021).
  * Heveling  _et al._ (2022) R. Heveling, J. Wang, and J. Gemmer, “Numerically probing the universal operator growth hypothesis,” (2022), arXiv:2203.00533 .
  * Anderson (1958) P.W. Anderson, “Absence of diffusion in certain random lattices,” Phys. Rev. 109, 1492–1505 (1958).
  * (18) In practice, we define the equilibration time as follows. We choose a threshold value of $0.005$. If the absolute value of $C(t)$ is smaller than the threshold value over the course a certain time interval [we choose three time units], the end-time corresponds to the upper bound in Eq. (11) and Eq. (12). The time step is $\delta t=0.01$. If the correlation function does not equilibrate at all, $N_{\text{eq}}$ is set to the maximal time before an echo due to the finiteness of the system occurs .
  * Joslin and Gray (1986) C.G. Joslin and C.G. Gray, “Calculation of transport coefficients using a modified Mori formalism,” Molecular Physics 58, 789–797 (1986).
  * Breuer and Petruccione (2006) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2006).
  * Joos  _et al._ (2003) E. Joos, H.D. Zee, C. Kiefer, D.J.W. Giulini, J. Kupsch, and I.-O. Stamatescu,  _Decoherence and the Appeareance of a Classical World from Quantum Theory_ (Springer, 2003).
  * Heveling  _et al._ (2020) R. Heveling, L. Knipschild, and J. Gemmer, “Modeling the impact of Hamiltonian perturbations on expectation value dynamics,” Zeitschrift für Naturforschung A 75, 475–481 (2020).
  * Richter  _et al._ (2020) J. Richter, F. Jin, L. Knipschild, H. De Raedt, K. Michielsen, J. Gemmer, and R. Steinigeweg, “Exponential damping induced by random and realistic perturbations,” Phys. Rev. E 101, 062133 (2020).

##  Appendix A Reverse-engineering $b_{n}$ from $C(t)$

In this section, we present the procedure to reverse-engineer the Lanczos
coefficients $b_{n}$ from a given correlation function $C(t)$. The idea is to
employ the Lanczos algorithm from Sec. II, but on the level of spectral
functions $\Phi(\omega)$ rather than on the level of observables
$|\mathcal{O})$. The inner product of operators turns into an inner product of
functions, i.e.,

$(\Phi_{1}|\Phi_{2})=\int\Phi_{1}^{\star}(\omega)\Phi_{2}(\omega)\,\text{d}\omega\,.$
(15)

Further, the application of the Liouvillian $\mathcal{L}$ to an operator
$\mathcal{O}$ corresponds to a multiplication of the respective function in
Fourier space with $-\omega$ [since the commutator with $\mathcal{H}$
corresponds to a time derivative, which is equivalent to a multiplication with
$\mathrm{i}\omega$ in Fourier space],

$\mathcal{L}|\mathcal{O})\rightarrow-\omega\Phi(\omega)\,.$ (16)

In practice, we choose a specific correlation function $C(t)$ and calculate
its Fourier transform $\Phi(\omega)$. The initial “seed” function is set to
the (normalized) $\sqrt{\Phi(\omega)}$. Then, the Lanczos algorithm operates
as laid out in Sec. II. In this manner, about $50$ coefficients can be
obtained before numerical instabilities become too pronounced.

##  Appendix B Link to spectral width

The eigenvalue equation for the Hamiltonian reads
$\mathcal{H}|E_{i}\rangle=E_{i}|E_{i}\rangle$, where $E_{i}$ and
$|E_{i}\rangle$ denote eigenvalues and eigenstates, respectively. The
corresponding eigenvalue equation for the Liouvillian superoperator is given
by $\mathcal{L}\mathcal{M}_{\beta}=\mathcal{E}_{\beta}\mathcal{M}_{\beta}$
with eigenvalues $\mathcal{E}_{\beta(i,j)}=E_{i}-E_{j}$ and eigenoperators
$\mathcal{M}_{\beta(i,j)}=|E_{i}\rangle\langle E_{j}|$. Without loss of
generality, we can set $\text{Tr}[\mathcal{H}]=\sum_{i}E_{i}=0$. Further, we
denote the dimension of the Hilbert space by $d_{\mathcal{H}}$, i.e.,
$d_{\mathcal{H}}^{\hskip 1.00374pt2}=d$. Then we have

$\displaystyle\sum_{n}b_{n}^{2}$
$\displaystyle=\dfrac{1}{2}\text{Tr}[\mathcal{L}^{2}]=\dfrac{1}{2}\sum_{i,j}(E_{i}-E_{j})^{2}$
(17) $\displaystyle=\dfrac{1}{2}\sum_{i,j}(E_{i}^{2}-2E_{i}E_{j}+E_{j}^{2})$
$\displaystyle=\dfrac{d_{\mathcal{H}}}{2}\big{(}\sum_{i}E_{i}^{2}+\sum_{j}E_{j}^{2}\big{)}$
$\displaystyle=d_{\mathcal{H}}\text{Tr}[\mathcal{H}^{2}]$

Thus, we see that the quantity in question is indeed linked to the spectral
width of the Hamiltonian.

##  Appendix C Scaling with perturbation strength

Consider a Hamiltonian $\mathcal{H}=\mathcal{H}_{0}+\lambda\mathcal{V}$
consisting of an unperturbed part $\mathcal{H}_{0}$ and a perturbation
$\lambda\mathcal{V}$ with $\text{Tr}[\mathcal{H}_{0}\mathcal{V}]=0$. For the
spectral variance of $\mathcal{H}$ it holds true that

$\text{Tr}[\mathcal{H}^{2}]=\text{Tr}[\mathcal{H}_{0}^{2}]+\lambda^{2}\text{Tr}[\mathcal{V}^{2}]\,.$
(18)

On the other hand, with the particular choice of perturbation in Eq. (9), we
get that

$\displaystyle\sum_{n}\tilde{b}_{n}^{2}$ $\displaystyle=\sum_{n}(b_{n}+\lambda
v_{n})^{2}$ (19) $\displaystyle=\sum_{n}(b_{n}^{2}+2\lambda
b_{n}v_{n}+\lambda^{2}v_{n}^{2})$
$\displaystyle\approx\sum_{n}b_{n}^{2}+\lambda^{2}\sum_{n}v_{n}^{2}\,.$

The term $2\lambda b_{n}v_{n}$ is negligible due to the specific choice of
$v_{n}$. Recalling the relation in Eq. (17), we see from Eq. (19) that the
scaling in Eq. (9) reproduces the scaling in Eq. (18) and is therefore
appropriate.


#  Properties of moments of density for nonlocal mean field game equations
with a quadratic cost function

Olga S. Rozanova*, Mikhail V. Inyakin  Mathematics and Mechanics Department,
Lomonosov Moscow State University, Leninskie Gory, Moscow, 119991, Russian
Federation<EMAIL_ADDRESS>
###### Abstract.

We consider mean field game equations with an underlying jump-diffusion
process $X_{t}$ for the case of a quadratic cost function and show that the
expectation and variance of $X_{t}$ obey second-order ordinary differential
equations with coefficients depending on the parameters of the cost function.
Moreover, for the case of pure diffusion, the characteristic function and the
fundamental solution of the probability density can be expressed in terms of
the expectation ${\mathbb{E}}$ and the variance ${\mathbb{V}}$ of the process
$X_{t}$, so that the moments of any order depend only on ${\mathbb{E}}$ and
${\mathbb{V}}$.

###### Key words and phrases:

mean field game, quadratic Hamiltonian, quadratic cost, moments of density

###### 1991 Mathematics Subject Classification:

Primary 35Q89; Secondary 60H30

Mean field game theory (MFG) has been intensively developed in the last decade
in terms of various mathematical disciplines such as partial differential
equations, control theory, probability theory and numerical methods. It also
has numerous practical applications in all areas where the behavior of a large
number of agents is critical, e.g. in economics, mobile network design,
logistics, etc. [11], [10], [6], [3].

In this paper we consider the forward-backward MFG problem for a particular
case of a coupled PDE system, generally speaking, nonlocal, and show that,
despite the impossibility of obtaining an explicit form of the density, the
behavior of its moments is relatively simple and often amenable to analytical
study.

1\. Let $X$ be a stochastic process with dynamics given by

$dX_{s}=\alpha_{s}\,ds+\delta\,dW_{s}+\lambda d\Gamma_{s},\quad X_{0}=x_{0},$
(1)

$x_{0}\in\mathbb{R}^{n}$ is a point in the space of states, $0\leq s\leq T$,
$W_{s}$ is a standard vectorial Brownian motion, $\Gamma_{s}$ is a purely
discontinuous vector process with jump size given by an independent random
variable, $\delta\geq 0$, $\lambda\geq 0$ are constants. Let $m(t,x)\geq 0$ be
a smooth probability density of underlying stochastic process $X_{s}$,
$\int\limits_{{\mathbb{R}}^{n}}m(t,x)=1$. We denote $p(z)$ a probability
density of jumps, $\int\limits_{{\mathbb{R}}^{n}}\,p(z)dz=1$, and assume
$\mathcal{M}_{1}=\int\limits_{{\mathbb{R}}^{n}}\,zp(z)dz<\infty$,
$\mathcal{M}_{2}=\int\limits_{{\mathbb{R}}^{n}}\,|z|^{2}p(z)dz<\infty$.

The mean field game optimization problem is to define maximum over all the
admissible Markov controls $\alpha_{s}\in
L_{\mathcal{F}}^{2}(0,T;{\mathbb{R}}^{n}$ of the functional

$J(\alpha_{s})={\mathbb{E}}\Bigg{[}\int_{0}^{T}\left(-\frac{|\alpha_{s}|^{2}}{2}+g(s,X_{s},m)\right)\,ds+K(X_{T})\Bigg{]},$

$g:{\mathbb{R}}_{+}\times{{\mathbb{R}}^{n}}\times{\mathbb{R}}_{+}\to\mathbb{R}$
and $K:{\mathbb{R}}^{n}\to{\mathbb{R}}$ are prescribed continuous functions,
the process $X_{s}$ obeys (1),
$\mathcal{F}_{t}=\sigma(x_{0},W_{s},\Gamma_{s},s\leq t)$.

Analytically this problem reduces to the solution of the following system
initial-terminal problem for coupled Hamilton-Jacobi-Bellman (HJB) and
Kolmogorov-Feller-Fokker-Planck (KFFP) equations for the density $m(t,x)$ and
the pay-off function $\Phi(t,x)$ [11], [12], [13]:

$\displaystyle\partial_{t}\Phi+\frac{1}{2}(|\nabla\Phi|)^{2}+\frac{\delta^{2}}{2}\Delta\Phi+\lambda\left(\int\limits_{\mathbb{R}}\Phi(t,x+z)p(z)dz-\Phi(x)\right)=-g(t,x,m),$
(2) $\displaystyle\partial_{t}m+{\rm
div}(m\nabla\Phi)-\frac{\delta^{2}}{2}\Delta
m-\lambda\left(\int\limits_{\mathbb{R}}m(t,x-z)p(z)dz-m(x)\right)=0,$ (3)
$\displaystyle m(0,x)=m_{0}(x),\qquad\Phi(x,T)=K(x).$ (4)

The main interest is the study of density evolution in the process of control.

We are going to show that, with a special choice of control, the first two
density moments can be found as solutions to second-order ODEs, and the
characteristic function is expressed in terms of them (and, perhaps, in terms
of other parameters of the problem). Therefore, for any initial density, its
full behavior can be restored.

The main assumption that allows us to prove our main result is as follows:

$g=a(t)|x|^{2}+b(t)\cdot x+c(t),$ (5)

where $a,$ $b=(b_{1},\dots,b_{n}),$ $c$ are smooth functions ($b$ is a vector)
that can depend on the expectation ${\mathbb{E}}(X_{s})$ and variance
${\mathbb{V}}(X_{s})$, and in this sense the control depends on the density,
and

$K(x)={A}_{T}|x|^{2}+{B}_{T}\cdot x+{C}_{T},$ (6)

with constant $A_{T}$, ${B}_{T}$, $C_{T}$.

Similar assumptions about the cost function were made, for example, in [9],
[2], [17], [5] for the case of pure diffusion, where it was shown that the
solution can be found by solving the matrix Riccati equation. Our results are
also based on the fact that the solution can be partially expressed in terms
of the solution of the Riccati equation.

2\. First we find a solution to the HJB equation.

######  Lemma 1.

Assume that $g$ has the form (5). Then (2) has a solution of the form

$\Phi={A}(t)|x|^{2}+{B}(t)\cdot x+{C}(t).\quad$ (7)

The coefficients $A(t)$, $B(t)$, $C(t)$ are uniquely defined by the terminal
conditions (6).

###### Proof.

We substitute (7) to (2) and get a polynomial of the second order with respect
to $x$. Then we combine the coefficients at $x^{2}$, $x$ and $1$ and obtain a
system of ODEs for $A,B,C$, subject to the terminal conditions,

$\displaystyle\dot{A}=-2A^{2}-{a},\quad\dot{B}+2AB=-2{\lambda}\mathcal{M}_{1}\,A-{b},$
(8)
$\displaystyle\dot{C}=-\frac{1}{2}B^{2}-\delta^{2}A-\lambda(\mathcal{M}_{2}A+\mathcal{M}_{1}B),$
$\displaystyle{A(T)=A_{T}},\quad B(T)=B_{T},\quad C(T)=C_{T}.$

The equation for $A$ splits from the rest of system and defines the entire
dynamics. The equation for $B$ is linear and contain $A$ in coefficients.
$\Box$

It what follows we assume that $a(t)$ and $b(t)$ are such that

$e^{2\int\limits_{0}^{T}A(\tau)d\tau}<\infty,\qquad\int_{0}^{T}e^{2\int\limits_{\eta}^{T}A(\tau)d\tau}B(\eta)\,d\eta\,<\infty.$
(9)

3\. Now we find the fundamental solution of the KFFP equation. If $\Phi(t,x)$
is known, (3) takes a more specific form,

$\displaystyle\partial_{t}m+{(2A\,x+B)}\cdot\nabla
m+{2A}m-\frac{\delta^{2}}{2}\Delta
m-\lambda\left(\int_{-\infty}^{\infty}m(t,x-z)p(z)dz-m\right)=0,$ (10)

with initial condition

$m_{0}(x)=\delta(x).$ (11)

We denote the solution to (10), (11) as $\mathcal{G}(t,x)$.

Applying the Fourier transform $x\to\omega$, we obtain the Cauchy problem for
the characteristic function $\hat{m}=\hat{m}(t,\omega)$:

$\displaystyle\partial_{t}\hat{m}-2{A(t)}\omega\partial_{\omega}\hat{m}+\left(\frac{\delta^{2}}{2}\omega^{2}+iB(t)\omega-(\hat{p}(\omega)-1)\lambda\right)\hat{m}=0,\quad\hat{m}(0,\omega)=1,$
(12)

where it is assumed that $A$ and $B$ are known from the previous step. The
solution to (12) can be explicitly found, namely,

$\displaystyle\hat{G}(t,\omega)={\exp}\left[-\int_{0}^{t}\left(\frac{\delta^{2}}{2}{\mathcal{R}}^{2}+iB(\eta){\mathcal{R}}-(\hat{p}({\mathcal{R}})-1)\lambda\right)d\eta\right],$
(13) $\displaystyle\quad{\mathcal{R}}={\mathcal{R}}(t,\eta,\omega)=\omega
e^{2R(t,\eta)},\quad R(t,\eta)=\int\limits_{\eta}^{t}A(\tau)d\tau.$

4\. We find the expectation and variance in terms of $A$ and $B$.

######  Lemma 2.

Let $X_{t}$ be a random process with an initial probability density
$m_{0}(x)$, $x_{i}^{k}\,m_{0}(x)\in L_{1}({\mathbb{R})}$, $k=1,2$,
$i=1,\dots,n$. The expectation and variance of $X_{t}$ can be found as

$\mathbb{E}(t)=\int_{0}^{t}e^{2\int\limits_{\eta}^{t}A(\tau)d\tau}(B(\eta)+\lambda\mathcal{M}_{1})\,d\eta+\mathbb{E}(0),$
(14)
$\mathbb{V}(t)=(\delta^{2}+\lambda\mathcal{M}_{2})\,\int_{0}^{t}e^{4\int\limits_{\eta}^{t}A(\tau)d\tau}\,d\eta+\mathbb{V}(0),$
(15)

where $A(t),B(t)$ is a solution to (8), $t\in[0,T].$

###### Proof.

First assume $m_{0}(x)=\delta(x)$. The solution of (10) in this case is the
fundamental solution $\mathcal{G}(t,x)$. Then, by means of standard
computations we get from (13)

$\mathbb{E}(X^{0}_{t})=i\frac{\partial\hat{\mathcal{G}}}{\partial\omega}(t,0)=\int_{0}^{t}e^{2R(t,\eta)}\left(B(\eta)+{\lambda}{\mathcal{M}}_{1}\right)\,d\eta,$
$\mathbb{V}(X^{0}_{t})=(\mathbb{E}(X^{0}_{t}))^{2}-\frac{\partial^{2}\hat{\mathcal{G}}}{\partial\omega^{2}}(t,0)=(\delta^{2}+\lambda\mathcal{M}_{2})\,\int_{0}^{t}e^{4\int\limits_{\eta}^{t}A(\tau)d\tau}\,d\eta.$

The solution to (10) with initial data $m_{0}(0)$, the density of the process
$X_{t}\|_{t=0}=X_{0}$, can be found as a convolution of $\mathcal{G}$ with
$m_{0}(x)$, therefore
$\hat{m}(t,\omega)=\hat{G}(t,\omega)\,\hat{m}_{0}(\omega)$,

$\mathbb{E}(X_{t})=i\frac{\partial\hat{m}(t,0)}{\partial\omega}=i\frac{\partial\hat{\mathcal{G}}(t,0)}{\partial\omega}\hat{m}_{0}(0)+i\frac{\partial\hat{m}_{0}(0)}{\partial\omega}\hat{\mathcal{G}}(t,0)=\mathbb{E}(X_{t}^{0})+\mathbb{E}(X_{0}),$
$\mathbb{V}(X_{t})=\left(\frac{\partial\hat{m}(t,0)}{\partial\omega}\right)^{2}-\frac{\partial^{2}\hat{m}(t,0)}{\partial\omega^{2}}=\mathbb{V}(X^{0}_{t})+\mathbb{V}(X_{0}),$

so we get (14), (15). Here we take into account that
$\hat{\mathcal{G}}(t,0)=1$ and
$\hat{m}_{0}(0)=\int\limits_{{\mathbb{R}}^{n}}m_{0}(y)dy=1$. $\Box$

5\. Now we express $\hat{\mathcal{G}}$ in terms of $\mathbb{E}$ and
$\mathbb{V}$.

######  Theorem 1.

Assume that condition (9) holds. The characteristic function
$\hat{\mathcal{G}}$ can be expressed in terms of $\mathbb{E}$, $\mathbb{V}$
and $A$ as follows:

$\displaystyle\hat{\mathcal{G}}(t,w)=\exp\\{-\frac{1}{2}|w|^{2}\mathbb{V}(t)-iw\cdot\mathbb{E}(t)+\lambda\mathcal{Q}(t,w)\\},$

where
$Q(t,w)=\int\limits_{0}^{t}(\hat{p}({\mathcal{R}})-1+I\mathcal{M}_{1}\cdot{\mathcal{R}})d\eta,\quad$
${\mathcal{R}}(t,\eta,\omega)=\omega e^{2\int\limits_{\eta}^{t}A(\tau)d\tau}$
($\mathcal{R}$ and $M_{1}$ are vectors).

For $\lambda=0$ the characteristic function depends on $\mathbb{E}$,
$\mathbb{V}$ only, therefore every moment of $X_{t}$ can be expressed in terms
of its expectation and variance.

The proof follows from (13) and Lemma 2 directly. Note that $A$ is defined by
$a$ and the terminal condition, see (8). $\Box$

######  Corollary 1.

In the case $\lambda=0$ the fundamental solution is

$\displaystyle{\mathcal{G}}(t,x)=\frac{1}{(\sqrt{2\pi\mathbb{V}(t)})^{n}}\,e^{-\frac{|x-\mathbb{E}(t)|^{2}}{2\mathbb{V}(t)}}.$

6\. Now we derive equations for the expectation and variance.

######  Theorem 2.

Assume that condition (9) holds. The expectation of $\mathbb{E}(t)$ and the
variance $\mathbb{V}(t)$ satisfy the following boundary value problem

$\displaystyle\mathbb{E}^{{}^{\prime\prime}}(t)+2a(t)\mathbb{E}(t)=-b(t),$
(16)
$\displaystyle\mathbb{E}(0)=x_{0},\quad\mathbb{E}(T)=\int_{0}^{T}e^{2\int\limits_{\eta}^{T}A(\tau)d\tau}\left(B(\eta)+{\lambda}{\mathcal{M}}_{1}\right)\,d\eta+x_{0},$
$\displaystyle\mathbb{V}^{{}^{\prime\prime}}(t)+4a(t)\mathbb{V}(t)-\displaystyle\frac{(\mathbb{V}^{{}^{\prime}}(t))^{2}-K^{2}}{2\mathbb{V}(t)}=0,\quad
K=\delta^{2}+\lambda{\mathcal{M}}_{2},$ (17)
$\displaystyle\mathbb{V}(0)=v_{0},\quad\mathbb{V}(T)=(\delta^{2}+\lambda\mathcal{M}_{2})\,\int_{0}^{T}e^{4\int\limits_{\eta}^{T}A(\tau)d\tau}\,d\eta+v_{0},\quad
t\in[0,T].$

###### Proof.

The statement follows from (14) and (15) by direct computation. The equation
(16) was obtained in [14]. $\Box$

7\. So far, we have not indicated the dependence of $a$ and $b$ on
$\mathbb{E}$ and $\mathbb{V}$. We can use various hypotheses about this
dependence and obtain a system (16), (17), which is non-linear and, generally
speaking, cannot be solved analytically. However, in the simplest cases, the
solution can be found explicitly. Corollary 2 [14] gives one of these
possibilities.

######  Corollary 2.

If $a(t)=\rm const$ and
$b(t)=b_{0}+b_{1}\mathbb{E}(t)+b_{2}\mathbb{E}^{\prime}(t)$,
$b_{0},b_{1},b_{2}=\rm const$, then (16) transforms into a linear second order
ODE

$\mathbb{E}^{{}^{\prime\prime}}+b_{2}\mathbb{E}^{\prime}+(2a+b_{1})\mathbb{E}=-b_{0},$

which can be explicitly solved.

However, the simplest case of an explicit solution is constants $a$ and $b$.
Namely,

i) For $a>0$

$\displaystyle\mathbb{E}(t)=C_{1}\sin\sqrt{2t}+C_{2}\cos\sqrt{2a}t-\frac{b}{2a},$
(18)
$\displaystyle\mathbb{V}(t)=\pm\frac{1}{a}\,\sqrt{a\left((C_{1}^{2}+C_{2}^{2})a+\frac{K^{2}}{8}\right)}+C_{1}\sin(2\sqrt{2at})+C_{2}\cos(2\sqrt{2at}),$
(19)

the sign is chosen for reasons of non-negativity $\mathbb{V}$.

In this case

$\displaystyle A(t)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{a}{2}}\,\tan\theta(t,T),\qquad\theta(t,T)=\arctan\sqrt{\frac{2}{a}}A_{T}+\sqrt{2a}(T-t),$
$\displaystyle B(t)$ $\displaystyle=$
$\displaystyle-\frac{b}{\sqrt{2a}}\,\tan\theta(t,T)+\frac{bA_{T}+aB_{T}}{\sqrt{a(a+2A_{T}^{2})}\cos\theta(t,T)}.$

Note that the solution of (8) may have singularities inside $(0,T)$,
nevertheless condition (9) holds, so the singularities of $A$ do not interfere
with the calculation of the boundary conditions.

It is possible to show (see [15]), that $\mathbb{E}(t)$ is differentiable in
the points $t$, where $A$ and $B$ fail to exist. In turn, $\mathbb{V}(t)$ is
continuous but not differentiable when $\mathbb{V}(t)=0$.

ii) For $a<0$

$\displaystyle\mathbb{E}(t)=C_{1}\sinh\sqrt{2t}+C_{2}\cosh\sqrt{2a}t-\frac{b}{2a},$
(20)
$\displaystyle\mathbb{V}(t)=C_{1}+\frac{8aC_{1}^{2}+K^{2}}{32aC_{2}}\exp\left(2\sqrt{-2a}t\right)+C_{2}\exp\left(-2\sqrt{-2a}t\right).$
(21)

Functions $A$ and $B$, necessary for boundary conditions, can be found from
(8), the problem of singularities does not arise.

iii) For $a=0$ both $\mathbb{E}(t)$ and $\mathbb{V}(t)$ are polynomials of the
second order with respect to $t$.

8\. Choosing $g$ in the form (5) seems artificial. However, it appears
naturally in a number of problems related to behavioral economics, e.g. [4],
[15], [16]. Note that in [4] a second-order equation was obtained for the
expectation under different assumptions. The MPG equations with nonlocal terms
were considered in [7].

It is well known that in the case of pure diffusion and a quadratic cost
function, one can find an explicit solution for the density in the form of a
Gaussian function [9]. In [15] an equation for the position of the maximum of
density was obtained, which coincides for this case with the expectation. For
the case of a jump diffusion, these two characteristics of the density
function do not coincide.

9\. Since $\mathbb{E}(t)$ and $\mathbb{V}(t)$ are observable for real
processes and can be estimated by statistical methods, we can propose a
procedure for restoring the parameters of the cost functional under the
assumption that the penalty term $g$ is quadratic. For example, assuming
constant $a$ and $b$, one can use the formulas (18), (19) or (20), (21)
depending on the oscillatory or non-oscillatory character of $\mathbb{E}(t)$
and $\mathbb{V}(t)$ in time. We can refer to [8], [1] and their references to
prove that both behaviors are possible.

10\. Acknowledgments. Supported by the Moscow Center for Fundamental and
Applied Mathematics under the agreement №075-15-2019-1621.

## References

  * [1] Agliari A., Hommes C.H., Pecora N., Path dependent coordination of expectations in asset pricing experiments: A behavioral explanation, Journal of Economic Behavior and Organization, 121 (2016), 15–28.
  * [2] Bensoussan A., Sung K. C. J. , Yam S. C. P., Yung S. P., Linear-quadratic mean field games, J. Optim. Theory Appl., 169 (2016), 496–529.
  * [3] Caines P.E., Huang M., Malhamé R.P., Mean Field Games. Handbook of Dynamic Game Theory (Eds. Basar T., Zaccour G.) Springer, 2017.
  * [4] Cardaliaguet P., Lehalle C.-A.,  Mean Field Game of Controls and An Application To Trade Crowding. Mathematics and Financial Economics, 12 (2018), no 3, 335-363.
  * [5] Fatone L., Mariani F., Recchioni M., Zirilli F., A trading execution model based on mean field games and optimal control, Appl. Math. 5 (2014), 3091–3116.
  * [6] Gomes D. A., Saede J.,  Mean Field Games Models – A Brief Survey. Dynamic Games and Applications, 4 (2013), 110–154.
  * [7] Graber, P.J., Ignazio V., Neufeld A.,  Nonlocal Bertrand and Cournot mean field games with general nonlinear demand schedule, Journal de Mathématiques Pures et Appliquées, 148 (2021) 150–198.
  * [8] Hommes C.H., Sonnemans J., Tuinstra J., Velden H. van de, Coordination of expectations in asset pricing experiments. Rev. Financ. Stud., 18 (3) (2005), 955–980.
  * [9] Guéant, O., A reference case for mean field games models, J. Math. Pures Appl. 92 (2009), 276–294.
  * [10] Guéant O., Lasry J.M., Lions P.L.,Mean Field Games and applications. Paris-Princeton lectures on mathematical finance, Springer, 2010, 205–266.
  * [11] Lasry, J.-M., Lions, P.-L,  Mean field games. Jpn. J. Math. 2 (2007) 229–260.
  * [12] Øksendal B., Stochastic differential equations. Introduction in theory and applications. 5th Ed., Springer, Heidelberg–New-York, 2000.
  * [13] Øksendal B., Sulem A., Applied Stochastic Control of Jump Diffusions. Springer, Berlin–Heidelberg–New-York, 2005.
  * [14] Rozanova O., Manapov I. Mean field game equations with underlying jump-diffusion process, AIP Conference Proceedings 2522, 060011 (2022); https://doi.org/10.1063/5.0100745.
  * [15] Rozanova O.S., Nikulin S.I., On certain analytically solvable problems of mean field games theory, Moscow University Mathematics Bulletin, 75 (2020) 139–148.
  * [16] Trusov N.,  Application of mean field games approximation to economic processes modeling (in Russian), Tr. ISA RAN 68 (2018), no. 2, 88–91.
  * [17] Yong J., Linear-Quadratic Optimal Control Problems for Mean-Field Stochastic Differential Equations, SIAM J. Control Optim., 51 (2013), 2809–2838.