# Femtosecond atomic-scale X-ray photon correlation spectroscopy enables
direct probing of relaxation dynamics in fragile liquid

Tomoki Fujita1 Yanwen Sun2∗ Haoyuan Li2 Thies J. Albert3 Sanghoon Song2
Takahiro Sato2 Jens Moesgaard1 Antoine Cornet4 Peihao Sun5 Ying Chen2 Mianzhen
Mo2 Narges Amini1 Fan Yang6 Arune Makareviciute1 Garrett Coleman7 Pierre
Lucas7 Jan Peter Embs8 Michael Sprung9 Vincent Esposito2 Joan Vila-Comamala8
Nan Wang2 Talgat Mamyrbayev8 Christian David8 Peter Zalden10 Jerome Hastings2
Beatrice Ruta4 Paul Fuoss2 Klaus Sokolowski-Tinten3 Diling Zhu2 Shuai Wei1,11∗
1Department of Chemistry, Aarhus University, 8000 Aarhus, Denmark
2Linac Coherent Light Source, SLAC National Accelerator Laboratory, Menlo
Park, CA 94025, USA
3Faculty of Physics and Center for Nanointegration Duisburg-Essen (CENIDE),
University of Duisburg-Essen, Lotharstrasse 1, 47048 Duisburg, Germany
4Institut Néel, CNRS, 38042 Grenoble, France
5Department of Physics, Università degli Studi di Padova, 35122 Padova, Italy
6Institute of Materials Physics, Deutsches Zentrum für Luft- und Raumfahrt
(DLR), Köln, Germany
7Department of Materials Science and Engineering, University of Arizona,
Tucson AZ 85721, USA
8Paul Scherrer Institut, 5232 Villigen PSI, Switzerland
9Deutsches Elektronen-Synchrotron DESY, 22607 Hamburg, Germany
10European XFEL, 22869 Schenefeld, Germany
11Centre for Integrated Materials Research, Aarhus University, 8000 Aarhus,
Denmark
∗To whom correspondence should be addressed; E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
> Structural relaxations are intrinsic to all liquids and glasses. They are
> the fundamental processes underlying many material properties. However, the
> mechanism of structural relaxations is still unclear due to the experimental
> difficulties in directly characterizing the atomic-scale relaxations. The
> latter exhibit a wide range of timescales from geological times to sub-
> picoseconds. Here we show that, using split-delay coherent pulses from X-ray
> free electron lasers (XFEL), we can achieve the femtosecond temporal and
> atomic-scale spatial resolutions of the X-ray photon correlation
> spectroscopy (fs-XPCS). The latter allows us to directly observe atomic-
> scale relaxation dynamics and determine the intermediate scattering function
> in a high-fluidity fragile liquid. We find that the beam characteristics of
> the double pulses from the split-delay drift over time, the correction of
> which is essential to obtain the accurate intermediate scattering functions
> from speckle contrasts. Our finding demonstrates a new avenue of direct
> experimental studies of relaxation dynamics with the femtosecond resolution,
> a 12 orders of magnitude extension from the sub-second regime of the
> synchrotron-based XPCS.

##  1 Introduction

Liquids and glasses exhibit structural relaxations with atomic rearrangements
toward their energetically favorable positions. The relaxation dynamics are
extremely diverse with the characteristic timescales ranging from millions of
years for glass aging to some sub-picoseconds for fast relaxation processes in
high-fluidity liquids. They are decisive to many properties of glass forming
systems, such as viscosity, vitrification, amorphous stability, and
crystallization (?, ?, ?, ?, ?). Relaxation dynamics are usually characterized
by measuring time (or frequency)-dependent changes of macroscopic properties
using techniques such as dynamic mechanical spectroscopy, dielectric
spectroscopy, dynamic light scattering, calorimetry and rheology (?, ?, ?, ?,
?). However, the atomic-level mechanism of structural relaxations is often
debated, as few experimental technique can directly probe the atomic-scale
structural relaxations in disordered systems. In the past decades, X-ray
photon correlation spectroscopy (XPCS) based on synchrotron X-ray sources has
been developed to probe the atomic-scale structural relaxations (?, ?, ?). The
intensity autocorrelation functions of measured speckle patterns can be
related to intermediate scattering functions (ISF), which reveal the atomic-
scale dynamics (?). However, XPCS has been limited to slow dynamics (i.e.
typically from milliseconds to hours) near and below the glass transition
temperature $\mathrm{T_{g}}$ due to the insufficient coherent photon flux at
synchrotron X-ray sources and the limited response time of detectors (?).

Direct experimental access to faster relaxation dynamics has been of
particular interest, because numerous glass forming liquids exhibit short
relaxation times (i.e., high fluidity) at high temperatures due to their high
liquid fragility. The fragility concept classifies the diverse variety of
liquids according to their temperature dependence of relaxation times (or
viscosity) (?). On approaching $\mathrm{T_{g}}$, some liquids exhibit a near-
Arrhenius rise in viscosity, classified as ‘strong’, while others, as ‘fragile
liquids’, show a range of super-Arrhenius behavior (?). While many liquids
follow a simple fragile or strong behavior, some anomalous liquids exhibit a
so-called fragile-to-strong transition (FST) (also referred as dynamic
crossover) (?, ?, ?). A FST is usually accompanied by thermodynamic response
function maxima (e.g. heat capacity and compressibility), as well as local
structural changes (?, ?, ?). Such a transition has been long debated in water
(?, ?, ?, ?) and suggested in silicon (?, ?, ?, ?), germanium (?), and many
others (?, ?, ?, ?, ?). A clear FST has been demonstrated in liquid
$\mathrm{Ge_{15}Te_{85}}$ near its eutectic melting point (?). In the related
systems (e.g. $\mathrm{Ge_{15}Sb_{85}}$ and $\mathrm{AgInSbTe}$), the FST has
been thought to play an important role in the functionality of phase-change
memory devices (?, ?, ?). Yet, the atomic-level understanding of relaxation
dynamics near these transitions presents a tremendous challenge, because 1)
FSTs typically occur far above $\mathrm{T_{g}}$ with a short relaxation
timescale from nanoseconds to sub-picoseconds, far beyond the capability of
synchrotron-based XPCS, and 2) FSTs are often hidden in the supercooled liquid
below the melting temperature $\mathrm{T_{m}}$ obscured by fast
crystallization.

X-ray free electron lasers (XFELs), delivering unprecedentedly high photon
fluxes within sub-100-femtosecond pulses with nearly full transverse
coherence, present the opportunities for developing the XPCS techniques being
capable to probe dynamics in the regime from nanoseconds to femtoseconds
timescales. The key idea is that a femtosecond X-ray pulse is split into two
nearly identical coherent pulses with an adjustable time delay in-between. As
the pulse pairs scatter from the sample in rapid succession, the summed
speckle patterns are collected by a 2D detector, which carry the information
of atomic dynamics on the timescale of the pulse separation. However,
implementing this concept has faced technical challenges including generating
identical double pulses, extremely low count rates at wide angles (e.g.,
$<10^{-4}$ photons/pixel), and lack of robust analytical methods of extracting
dynamics from noise or artifacts.

Pioneering work has explored the feasibility of split-delay optics of angstrom
wavelength to deliver the identical double pulses (?, ?, ?, ?, ?, ?, ?). The
latter is crucial for extracting the dynamics information of the sample. The
conventional division-of-wavefront split line, despite providing double pulses
with good efficiency, suffers from the instability of the crystal-optics-based
beam splitters and does not preserve sufficient mutual coherence between the
two pulses (?). A more recent approach, introduced by some co-authors of this
work, employs an amplitude-splitting delay line using transmission grating-
based splitters, which generated highly mutually coherent hard-X-ray pulse
pairs (?). In addition, the femto- to pico-seconds separation of pulse pairs
is beyond the time resolution of any X-ray detectors, which renders the data
analysis strategy different from the synchrotron-based XPCS. Theories, relying
on the analysis of photon statistics of the summed speckles, have been
developed for obtaining the ISF (?). Earlier studies reported measurements on
nanometer-scale dynamics with small angle scattering (?, ?). However, the
access to atomic-scale dynamics has been difficult due to the orders of
magnitude lower photon flux at wide angle scattering. A recent study of water
observed the speckle contrast decay with a split-delay line, possibly
reflecting water’s atomic-scale dynamics; however, it is unclear whether the
two-pulse overlap remains stable and how the overlap perturbs the measured
dynamics (?).

In this work, we show that, using the double X-ray pulses in rapid succession
delivered with the compact split-delay optics system, we can achieve
femtosecond temporal and atomic-scale spatial resolutions with X-ray photon
correlation spectroscopy (fs-XPCS). We demonstrate a direct determination of
the ISF in the sub-picosecond timescale in a model system
$\mathrm{Ge_{15}Te_{85}}$ above its FST in the fragile state. We have
developed a method of simultaneously measuring the SAXS and WAXS signals for
the overlap corrections of pulse pairs. Even with extremely low photon counts,
contrast extraction with photon statistics analysis allows for direct
observation of speckle contrast decay resulting from sample dynamics.

##  2 Results

Figure 1 shows the experimental schematic of the fs-XPCS setup. The output
pulse pairs from the compact split-delay system (?) with a photon energy of
9.5 keV were focused on the sample location by a set of compound refractive
beryllium lenses with a focal length of 1.5 m. The sample,
$\mathrm{Ge_{15}Te_{85}}$, was encapsulated in a quartz capillary with an
inner diameter of about 10 $\mu$m, and resistively heated to the fragile
liquid state at 550 $\mathrm{{}^{\circ}C}$ (well above the eutectic melting
point of 380 $\mathrm{{}^{\circ}C}$). The average focal spot size is
approximately 2 $\mu$m (FWHM) and the average pulse pair energy is
characterized to be 0.15 $\mu$J at the sample plane. Four ePix100 detectors
(?) were assembled in a 2 by 2 array and placed 2.5 m downstream to cover the
first main structure factor peak (at Q$=2~{}\mathrm{\AA^{-1}}$) of the wide-
angle X-ray scattering (WAXS) of the sample to probe its atomic-scale
relaxation dynamics. Another ePix100 detector was mounted 5 m downstream to
simultaneously measure the small angle scattering (SAXS), which came mostly
from the quartz capillary and allowed for in-situ quantifying of the spatial
overlap of the two pulses.

Figure 1: The schematic of the experimental setup of the femtosecond X-ray
photon correlation spectroscopy (fs-XPCS). The grating-based split-delay
optics are schematically shown on the left panel. An array of WAXS and SAXS
detectors, placed downstream of the sample, collects speckle patterns at wide
and small angles simultaneously.

Ensuring the satisfactory condition of pulse pairs is essential for
distinguishing intrinsic sample dynamics from many possible artifacts. Figure
1 (left panel) shows a schematic of the split-delay optics. Each XFEL pulse is
split into two pulses via a diamond transmission grating (?). They are
directed by Bragg crystal reflections to travel in different optical paths
with the path length difference determining their time delay $\Delta t$. The
two pulses are then recombined to the same trajectory through the other
downstream diamond grating. Transmissive intensity diagnostics consisting of a
diode collecting the scattering signal from a Kapton target are placed in the
optical paths of individual branches to measure their intensities on a shot-
to-shot basis. By installing shutters in the respective optical paths for the
two beams, we can constantly cycle between three modes of pulses: 1) the
single pulse per shot through the path length adjustable branch (i.e. delayed
branch), 2) the single pulse per shot through the fixed path branch (i.e. fix-
delay branch), and 3) the pulse pairs per shot through both branches. As a
result, speckle patterns were collected at both SAXS and WAXS detectors for
each mode (i.e., delayed, fix-delay, both). By comparing the speckle contrast
of SAXS-signals from the three measurements, we are able to monitor the degree
of transverse coherence as well as the spatial overlap of the two pulses
continuously. This information provides the feedback to the analysis of the
wide-angle visibility degradation, mandatory for the accurate extraction of
ISF for sample dynamics (as discussed below).

Figure 2: The small-angle coherent scattering from the static reference for
obtaining the real-time overlap coefficient of the pulse pairs. The speckle
patterns are averaged over 40 frames for the delayed (A), fix-delay (B) and
both (C) branch modes. They were measured while collecting the 0.1 ps wide-
angle data. The displayed region corresponds to the square highlighted in
magenta in the SAXS detector shown in Fig. 1. The gray scales show the count
rate and have a unit of photons/pixel. The scale bar in (A) is the same for
all 3 plots and corresponds to 0.0005 $\mathrm{\AA^{-1}}$. (D) The speckle
contrasts (upper panel), the branching ratio $r$ of the pulse pairs, and the
overlap coefficient $\mu$ (lower panel) obtained from the intensity diagnostic
and the SAXS data with the elapse of the measurement time. During the same
period, the wide-angle speckle data for 0.1 ps and 0.65 ps delays were
collected simultaneously. The magenta arrow points to the time when the
speckles in (A-C) are measured.

The contrast of the summed speckle pattern from both beams
($\beta_{\mathrm{both}}$) can be related to that of the single branch beam
($\beta_{\mathrm{fix-delay}}$ and $\beta_{\mathrm{delayed}}$) with

$\beta_{\mathrm{both}}=r^{2}\beta_{\mathrm{fix-
delay}}+(1-r)^{2}\beta_{\mathrm{delayed}}+r(1-r)(\beta_{\mathrm{fix-
delay}}+\beta_{\mathrm{delayed}})\mu|f(\mathrm{Q},\Delta t)|^{2},$ (1)

where $f(\mathrm{Q},\Delta t)$ is the ISF corresponding to the scattering
wavevector Q and a time delay $\Delta t$. $r=i_{\mathrm{fix-
delay}}/(i_{\mathrm{fix-delay}}+i_{\mathrm{delayed}})$ denotes the branching
ratio, i.e., the fix-delay branch beam intensity fraction in the pulse pair
(see Supplementary Material Section 1 for details regarding the intensity
distribution of the split-delay system). $\mu$ ($0\leq\mu\leq 1$), known as
the beam overlap coefficient (?), characterizes the effect due to the
imperfect spatial overlap between the two branches. $\mu=1$ corresponds to an
ideal overlap, while $\mu=0$ instantly decouples contrast values from the ISF
determination. Since both $\mu$ and ISF contribute to the change in
$\beta_{\mathrm{both}}$, extracting the contrast changes alone does not yield
the ISF. As a result, a sufficiently high overlap level and a real-time
diagnostic of the overlap coefficient are a prerequisite for accurately
extracting the ISF.

In the experiment, the SAXS intensity originates from the scattering of the
capillary quartz glass, which serves as a static reference, i.e.,
$|f(\mathrm{Q},\Delta t)|=1$, assuming no picosecond dynamics of the quartz
glass at this low-Q region on such a short timescale. We are therefore able to
quantify $\mu$ directly from the single and both beam SAXS contrasts. Figure 2
(A-C) show a region in the averaged speckle patterns of the SAXS data
corresponding to the three measurement modes. The high similarity of the three
patterns indicate an excellent spatial overlap between the two pulses. Their
contrasts as a function of the measurement time are shown in the upper panel
of Fig. 2 (D). Since the transmissive intensity diagnostic provides a real
time branching ratio measurement as shown in the middle panel, the overlap
coefficient $\mu$ can be derived using Equation 1 and is shown in the lower
panel. One sees $\mu$ drifts on a minute timescale. The value of $\mu$ remains
mostly above 0.6, indicating a satisfactory overlap, thanks to the real-time
monitoring and re-alignment of the beams, although $\mu$ of some frames may
occasionally drift well below 0.6 (e.g. on Jun-05 9:00-9:30 as in Fig. 2). For
the consistency of analyses, we have excluded the WAXS data frames with
$\mu<0.6$. Furthermore, an effective overlap value $\bar{\mu}$ is derived from
the averaged contrast values of the 3 modes during the exact time span of the
WAXS measurements. This value is then used to benchmark the WAXS contrast
reduction and compute the ISF for each time delay points, as detailed in Table
1.

$\Delta t$ (ps) |  $\bar{k}$ (photons/pixel) | $N_{\mathrm{frame}}$ | $N_{\mathrm{pixel}}$ | $\hat{\beta}$ | $\delta\beta_{\mathrm{MLE}}$ | $\delta\beta_{\mathrm{Poisson}}$ | $\bar{\mu}$ | $\bar{r}$ | $|f|^{2}$
---|---|---|---|---|---|---|---|---|---
0.1 | 5.8 $\times 10^{-5}$ | 257250 | 657598 | 0.174 | 0.046 | 0.046 | 0.80 (0.03) | 0.51 (0.05) | 1.3 (0.9)
0.65 | 4.6 $\times 10^{-5}$ | 283067 | 723038 | 0.042 | 0.049 | 0.049 | 0.76 (0.08) | 0.59 (0.07) | -0.7 (0.8)

Table 1: Measurement statistics, parameters and results including the count
rate $\bar{k}$, the number of frames $N_{\mathrm{frame}}$ and pixels
$N_{\mathrm{pixel}}$, the speckle contrast levels $\hat{\beta}$, their error
bars from both the maximum likelihood fitting $\delta\beta_{\mathrm{MLE}}$ and
the photon statistics $\delta\beta_{\mathrm{Poisson}}$ for the two delays,
i.e., 0.1 ps and 0.65 ps. The effective overlap $\bar{\mu}$ from the SAXS
monitoring and the averaged branching ratio $\bar{r}$ from the intensity
diagnostic, corresponding to the same period of WAXS data collection, are used
to calculate the ISF, $|f|^{2}$.

The wide-angle speckle patterns collected at the first structure factor
maximum ($\mathrm{Q=2~{}\AA^{-1}}$) contain the information of relaxation
dynamics of the fragile liquid state. Since the scattering cross-section from
the atomic-scale order is limited, the count rate is on average about $5\times
10^{-5}$ photons/pixel with both beams illuminating the sample. It therefore
requires a large number of frames from detectors containing multi-mega pixels
and the careful treatment of artifacts including abnormal detector pixels,
background radiation from impurities in concrete and cosmic rays (?), etc, to
accurately extract the contrasts from photon statistics, i.e., from the
probability of multiple photon per pixel events (see Supplementary Material
Section 2). The method of obtaining contrast from speckle patterns with
discrete photon events has been demonstrated by Hruszkewycz et, al in (?),
albeit at a much larger count rate, on average more than 0.01 photons/pixel.
At this count rate, beam-induced permanent changes are visible when the sample
was evaluated with a scanning electron microscope. It remains experimentally
unclear whether the visibility extraction method can be further applied to
even lower count rates in the non-perturbed regime, by simply extending the
measurement time. Moreover, since monochomatic SASE (self-amplified
spontaneous emission) pulses from an XFEL exhibits large intensity
fluctuations (?), it is essential to appropriately weight the contribution of
each individual frame to the overall photon statistics. Therefore, a maximum
likelihood based contrast estimator (MLE) is introduced. The likelihood ratio
$\chi^{2}$ is defined as (?, ?)

$\chi^{2}(\beta)=-2\sum_{f=1}^{N_{\mathrm{frame}}}\sum_{k=0}^{n_{k}}p_{f,k}N_{\mathrm{pixel}}\ln(\frac{P_{k}(\beta,\bar{k_{f}})}{p_{f,k}})$
(2)

For each frame $f$ with $N_{\mathrm{pixel}}$ available pixels, the
probabilities of each pixel having photons $k$ = 0, 1, 2 (we only consider up
to $n_{k}=2$ photon/pixel events) are calculated and denoted as $p_{f,k}$,
which are then compared with the projected theoretical probability $P_{k}$
given the count rate $\bar{k_{f}}$. The optimal estimate of the contrast,
$\hat{\beta}$, is obtained by finding the minimum $\chi^{2}$ values from
numerical calculations. The error $\delta\beta_{\mathrm{MLE}}$ can be
retrieved by computing the second derivative of $\chi^{2}$ at $\hat{\beta}$

$\delta^{2}\beta_{\mathrm{MLE}}=(\frac{\partial\chi^{2}(\beta)}{2!\partial\beta^{2}})^{-1}|_{\beta=\hat{\beta}}.$
(3)

This method is used to determine the contrast values for all WAXS
measurements, including the two single-pulse modes, (i.e., fix-delay and
delayed branch) as well as the both-beam mode at 0.1 ps and 0.65 ps, as shown
in Fig. LABEL:fig:fig3 (A-D). Strikingly, as shown in Table 1, the MLE error
$\delta\beta_{\mathrm{MLE}}$ agrees well with the error derived from the
photon counting statistics,

$\delta\beta_{\mathrm{Poisson}}=\frac{1}{\bar{k}}\sqrt{\frac{2(1+\beta)}{(1+M)N_{\mathrm{pixel}}N_{\mathrm{frame}}}},$
(4)

where $M\approx 1$ characterizes the intensity fluctuations of the double
pulses incident onto the sample (see Supplementary Material Section 2 for the
derivation). Such a high level of agreement shows that Equation 2 is an
efficient estimator and the shot noise is the main source of error even at
this extremely low count rate. This non-trivial finding suggests that the
measurement accuracy of fs-XPCS can be improved by increasing data volume. It
is also worth noting that each individual branch contrast is calculated using
more than half a million frames (Fig. LABEL:fig:fig3 (A, B)), although their
error bars are still large as the count rate is about half of that of the
both-beam mode.

With the contrasts of individual and both beams at the wide angle, the
effective overlap coefficient $\bar{\mu}$ from the SAXS monitor, and the
averaged branching ratio $\bar{r}$ from the intensity diagnostics, the ISF for
each delay can be determined using Equation 1. Figure LABEL:fig:fig3 (F) shows
the normalized ISF $|f|^{2}$ at two time delays, 0.1 and 0.65 ps (also in
Table 1). At the delay of 0.1 ps, $|f|^{2}=1.3\pm 0.9$ indicates a nearly full
correlation, as one would expect at this temperature. At 0.65 ps, the
$|f|^{2}$ falls to almost zero within the error bar, suggesting a complete
decorrelation. This indicates that the relaxation timescale is on the order of
hundreds of femtoseconds. The possibility of first-pulse induced dynamics can
be ruled out by a contrast analysis with varying pulse influx. Since the pulse
influx is approximately proportional to the count rate, the contrast data for
0.65 ps are grouped with respect to the count rate by a threshold value
$\bar{k_{0}}=2\times 10^{-4}$ photons/pixel (Fig. LABEL:fig:fig3 (E)). Between
the two groups, the average incident pulse flux differs approximately sixfold,
but the contrast values of the two groups show little change within the error
margin. This indicates that the reduction in contrast observed at the delay of
0.65 ps is not related to the X-ray pulse itself. To validate the ISF
determined by the fs-XPCS, we performed the quasi-elastic neutron scattering
experiment (QENS) at the same Q and the same temperature at the Swiss
spallation neutron source SINQ (?), Paul Scherrer Institute (see Supplementary
Material Section 4). The quasi-elastic contribution (Lorentzian) of the
dynamic structure factors is Fourier transformed into the time domain, which
is an exponential decay in the ISF (blue curve) as shown in Fig. 3 (F). The
characteristic relaxation time $\tau=0.67\pm 0.12$ ps from the QENS is well in
line with the drop of the ISF at 0.65 ps from the fs-XPCS, which confirms the
results of fs-XPCS. Given the exponential shape of the ISF, a fit to the fs-
XPCS data leads to a rough estimation of $\tau$ around 0.5 ps, which is
plotted as a red circle in Fig. LABEL:fig:fig3 (G).

Figure LABEL:fig:fig3 (G) shows the relaxation times $\tau$ (in ps) over 16
orders of magnitude against the inverse temperature $\mathrm{1000/T}$
$(\mathrm{K^{-1}})$, where the blue curve represents the predicted $\tau$ from
the Adam-Gibbs equation, $\eta=\eta_{0}\exp(C/TS_{c})$ (where $\eta$ is
viscosity, $C$ and $\eta_{0}$ constants, $S_{c}$ configurational entropy),
based on heat capacity and high-temperature viscosity data (see Supplementary
Material Section 3). A characteristic ‘double-kink’ in the temperature
dependence of viscosity (yellow area) has been identified as the FST in the
liquid (?). Above 500 $\mathrm{{}^{\circ}C}$, the liquid is in a fragile
liquid state (with a fragility index $m\sim 90$, similar to
$\mathrm{Ge_{2}Sb_{2}Te_{5}}$). The relaxation time from the fs-XPCS
experiment agrees surprisingly well with the theoretical predictions (the red
circle on the blue line).

For a comparison, we have additionally conducted synchrotron-based XPCS at
DESY/PETRA III (beamline P10) within its limit on the timescale of
milliseconds or slower (Supplementary Material Section 5). In contrast to the
XFEL-based fs-XPCS, the sample dynamics are determined by varying the
integration window of different detector exposure time with synchrotron X-ray
beams. The temporal evolution of intensity fluctuations $g_{2}(\mathrm{Q},t)$,
where $\mathrm{Q}=2~{}\mathrm{\AA^{-1}}$, is monitored, which can be related
to the ISF, $f(\mathrm{Q},t)$, via the Siegert relation
$g_{2}(\mathrm{Q},t)=1+B\cdot|f(\mathrm{Q},t)|^{2}$ with $B$ denoting the
optical contrast for the experimental setup (?). We can characterize the
relaxation dynamics in the ultraviscous liquid slightly above $\mathrm{T_{g}}$
(130 $\mathrm{{}^{\circ}C}$) by modelling the data with the
Kohlrausch–Williams–Watts equation (?),
$|f(\mathrm{Q},t)|^{2}=c\cdot\exp[(-2(t/\tau)^{\gamma}]$, where $\tau$ denotes
the relaxation time, and $\gamma$ is the shape parameter describing the nature
of the microscopic relaxation process. The resulting relaxation times $\tau$
on the second timescale are plotted as magenta crosses near $\mathrm{T_{g}}$
in Fig. LABEL:fig:fig3 (G). They are in agreement with the dynamic mechanical
spectroscopy (DMA) data (diamond symbols) (?), and consistent with the
prediction of the Adam-Gibbs theory. The $\gamma$ ranging from 0.6 to 0.9
indicates a stretched exponential decay and dynamical heterogeneity in the
strong liquid phase (Supplementary Material Section 5). Clearly, the XFEL-
based fs-XPCS extends the capability of XPCS by about 12 orders of magnitude
in timescales, enabling us to cover the entire temperature range from the
ultraviscous strong liquid near the glass formation to the fluid highly
fragile liquid well above its melting point.

##  3 Discussion and Conclusion

While XPCS has become a major experimental tool for understanding the atomic-
level relaxation dynamics near or below $\mathrm{T_{g}}$ in glass science, the
split-delay based fs-XPCS extends the boundary of the time resolution to the
femtosecond to picosecond regimes. It is worth noting that the pulse pairs can
be separated by tens of nanoseconds using the so-called nanosecond double-
bunch mode, allowing for accessing the dynamics on the nanosecond timescales
(?), as demonstrated in a colloidal system (?). Unlike the frequency-domain
analysis performed in inelastic X-ray or neutron scattering experiments, the
direct measurements in the time-domain with XPCS allows for obtaining the
information of higher order time-correlations (?). The latter is important for
temporally heterogeneous relaxations of disordered systems (e.g. intermittent
relaxations) (?, ?). In these contexts, the XFEL-based fs-XPCS using separated
pulse pairs will open a new era for studying the atomic-scale structural
relaxations in uncharted territories (obscured by extremely short time
relaxation or fast crystallization) in the glass and liquid sciences. For
instance, one could implement a laser-pulse pump shortly before the double
X-ray pulses to probe the dynamics during the rapid melt-and-quench process of
a liquid (e.g. with a cooling rate of $10^{11}$ K/s) before fast
crystallization interferes. Such a pump-and-double-probe scheme might be used
to access the relaxation dynamics of supercooled water below its homogeneous
nucleation limit 231 K (?) (termed as ‘no-man’s land), where the origin of the
well-known thermodynamic anomalies has been debated for decades (?). Poor
glass-forming phase-change materials (e.g. $\mathrm{Ge_{15}Sb_{85}}$ and
others (?, ?)) undergo a metal-to-semiconductor fragile-strong transition
during supercooling (20-30% below $\mathrm{T_{m}}$) before it crystallizes
within a few nanoseconds. In this context, the atomic-scale relaxation
dynamics near these ‘short-lived’ transitions is of particular relevance for
understanding the functionality of these materials in phase-change memory
devices.

It is worth noting that approximately 1.7 million frames are analyzed in this
study to determine the ISF at two delay points. The number of delay points is
constrained by the pulse repetition rate of 120 Hz presently available at the
LCLS, and the precision of our measurements is primarily limited by shot
noise. Thus, data collection for each delay point requires several hours.
However, this scenario will undergo a transformative change with the advent of
LCLS-II-HE, where an X-ray source with a MHz repetition rate will be available
(?). Together with the cutting-edge multi-mega-pixel X-ray detectors such as
CITIUS (?) and SParkPix-S (?), which operate at tens of kHz to MHz frame
rates, it naturally allows us to collect data at substantially higher rates
($\sim$ four orders of magnitude faster), potentially reducing the measurement
time per delay point to tens of seconds. The methodology demonstrated in this
work can then be efficiently employed to investigate the structural relaxation
over the entire temperature range where the dynamic crossovers may occur.

## Acknowledgement

Use of the Linac Coherent Light Source (LCLS), SLAC National Accelerator
Laboratory, is supported by the U.S. Department of Energy, Office of Science,
Office of Basic Energy Sciences under Contract No. DE-AC02-76SF00515. The QENS
results are based on experiments performed at the Swiss spallation neutron
source SINQ, Paul Scherrer Institute, Villigen, Switzerland. The synchrotron-
based XPCS was performed at DESY PETRAIII P10 (Hamburg, Germany). Y.S. would
like to acknowledge many inspiring discussions with Mark Sutton, from
experimental design to data analysis. The authors thank Eloi Pineda for
assistance of the DESY experiment. This work was supported by a research grant
(42116) from VILLUM FONDEN (S.W., T.F.) and partially supported by Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) (WE 6440/1-1). We
thank the Danish Agency for Science, Technology, and Innovation for funding
the instrument center DanScatt. KST and TJA acknowledges financial support by
the DFG; Project No. 278162697-SFB 1242.

### Author contributions

YS and SW conceived and led the project. YS, HL, TJA, SS, TS, KST, DZ, JM, AC,
PS, SW, MM, CY, JVC, VE, JV, NW, TM, CD, PF designed and performed the fs-XPCS
experiment with samples from NA and JM, and with furnace by YC, YS, and FY; TF
and YS analyzed the fs-XPCS data. SW, TF, PS, JM, AM, and JE performed the
QENS with samples from GC and PL. The QENS data were analyzed by TF.
Synchrotron-XPCS experiment was carried out by AC, BR, NA, FY, and MS, and
data were analyzed by AC, NA, BR, and SW; YS, SW, and TF wrote the paper with
input from all the other authors.

## References and Notes

  * 1. W. H. Wang, Dynamic relaxations and relaxation-property relationships in metallic glasses, Prog. Mater. Sci. 106, 100561 (2019).
  * 2. C. A. Angell, K. L. Ngai, G. B. McKenna, P. F. McMillan, S. W. Martin, Relaxation in glassforming liquids and amorphous solids, J. Appl. Phys. 88, 3113 (2000).
  * 3. H. Tanaka, Roles of liquid structural ordering in glass transition, crystallization, and water’s anomalies, J. Non-Cryst. Solids: X 13, 100076 (2022).
  * 4. B. Ruta, E. Pineda, Z. Evenson, Relaxation processes and physical aging in metallic glasses, J. Phys. Condens. Matter 29, 503002 (2017).
  * 5. J. C. Dyre, Colloquium: The glass transition and elastic models of glass-forming liquids, Rev. Mod. Phys. 78, 953 (2006).
  * 6. X. Monnier, D. Cangialosi, B. Ruta, R. Busch, I. Gallino, Vitrification decoupling from $\alpha$-relaxation in a metallic glass, Sci. Adv. 6, eaay1454 (2020).
  * 7. Y. Cheng, et al., Highly tunable $\beta$-relaxation enables the tailoring of crystallization in phase-change materials, Nat. Commun. 13, 7352 (2022).
  * 8. E. R. Weeks, J. C. Crocker, A. C. Levitt, A. Schofield, D. A. Weitz, Three-dimensional direct imaging of structural relaxation near the colloidal glass transition, Science 287, 627 (2000).
  * 9. R. Richert, Supercooled liquids and glasses by dielectric relaxation spectroscopy, Adv. Chem. Phys. 156, 101 (2014).
  * 10. B. Ruta, et al., Atomic-scale relaxation dynamics and aging in a metallic glass probed by x-ray photon correlation spectroscopy, Phys. Rev. Lett. 109, 165701 (2012).
  * 11. O. G. Shpyrko, X-ray photon correlation spectroscopy, J. Synchrotron Rad. 21, 1057 (2014).
  * 12. B. Ruta, et al., Wave-vector dependence of the dynamics in supercooled metallic liquids, Phys. Rev. Lett. 125, 055701 (2020).
  * 13. G. Grübel, G. Stephenson, C. Gutt, H. Sinn, T. Tschentscher, Xpcs at the european x-ray free electron laser facility, Nucl. Instrum. Methods Phys. Res. B 262, 357 (2007).
  * 14. C. A. Angell, Formation of glasses from liquids and biopolymers, Science 267, 1924 (1995).
  * 15. C. Angell, Liquid fragility and the glass transition in water and aqueous solutions, Chem. Rev. 102, 2627 (2002).
  * 16. H. Tanaka, Liquid–liquid transition and polyamorphism, J. Chem. Phys. 153 (2020).
  * 17. S. Wei, Anomalous liquids on a new landscape: From water to phase-change materials, J. Non-Cryst. Solids: X 14, 100094 (2022).
  * 18. S. Wei, P. Lucas, C. A. Angell, Phase change alloy viscosities down to tg using adam-gibbs-equation fittings to excess entropy data: A fragile-to-strong transition, J. Appl. Phys. 118 (2015).
  * 19. P. Zalden, et al., Femtosecond x-ray diffraction reveals a liquid–liquid phase transition in phase-change materials, Science 364, 1062 (2019).
  * 20. S. Wei, P. Lucas, C. A. Angell, Phase-change materials: The view from the liquid phase and the metallicity parameter, MRS Bull. 44, 691 (2019).
  * 21. P. Gallo, et al., Water: A tale of two liquids, Chem. Rev. 116, 7463 (2016).
  * 22. K. H. Kim, et al., Maxima in the thermodynamic response and correlation functions of deeply supercooled water, Science 358, 1589 (2017).
  * 23. A. Dehaoui, B. Issenmann, F. Caupin, Viscosity of deeply supercooled water and its coupling to molecular diffusion, Proc. Natl. Acad. Sci. 112, 12020 (2015).
  * 24. S. Sastry, C. Austen Angell, Liquid–liquid phase transition in supercooled silicon, Nat. Mater. 2, 739 (2003).
  * 25. N. Jakse, A. Pasturel, Liquid-liquid phase transformation in silicon: evidence from first-principles molecular dynamics simulations, Phys. Rev. Lett. 99, 205702 (2007).
  * 26. G. Zhao, et al., Phase behavior of metastable liquid silicon at negative pressure: Ab initio molecular dynamics, Phys. Rev. B 93, 140203 (2016).
  * 27. V. V. Vasisht, S. Saw, S. Sastry, Liquid–liquid critical point in supercooled silicon, Nat. Phys. 7, 549 (2011).
  * 28. M. Bhat, et al., Vitrification of a monatomic metallic liquid, Nature 448, 787 (2007).
  * 29. L. Henry, et al., Liquid–liquid transition and critical point in sulfur, Nature 584, 382 (2020).
  * 30. L. Xu, I. Ehrenberg, S. V. Buldyrev, H. E. Stanley, Relationship between the liquid–liquid phase transition and dynamic behaviour in the jagla model, J. Phys. Condens. Matter 18, S2239 (2006).
  * 31. M. Stolpe, et al., Structural changes during a liquid-liquid transition in the deeply undercooled zr58.5cu15.6ni12.8al10.3nb2.8 bulk metallic glass forming melt, Phys. Rev. B 93, 014201 (2016).
  * 32. J. Orava, D. W. Hewak, A. L. Greer, Fragile-to-strong crossover in supercooled liquid ag-in-sb-te studied by ultrafast calorimetry, Adv. Funct. Mater. 25, 4851 (2015).
  * 33. W. Roseker, et al., Performance of a picosecond x-ray delay line unit at 8.39 kev, Opt. Lett. 34, 1768 (2009).
  * 34. T. Osaka, et al., A bragg beam splitter for hard x-ray free-electron lasers, Opt. Express 21, 2823 (2013).
  * 35. T. Osaka, et al., Wavelength-tunable split-and-delay optical system for hard x-ray free-electron lasers, Opt. Express 24, 9187 (2016).
  * 36. D. Zhu, et al., Development of a hard x-ray split-delay system at the linac coherent light source, Proc. SPIE 10237, 102370R (2017).
  * 37. W. Lu, et al., Development of a hard x-ray split-and-delay line and performance simulations for two-color pump-probe experiments at the european xfel, Rev. Sci. Instrum. 89, 063121 (2018).
  * 38. R. Rysov, W. Roseker, M. Walther, G. Grübel, Compact hard x-ray split-and-delay line for studying ultrafast dynamics at free-electron laser sources, J. Synchrotron Rad. 26 (2019).
  * 39. Y. Sun, et al., Compact hard x-ray split-delay system based on variable-gap channel-cut crystals, Opt. Lett. 44, 2582 (2019).
  * 40. Y. Sun, et al., Realizing split-pulse x-ray photon correlation spectroscopy to measure ultrafast dynamics in complex matter, Phys. Rev. Res. 2, 023099 (2020).
  * 41. H. Li, et al., Generation of highly mutually coherent hard-x-ray pulse pairs with an amplitude-splitting delay line, Phys. Rev. Res. 3, 043050 (2021).
  * 42. S. Hruszkewycz, et al., High contrast x-ray speckle from atomic-scale order in liquids and glasses, Phys. Rev. Lett. 109, 185502 (2012).
  * 43. W. Roseker, et al., Towards ultrafast dynamics with split-pulse x-ray photon correlation spectroscopy at free electron laser sources, Nat. Commun. 9, 1704 (2018).
  * 44. Y. Sun, et al., Nonuniform flow dynamics probed by nanosecond x-ray speckle visibility spectroscopy, Phys. Rev. Lett. 127, 058001 (2021).
  * 45. Y. Shinohara, et al., Split-pulse x-ray photon correlation spectroscopy with seeded x-rays from x-ray laser to study atomic-level dynamics, Nat. Commun. 11, 6213 (2020).
  * 46. M. Sikorski, et al., Application of an epix100 detector for coherent scattering using a hard x-ray free-electron laser, J. Synchrotron Rad. 23, 1171 (2016).
  * 47. A. Q. Baron, D. Ishikawa, Background rates in several hard x-ray photon counting pixel array detectors, Nucl. Instrum. Methods Phys. Res. A 1049, 168101 (2023).
  * 48. D. Zhu, et al., Performance of a beam-multiplexing diamond crystal monochromator at the linac coherent light source, Rev. Sci. Instrum. 85 (2014).
  * 49. S. S. Wilks, The large-sample distribution of the likelihood ratio for testing composite hypotheses, Ann. Math. Statist. 9, 60 (1938).
  * 50. S. Janßen, J. Mesot, L. Holitzner, A. Furrer, R. Hempelmann, FOCUS: a hybrid TOF-spectrometer at SINQ, Phys. B: Condens. Matter. 234, 1174 (1997).
  * 51. S.-X. Peng, et al., Uncovering $\beta$-relaxations in amorphous phase-change materials, Sci. Adv. 6, eaay6726 (2020).
  * 52. F. Decker, et al., Two bunches with ns-separation with lcls, Proc. FEL2015 WEP023, 634 (2015).
  * 53. A. Madsen, R. L. Leheny, H. Guo, M. Sprung, O. Czakkel, Beyond simple exponential correlation functions and equilibrium dynamics in x-ray photon correlation spectroscopy, New J. Phys. 12, 055001 (2010).
  * 54. Z. Evenson, et al., X-ray photon correlation spectroscopy reveals intermittent aging dynamics in a metallic glass, Phys. Rev. Lett. 115, 175701 (2015).
  * 55. H. Kanno, C. A. Angell, Water: Anomalous compressibilities to 1.9 kbar and correlation with supercooling limits, J. Chem. Phys. 70, 4008–4016 (1979).
  * 56. K. H. Kim, et al., Experimental observation of the liquid-liquid transition in bulk supercooled water under pressure, Science 370, 978–982 (2020).
  * 57. T. Raubenheimer, et al., The lcls-ii-he, a high energy upgrade of the lcls-ii, 60th ICFA Advanced Beam Dynamics Workshop on Future Light Sources MOP1WA02, 6 (2018).
  * 58. H. Nishino, et al., The false beat signal of a high-speed x-ray imaging detector for synchrotron radiation experiments and its elimination with a synchronized citius detector, Nucl. Instrum. Methods Phys. Res. A 1057, 168710 (2023).
  * 59. L. Rota, SparkPix-ED: a readout ASIC with 1 MHz frame-rate for rare event experiments at LCLS-II (2023). 24th International Workshops on Radiation Imaging Detectors.


# The Optimal Growth Mode in the Relaxation to Statistical Equilibrium

Manuel Santos Gutiérrez<EMAIL_ADDRESS>Department of
Earth and Planetary Sciences, Weizmann Institute of Science, Rehovot 76100,
Israel  Mickaël D. Chekroun  Department of Earth and Planetary Sciences,
Weizmann Institute of Science, Rehovot 76100, Israel  Department of
Atmospheric and Oceanic Sciences, University of California, Los Angeles, CA
90095-1565, USA

###### Abstract

Systems far from equilibrium approach stability slowly due to “anti-mixing”
characterized by regions of the phase-space that remain disconnected after
prolonged action of the flow. We introduce the Optimal Growth Mode (OGM) to
capture this slow initial relaxation. The OGM is calculated from Markov
matrices approximating the action of the Fokker-Planck operator onto the phase
space. It is obtained as the mode having the largest growth in energy before
decay. Important nuances between the OGM and the more traditional slowest
decaying mode are detailed in the case of the Lorenz 63 model. The
implications for understanding how complex systems respond to external forces,
are discussed.

††preprint: APS/123-QED

In many chaotic and complex systems, mixing plays a crucial role. It describes
how the system loses memory of its initial state and reaches a well-defined,
statistical steady state. Correlation functions decay then over time during
mixing. This decay rate determines two key properties: 1) how quickly the
system reaches a statistical steady state [1], and 2) its resilience against
external disturbances [1, 2].

For systems at equilibrium, the approach to a statistical steady state is
remarkably fast (exponential), establishing a direct link between mixing,
correlations, and the system’s response to external forces [3, 4]. However,
this connection is not guaranteed for systems out of equilibrium. This is
because these systems lack smooth invariant probability measures [5, 6].
Additionally, the introduction of non-conservative forces violates Liouville’s
theorem, a cornerstone of statistical physics [7, §I]. As a result, the
governing equation (Liouville equation [7]) may lose a key property (self-
adjointness), severing the link between the system’s long-term and short-term
behavior due to the presence of non-normal modes. This phenomenon as been well
documented in the fluid dynamics literature [8, 9, 10] but its consequences on
the action of Liouville propagators—also known as transfer or Frobenius-Perron
operators [11, 12]—in the relaxation to statistical equilibria, has been much
less discussed. This paper proposes to bridge this gap.

By hinging on the Ulam method [13], Markov matrices are instrumental in the
learning of transfer operator’s properties; see e.g. [11, 14, 15]. Markov
matrices are also instrumental in analyzing mixing [16, 17]. It is thus not
surprising that they are used for estimating the long-term behavior of the
system’s Liouville propagator characterizing how probability densities evolve
and revealing how different aspects of the system, including correlations,
change over time; see e.g. [18, 19, 20]. This approach can be justified for
both deterministic and stochastic systems [21, 22, 23].

Essentially, Markov-Ulam matrices aim at capturing the system’s transition
probabilities by identifying key regions in its phase space (possible states)
and the specific modes (behaviors) that are slowest to lose memory of their
initial state [24, 25]. They also pinpoint areas most susceptible to external
disturbances [26, 27, 28]. Within a Lagrangian setting, these matrices can
even identify persistent patterns and barriers to movement within the system’s
physical space itself [29, 30].

Despite their relevance to study nonlinear dynamics, the non-uniform
exponential decay to equilibrium of Markov-Ulam matrices has been poorly
studied. For systems out of equilibrium, these matrices, reflecting the
Liouville propagator [11, 19], are expected to be also non-selfadjoint. As a
consequence, non-normal modes are expected to occur: modes that are slow to
changes over finite time-intervals and which might not correspond to the
slowest (normal) modes to decay in the long run. This paper aims to
characterize these modes by introducing a new concept: the Optimal Growth Mode
(OGM). In simpler terms, this mode is the one that takes the longest to lose
memory and reach equilibrium. As shown below, the OGM is characterized by some
sort of “anti-mixing behavior.”

To avoid technical difficulties, we consider deterministic systems embedded
within a noise bath. This way, the existence of a unique (and smooth) ergodic
statistical steady state can be guaranteed for a broad class of systems [31,
23]. Thus, our dynamical systems, acting on the $d$-dimensional phase-space
$\mathbb{R}^{d}$, are described by a vector-field
$\mathbf{F}:\mathbb{R}^{d}\longrightarrow\mathbb{R}^{d}$ and the following
stochastic differential equation (SDE):

$\dot{\bm{x}}=\mathbf{F}(\bm{x})+\sqrt{2\epsilon}\dot{\bm{W}}_{t},$ (1)

where $\dot{\bm{W}}_{t}$ is an $d$-dimensional white-noise process and
$\varepsilon$ is a measure of its covariance. The evolution of probability
density functions associated with the SDE (1), is given by the Fokker-Planck
equation (FKE) [32]:

$\partial_{t}\rho=\mathcal{L}_{FP}\rho=-\nabla\cdot\left(\mathbf{F}\rho\right)+\epsilon\Delta\rho.$
(2)

Its solution operator is given by the Fokker-Planck propagator
$e^{t\mathcal{L}_{FP}}$ assumed to form a $C_{0}$-semigroup [33]. When
$\epsilon=0$, it reduces to the Liouville propagrator [19].

Under generic conditions on $\mathbf{F}$, Eq. (2) has an ergodic statistical
steady state $\rho_{s}$ satisfying $\mathcal{L}_{FP}\rho_{s}=0$ and the
Fokker-Planck operator $\mathcal{L}_{FP}$ has an isolated spectrum constituted
of damped eigenvalues separated from the eigenvalue $0$ (spectral gap); see
e.g. [23]. In this case, the underlying damped eigenvalues of the Fokker-
Planck operator $\mathcal{L}_{FP}$ are referred to as the Ruelle-Pollicot (RP)
resonances and imply an asymptotic exponential decay of correlations [23].
These resonances hold significant importance. They allow for decomposition
formulas for both temporal correlation functions and power spectral density
functions [23, Sec. 2.3] across a wide range of stochastic systems encountered
in physics; see also [34, 35] for deterministic systems. The spectral gap
ensures also the exponential convergence of $e^{t\mathcal{L}_{FP}}\rho_{0}$
towards $\rho_{s}$ (in some $L^{2}$ space) for a broad class of initial
distribution $\rho_{0}$; see e.g. [23, Theorem 6]. Exponential rate of
convergence and spectral gaps also hold in Wasserstein spaces for an even
broader class of models than Eq. (1); see [36, 37]. However the study of
transient behavior before reaching exponential convergence remains largely
unexplored from an ergodic viewpoint.

In classical fluid dynamics the role of transient dynamics has been thoroughly
examined. Certain flows experience a burst of transient energy growth before
settling into an exponential decay towards steady state, whose rate is
controlled by the damped eigenvalues of the linearized governing equations.
This temporary spike is not bounded by these eigenvalues and is caused instead
by non-normal modes of the underlying linearized operator [8, 9, 38, 10]. Non-
normal modes, when influenced by random fluctuations have shown to be highly
relevant. They can explain for instance the variability of large-scale climate
phenomena like the El Niño-Southern Oscillation [39, 40], ocean circulation
patterns [41, 42], or atmospheric flows [43, 44].

This paper focuses on non-normal modes associated with Fokker-Planck operators
$\mathcal{L}_{FP}$, for which exponential relaxation towards a (unique)
ergodic statistical equilibrium holds. We are particularly interested to
clarify how these modes affect this relaxation process at early times.

An operator, $\mathcal{L}$, is normal if it commutes with its adjoint,
$\mathcal{L}^{\ast}$ [45, §I.6], namely if
$\mathcal{L}^{\ast}\mathcal{L}=\mathcal{L}\mathcal{L}^{\ast}$. Note that self-
adjoint operators are a particular class of normal operators. From Liouville’s
theorem for Hamiltonian systems, the operator $\mathcal{L}_{FP}$ is normal and
more generally, for stochastic dynamics within a potential $V$, such that
$\mathbf{F}(\bm{x})=-\nabla V(\bm{x})$ in Eq. (1), it is also true that the
Fokker-Planck operator $\mathcal{L}_{FP}$ is normal [46]. Operators for which
$\mathcal{L}^{\ast}\mathcal{L}\neq\mathcal{L}\mathcal{L}^{\ast}$ are called
non-normal, for which an extensive literature exists [10] albeit not
specifically dedicated to the case of Fokker-Planck operators.

Analytical formulas of RP resonances and modes of the Fokker-Planck operator
is restricted to very few instances; see e.g. [47, 48, 49]. In practice, such
resonances are numerically estimated from (long) time series
$\\{\bm{x}_{k}\\}_{k=1}^{T}$ obtained by solving Eq. (1) [50, 23]. As
mentioned above, Ulam’s method grants a means to obtain estimates of the
dynamics’ propagator $e^{t\mathcal{L}_{FP}}$ via Markov matrices. The
eigenvalues of $\mathcal{L}_{FP}$ are then obtained through logarithm formulas
[50, 23]. This method consists of subdividing the phase space $\mathcal{X}$
(typically shadowing the dynamics’ attractor) into $N$ non-intersecting
_cells_ or _boxes_ $\\{B_{i}\\}_{i=1}^{N}$ and estimating the dynamics’
probability transitions across these boxes. Mathematically, the propagator
$e^{\tau\mathcal{L}_{FP}}$ for a given transition time $\tau$, is approximated
by a $N\times N$ Markov matrix $M_{\tau}$, whose entries are given by [23]:

$\left[M_{\tau}\right]_{i,j}=\frac{1}{\rho_{0}(B_{i})}\int_{B_{i}}e^{\tau\mathcal{L}_{FP}}\chi_{B_{j}}(\bm{x})\rho_{s}(\mathrm{d}\bm{x}),$
(3)

for $i,j=1,\ldots,N$, where $\chi_{B_{i}}$ denotes the indicator function of
the box $B_{i}$, $\rho_{s}$ the stationary statistical equilibrium, and
$\rho_{0}$ the initial distribution. The transition matrix is then estimated
by computing a classical maximum likelihood estimator given by [50, 11]:

$\left[M_{\tau}\right]_{i,j}=\frac{\\#\bigg{\\{}\Big{(}\bm{x}_{k}\in
B_{j}\Big{)}\land\Big{(}\bm{x}_{k+\ell}\in
B_{i}\Big{)}\bigg{\\}}}{\\#\Big{\\{}\bm{x}_{k}\in B_{j}\Big{\\}}},$ (4)

for $i,j=1,\ldots,N$. Given a sampling time $\delta t$ at which the time
series $\bm{x}_{k}$ is collected, the formula (4) counts the number of
observations ($\\#\\{\cdot\\}$) visiting the box $B_{i}$ after
$\ell=\lfloor\tau/\delta t\rfloor$ iterations, knowing that it already landed
in $B_{j}$. The resulting operation results into a coarse-graining of the
dynamics encoded by $M_{\tau}$ [19, 23] that incorporates artificial diffusion
[20, 51]. Note that for high-dimensional systems, time series constructed from
reduced coordinates (observables) are often used instead of the full-state
vector $\bm{x}_{k}$. In this case, transition matrices can still be estimated
in a reduced state space and their eigenvalues are still informative about the
genuine system’s transitions (in the full state space) through the notion of
reduced RP resonances; see [19, Theorem A] and [23, 52].

Either estimated in the full or a reduced state space, the transition matrix
$M_{\tau}$ satisfies the Perron-Frobenius theorem [16]. As such, its
eigenvalues are contained in the unit ball and there exists an eigenvector
$\mathbf{u}$ in $\mathbb{R}^{N}$ such that $M_{\tau}\mathbf{u}=\mathbf{u}$.
This vector approximates the ergodic invariant measure of the system. For many
stochastic systems, the resulting stochastic Matrix $M_{\tau}$ is irreducible
[23], ensuring its dominant eigenvalue $\mu_{0}$ equal to $1$ to be simple,
and the following limit to hold:

$\lim_{k\rightarrow\infty}\left\|(M_{\tau}-E_{1})^{k}\right\|_{2}=0,$ (5)

where the matrix $E_{1}$ is the spectral projector associated with the
eigenvalue $\mu_{0}$ given by $E_{1}=\mathbf{u}\mathbf{1}^{\mathrm{tr}}$,
where $\mathbf{1}$ is the $N$-dimensional (column) vector whose entries are
equal to $1$.

Then, the asymptotic rate of convergence to equilibrium is determined by the
spectral gap $\gamma$ defined as the distance between the subdominant
eigenvalue $\mu_{1}$ and the unit circle, namely $\gamma=1-|\mu_{1}|$, see
[16]. In particular, $\gamma$ controls the geometric rate of convergence of
the matrix $M_{\tau}$ towards the spectral projector $E_{1}$ in Eq. (5).

Such a geometric convergence towards zero of
$\delta(k)=\|(M_{\tau}-E_{1})^{k}\|_{2}$ is often characteristic of
exponential decay of correlations for the original system; see e.g. [23, Sec.
4]. This picture is nevertheless tempered by the potential non-normality of
the matrix $M_{\tau}-E_{1}$. In such a case, the relaxation towards zero of
$\delta(k)$ takes place only after a few iterations, with some possible
transient growth. During this transient, the system is slow to lose memory and
does not exponentially decorrelate in time.

It remains to identify what regions of the phase-space—as approximated by the
dynamics’ coarse-graining tied to Eq. (4)—display such a transient growth. The
resulting mode capturing such regions is expected to indicate sets of initial
conditions which remain disconnected for a prolonged interval of time before
mixing takes over.

To calculate such a mode, and inspired by the aforementioned work of [9] in
fluid dynamics, we seek the directions of maximum stretching upon the action
of the relaxation matrix $Q_{k}=(M_{\tau}-E_{1})^{k}$. In that respect, we
consider its singular value decomposition (SVD):

$(M_{\tau}-E_{1})^{k}=U_{k}\Sigma_{k}V_{k}^{\ast},$ (6)

where $U_{k}$ and $V_{k}$ are unitary $N\times N$ matrices and $\Sigma_{k}$ is
an $N\times N$ diagonal matrix containing the singular values
$\\{\sigma_{k,i}\\}_{i=1}^{N}$, ordered as:
$\sigma_{k,1}\geq\sigma_{k,2}\geq\ldots\geq\sigma_{k,N}$. Thus,
$\|Q_{k}\|_{2}=\sigma_{k,1}$ and there exists a vector $\mathbf{v}_{k,1}$ in
$\mathbb{R}^{N}$ such that:

$\|Q_{k}\|_{2}=\|Q_{k}\mathbf{v}_{k,1}\|_{2}=\sigma_{k,1}.$ (7)

We define then the Optimal Growth Mode (OGM) as the mode that maximizes
$\|Q_{k}\|_{2}$ over $k$. This mode, denoted by $\mathbf{w}_{OGM}$ below, is
actually given by the vector $\mathbf{v}_{k_{0},1}$ that satisfies:

$\max_{k\geq
0}\|Q_{k}\|_{2}=\|(M_{\tau}-E_{1})^{k_{0}}\mathbf{v}_{k_{0},1}\|_{2}=\sigma_{k_{0},1}.$
(8)

In other words the OGM maximizes the relaxation curve

$k\mapsto\|Q_{k}\|_{2}.$ (9)

The OGM must not be confused though with the slowest decaying mode (SDM),
defined as the eigenvector $\mathbf{w}_{SDM}$ of $\mathbb{R}^{N}$ that
satisfies:

$M_{\tau}\mathbf{w}_{SDM}=\mu_{1}\mathbf{w}_{SDM}.$ (10)

Since
$\lim_{k\rightarrow\infty}\|M_{\tau}^{k}\mathbf{w}_{SDM}\|_{2}=\lim_{k\rightarrow\infty}\mu_{1}^{k}=0$
and $\mu_{1}$ is the subdominant eigenvalue of $M_{\tau}$, the SDM is the
eigenmode of $M_{\tau}$ that has the slowest decay at asymptotic times. It
indicates regions of phase-space of slow asymptotic mixing. It remains to
illustrate how these two modes, the OGM and the SDM, reflect different
dynamical properties of chaotic systems.

_The Lorenz 63 system.—_ The Lorenz 63 (L63) model [53] is a low-dimensional
truncation of the equations of atmospheric convection [54], whose dynamical
properties have been extensively studied [55]. In its standard chaotic regime,
statistical properties have been also analyzed [56, 57]. Linear response has
been shown to hold (on numerical grounds), i.e. small changes in parameter
forcing equate to proportional changes in the statistical moments [58, 59,
60]. The L63 system writes as follows:

$\dot{\mathbf{x}}(t)=\mathbf{F}(\mathbf{x})=\begin{cases}s(y-x)\\\ x(r-z)-y\\\
xy-bz\end{cases},$ (11)

where the model parameters are set to be $s=10$, $r=28$ and $b=8/3$, for which
the system displays its famous “butterfly strange attractor” [61].

To estimate the Markov matrix $M_{\tau}$ associated with Eq. (11), we
integrate it over $10^{5}$ time units (after spinup) with a time step of
$\delta t=10^{-3}$ time units, using a fourth order Runge-Kutta algorithm.
Thus, a total of $T=10^{8}$ data points in phase-space is recorded. The
trajectory $\\{(x_{k},y_{k},z_{k})\\}_{k=0}^{T}$ shadows an attractor
contained in a compact subset of $\mathbb{R}^{3}$ [55], which lies in the
domain $[-20,20]\times[-30,30]\times[0,50]$. We work then within the reduced
state consisting of the $(x,z)$-domain, $[-20,20]\times[0,50]$. This domain is
divided into $N=2^{12}$ equally sized, non-intersecting boxes
$\\{B_{i}\\}_{i=1}^{N}$ of dimensions $40/2^{6}\times 50/2^{6}$. Then, the
transition time $\tau=10\cdot\delta t$ is selected and the transition matrix
$M_{\tau}$ given in Eq. (4) is computed with $\bm{x}_{k}=(x_{k},z_{k})$. While
the L63 system is deterministic, the coarse box-resolution yields artificial
diffusion [51] analogous to the addition of small noise to Eq. (11).

In the reduced state space focusing only on variables $x$ and $z$, the
eigenvalues of the estimated matrix $M_{\tau}$ correspond to simplified
versions of resonances of the full transfer operator (reduced RP resonances
[23]). These reduced resonances are linked to the conditional expectation [23]
in the equation closure of Eq. (11) within this reduced state space [62, 63].
This conditional expectation represents the behavior of the system after
averaging out unobserved variables. However, it may not capture all important
aspects of the system’s variability [64]. As shown in the appendix, the
spectrum of the reduced Markov matrix $M_{\tau}$ allows us here to approximate
the L63 model’s temporal variability, as reflected by its correlation
functions, with reasonable accuracy. This suggests that the modes obtained
from the reduced transition matrix $M_{\tau}$ capture, in the reduced state
space, actual dynamical properties of the L63 model. This is confirmed as
explained below.

Figure 1: Panel (A): the black curve shows the relaxation curve (9). The red
and blue curves are the successive evaluation of the matrix $M_{\tau}-E_{1}$
applied on the OGM and SDM, respectively. Panel (B): OGM calculated using the
SVD of the reduced transition matrix $M_{\tau}$. Panel (C): SDM (modulus)
calculated from the eigen-decomposition of the reduced transition matrix
matrix $M_{\tau}$. Figure 2:  Evolution of the OGM. The A-panels show
$(M_{\tau}-E_{1})^{k}\mathbf{w}_{OGM}$, for $k=0,14,\ldots,98$. The B-panels
show the evolution of initial conditions seeded around the positive and
negative areas of the OGM and evolved in time for the intervals indicated in
the titles.

The relaxation curve (Eq. (9)) is shown as the solid black line in Fig. 1 and
clearly shows the signature of non-normal transient growth and sub-exponential
decay, reaching a peak at $t\approx 0.3$ time units. Only for larger times the
leading eigenvalue $\mu_{1}$ takes over, allowing for the expected exponential
rate of decay to take place. The OGM (Fig. 1B) follows an evolution in time
(red curve in Fig. 1A) that reaches its maximum at around $0.3$ time units
while the SDM (Fig. 1C) does not show transient growth (blue curve in Fig. 1A)
although it corresponds to the slowest decaying mode at large times.

Figure 2 shows the evolution of
$\left(M_{\tau}-E_{1}\right)^{k}\mathbf{w}_{OGM}$ over time. The panels (A)
depict this for different time instants (see subtitles) corresponding to
different values of $k$ $(0,14,\ldots,98)$. Panels (B) show the fate of
$10^{3}$ red and $10^{3}$ blue trajectories over time. The red trajectories
start in areas where normalized OGM is close to $0.05$, and the blue
trajectories start in areas where it is close to $-0.05$. We can see that the
OGM, which initially concentrates near the top of the attractor, moves towards
the center and collapses before expanding outwards again. Interestingly, the
first signs of this expansion occur between $t=0.3$ and $t=0.35$, which
coincides with the peak of the relaxation curve. The trajectories (red and
blue patches in panel B) also contract and expand with the attractor, but
importantly, they do not mix with each other until $t=0.98$. This suggests
that the OGM can effectively identify regions of the attractor that mix
slowly.

In contrast to OGM, the SDM (whose absolute value is shown in Fig. 3) appears
to remain unchanged over time. Its shape is reminiscent to that of the almost
invariant set calculated in [65] for the full 3D model. This is evident not
only in panels A of Figure 3, but also in the trajectories (red and blue
patches in panel B). These trajectories mostly stay within their respective
positive and negative SDM patches. Additionally, unlike the OGM, the
trajectories colored blue and red mix much faster. By $t=0.14$, we observe for
instance particles of both colors crossing the $x=0$ line in both directions.

Both OGM and SDM exhibit properties that slow down mixing in the system, but
in different ways. Initially, the OGM experience energy growth at early times.
However, as time progresses, its convergence to a stable, long-term state
(ergodic steady state) might become faster. The SDM experiences a steady
decrease over time, following an exponential function with a rate determined
by a specific time scale given as $\tau^{-1}\log(|\mu_{1}|)$ [23, Eq. (3.40)].
However, by design, the SDM is guaranteed to decay the slowest compared to
other modes at very long times. This behavior is evident in Figure 1A for the
L63 system. Furthermore, trajectories starting in regions with opposite OGM
values mix much slower than those starting in opposite SDM regions. This is
illustrated by comparing the colored patches in Figures 2B and 3B.

Figure 3:  Evolution of the SDM. In the A-panels we show
$(M_{\tau}-E_{1})^{k}\mathbf{w}_{SDM}$, for $k=0,14,\ldots,98$. The B-panels
show the evolution of initial conditions seeded around the positive and
negative areas of the SDM and evolved in time for the intervals indicated in
the titles.

_Discussion._ — Systems out of equilibrium, characterized by non-normal
equations of motion, exhibit uneven mixing in phase space. This prevents rapid
convergence to a statistical steady state initially, even if the system
eventually decorrelates exponentially fast. The SDM (Eq. (10)) captures the
slowest asymptotic decay, but it misses the crucial short-term ”anti-mixing”
behavior common in non-normal Fokker-Planck operators (Eq. (2)). On the
contrary, the OGM captures transient dynamics. This work introduces an
algorithm based on SVD decomposition of Markov-Ulam matrices to compute the
OGM. Unlike the SDM, the OGM exhibits growth in its norm at early times,
peaking at the maximum point in the relaxation curve (Eq. (9)). This anti-
mixing character of the OGM manifests in two ways: (i) sharp increase— the OGM
norm rises significantly before decaying, and (ii) slow mixing
trajectories—trajectories starting in regions with opposite OGM values mix
much slower compared to, for example, those in the near-invariant regions
identified by the SDM (Figures 2B and 3B).

Transient slow mixing has consequences for the system’s response to
perturbations. A slow approach to equilibrium translates to increased
sensitivity to perturbations due to the slow decay of the Green function [2,
19, 66, 67]. We expect the OGM to be associated with perturbations that take
longer to enter the exponentially decaying regime, as suggested by the L63
system’s relaxation curve (Fig. 1 A) and the system’s non-equilibrium nature.
The OGM is optimized over time for transient behavior, so it may not be ideal
for optimizing the steady-state response to perturbations. Prior work has
addressed finding optimal perturbations in the linear response regime ([27,
68, 59]). Our future work will explore the implications of non-normality for
linear response in more detail.

_Acknowledgements._ —MSG is grateful to the Feinberg Graduate School for their
support through the Dean of Faculty Fellowship. This work has been supported
by the Office of Naval Research (ONR) Multidisciplinary University Research
Initiative (MURI) grant N00014-20-1-2023, and by the National Science
Foundation grant DMS-2407484. This work has been also supported by the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program (grant agreement No. 810370).

The authors declare that there are no competing interests.

_Appendix: Correlations approximation._ — In computing $M_{\tau}$ in Eq. (4),
there are two main parameters that affect the computations: the number of
boxes $N$, and the transition time $\tau$. While increasing $N$ improves the
approximation [22], in this appendix, we shall illustrate on the L63 system
how $\tau$ influences the approximation of correlation functions and power-
spectra using the decomposition formulas of [23, Sec. 2.3]. We furthermore
compare the performance when the underlying resonances are estimated from the
(3D) transition matrices $M_{\tau}$ (Eq. (4) with
$\bm{x}_{k}=(x_{k},y_{k},z_{k})$) versus those obtained from the reduced
transition matrix (Eq. (4) with $\bm{x}_{k}=(x_{k},z_{k})$), estimated in the
$(x,z)$-reduced state space. For the computations in 3D, the $(x,y,z)$-domain
is $[-20,20]\times[-30,30]\times[0,50]$ that is divided into $N=2^{15}$
equally sized, non-intersecting boxes $\\{B_{i}\\}_{i=1}^{N}$ of dimensions
$40/2^{5}\times 60/2^{5}\times 50/2^{5}$. Then, the full transition matrix
$M_{\tau}$ is computed for different values of $\tau$ as detailed below.

In Figure 4A we show the autocorrelation functions (ACFs) of the $z$ variable
of the L63 model in a number of ways. The black solid line is calculated as
the sample autocorrelation function of a $10^{5}$ time unit time series, and
the dotted lines are computed using Markov matrices over the three-dimensional
attractor, with values of $\tau$ indicated in the legend. In Panel B, dotted
lines are computed in the $(x,z)$-reduced state space. In Panels C and D, we
show the corresponding power spectral densities (PSDs). The results show that
the oscillatory component of the correlation function is well captured
regardless of the value of $\tau$; see the peaks in such panels. The readouts
also suggest that an increasing value of $\tau$ improves the correlation
function approximation, yet it compromises the approximation of the low-
frequency variability portrayed by the PSDs at the smallest $\omega$-values in
panels C and D. This surprising ability of the Markov chain method to capture
the frequency of oscillations was also reported for stochastic differential
equations in [23, Sec. 4].

Figure 4:  Panels A and C show the autocorrelation function and power spectral
density of the variable $z$ for the L63 model. The black lines are obtained
from the integration of the full model. The dotted lines correspond to the
correlation functions as computed from the spectrum of the (full) transition
matrix $M_{\tau}$ (Eq. (4) with $\bm{x}_{k}=(x_{k},y_{k},z_{k})$), for values
of $\tau$ as indicated in the legend. Panels B and D, show correlation
functions and PSDs as computed from the spectrum of the reduced transition
matrix (Eq. (4) with $\bm{x}_{k}=(x_{k},z_{k})$), estimated in the
$(x,z)$-reduced state space.

## References

  * Lasota and Mackey [1994] A. Lasota and M. C. Mackey,  _Chaos, fractals and noise_ (Springer, New York, 1994).
  * Ruelle [2009] D. Ruelle, Nonlinearity 22, 855 (2009).
  * Kubo [1957] R. Kubo, Journal of the Physical Society of Japan 12, 570 (1957).
  * Kubo [1966] R. Kubo, Reports on Progress in Physics 29, 255 (1966).
  * Young [2002] L.-S. Young, Journal of Statistical Physics 108, 733 (2002).
  * Ruelle [1998] D. Ruelle, Nonlinearity 11, 5 (1998).
  * Landau and Lifshitz [2013] L. D. Landau and E. M. Lifshitz,  _Statistical Physics: Volume 5_ , Vol. 5 (Elsevier, 2013).
  * Reddy and Henningson [1993] S. C. Reddy and D. S. Henningson, Journal of Fluid Mechanics 252, 209 (1993).
  * Farrell and Ioannou [1996] B. F. Farrell and P. J. Ioannou, Journal of Atmospheric Sciences 53, 2025 (1996).
  * Trefethen and Embree [2005] L. N. Trefethen and M. Embree,  _Spectra and Pseudospectra_ (Princeton University Press, 2005).
  * Schütte  _et al._ [1999] C. Schütte, A. Fischer, W. Huisinga, and P. Deuflhard, J. Stat. Phys. 151, 146 (1999).
  * Klus  _et al._ [2018] S. Klus, F. Nüske, P. Koltai, H. Wu, I. Kevrekidis, C. Schütte, and F. Noé, Journal of Nonlinear Science 28, 985 (2018).
  * Ulam [1960] S. M. Ulam, A Collection of Mathematical Problems (1960).
  * Junge and Koltai [2009] O. Junge and P. Koltai, SIAM Journal on Numerical Analysis 47, 3464 (2009).
  * Froyland  _et al._ [2021] G. Froyland, D. Giannakis, B. R. Lintner, M. Pike, and J. Slawinska, Nature Communications 12, 6570 (2021).
  * Meyn and Tweedie [2012] S. P. Meyn and R. L. Tweedie,  _Markov Chains and Stochastic Stability_ (Springer Science & Business Media, 2012).
  * Hairer [2010] M. Hairer, Lecture notes, University Warwick 18, 13 (2010).
  * Froyland and Dellnitz [2003] G. Froyland and M. Dellnitz, SIAM Journal on Scientific Computing 24, 1839 (2003).
  * Chekroun  _et al._ [2014] M. D. Chekroun, J. D. Neelin, D. Kondrashov, J. C. McWilliams, and M. Ghil, Proc. Natl. Acad. Sci. USA 111, 1684 (2014).
  * Froyland  _et al._ [2013] G. Froyland, O. Junge, and P. Koltai, SIAM J. Numerical Analysis 51, 223 (2013), arXiv:1101.4166 .
  * Froyland [1998] G. Froyland, Nonlinear Analysis 32, 831 (1998).
  * Fishman and Rahav [2002] S. Fishman and S. Rahav, in  _Dynamics of Dissipation_ (Springer, Berlin, 2002) pp. 165–192, arXiv:0204068 [nlin] .
  * Chekroun  _et al._ [2020] M. Chekroun, A. Tantet, H. Dijkstra, and J. D. Neelin, J. Stat. Phys. 179, 1366 (2020).
  * Froyland [2005] G. Froyland, Physica D: Nonlinear Phenomena 200, 205 (2005).
  * Froyland and Padberg [2009a] G. Froyland and K. Padberg, Physica D: Nonlinear Phenomena 238, 1507 (2009a).
  * Froyland  _et al._ [2016] G. Froyland, C. González-Tokman, and T. Watson, SIAM Review 58, 494 (2016).
  * Antown  _et al._ [2018] F. Antown, D. Dragičević, and G. Froyland, Journal of Statistical Physics 170, 1051 (2018).
  * Froyland  _et al._ [2020] G. Froyland, P. Koltai, and M. Stahn, SIAM Journal on Applied Dynamical Systems 19, 1659 (2020).
  * Froyland  _et al._ [2007] G. Froyland, K. Padberg, M. H. England, and A. M. Treguier, Physical Review Letters 98, 1 (2007).
  * Froyland [2013] G. Froyland, Physica D 250, 1 (2013).
  * Da Prato and Zabczyk [1996] G. Da Prato and J. Zabczyk,  _Ergodicity for infinite dimensional systems_ , Vol. 229 (Cambridge university press, 1996).
  * Risken [1989] H. Risken,  _The Fokker-Planck Equation_ , 2nd ed. (Springer, 1989).
  * Engel and Nagel [2000] K.-J. Engel and R. Nagel,  _One-Parameter Semigroups for Linear Evolution Equations_ (Springer-Verlag, New York, 2000).
  * Ruelle [1986] D. Ruelle, Physical Review Letters 56, 405 (1986).
  * Baladi [2000] V. Baladi,  _Positive Transfer Operators and Decay of Correlations_ , Vol. 16 (World scientific, 2000).
  * Hairer [2002] M. Hairer, Probability Theory and Related Fields 124, 345 (2002).
  * Hairer and Mattingly [2008] M. Hairer and J. C. Mattingly, The Annals of Probability 36, 2050 (2008).
  * MacTaggart and Stewart [2017] D. MacTaggart and P. Stewart, Sol. Phys. 292 (2017).
  * Penland and Magorian [1993] C. Penland and T. Magorian, Journal of Climate 6, 1067 (1993).
  * Penland and Sardeshmukh [1995] C. Penland and P. D. Sardeshmukh, Journal of climate 8, 1999 (1995).
  * Tziperman and Ioannou [2002] E. Tziperman and P. J. Ioannou, Journal of Physical Oceanography 32, 3427 (2002).
  * Zanna and Tziperman [2005] L. Zanna and E. Tziperman, Journal of physical oceanography 35, 1593 (2005).
  * Farrell and Ioannou [1993] B. F. Farrell and P. J. Ioannou, Journal of Atmospheric Sciences 50, 4044 (1993).
  * Bousquet  _et al._ [2015] A. Bousquet, M. D. Chekroun, Y. Hong, R. M. Temam, and J. Tribbia, Mathematics of Climate and Weather Forecasting 1 (2015).
  * Kato [1966] T. Kato,  _Perturbation Theory for Linear Operators_ (Springer-Verlag, 1966).
  * Pavliotis [2014] G. A. Pavliotis,  _Book_ , Vol. 60 (Springer, New York, 2014).
  * Metafune  _et al._ [2002] G. Metafune, D. Pallara, and E. Priola, Journal of Functional Analysis 196, 40 (2002).
  * Gaspard  _et al._ [1995] P. Gaspard, G. Nicolis, A. Provata, and S. Tasaki, Physical Review E 51, 74 (1995).
  * Tantet  _et al._ [2020a] A. Tantet, M. D. Chekroun, H. A. Dijkstra, and J. D. Neelin, Journal of Statistical Physics 10.1007/s10955-020-02526-y (2020a).
  * Crommelin and Vanden-Eijnden [2011] D. Crommelin and E. Vanden-Eijnden, Multiscale Modeling & Simulation 9, 1588 (2011).
  * Tantet  _et al._ [2018] A. Tantet, V. Lucarini, and H. A. Dijkstra, Journal of Statistical Physics 170, 584 (2018), arXiv:1705.08178 .
  * Tantet  _et al._ [2020b] A. Tantet, M. D. Chekroun, H. A. Dijkstra, and J. D. Neelin, J. Stat. Phys. 179, 1449 (2020b).
  * Lorenz [1963] E. N. Lorenz, Journal of the Atmospheric Sciences 20, 130 (1963).
  * Ghil and Childress [2012] M. Ghil and S. Childress,  _Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dynamo Theory, and Climate Dynamics_ , Vol. 60 (Springer Science & Business Media, 2012).
  * Sparrow [1982] C. Sparrow,  _The Lorenz Equations_ (Springer, 1982).
  * Anishchenko  _et al._ [2005] V. S. Anishchenko, E. V. Tat’yana, G. A. Okrokvertskhov, and G. I. Strelkova, Physics-Uspekhi 48, 151 (2005).
  * Chekroun  _et al._ [2011] M. D. Chekroun, E. Simonnet, and M. Ghil, Physica D: Nonlinear Phenomena 240, 1685 (2011).
  * Reick [2002] C. H. Reick, Physical Review E 66, 36103 (2002).
  * Lucarini [2016] V. Lucarini, Journal of Statistical Physics 162, 312 (2016).
  * Santos-Gutiérrez and Lucarini [2020] M. Santos-Gutiérrez and V. Lucarini, Journal of Statistical Physics 179, 1572 (2020).
  * Tucker [2002] W. Tucker, Foundations of Computational Mathematics 2, 53 (2002).
  * Kondrashov  _et al._ [2015] D. Kondrashov, M. D. Chekroun, and M. Ghil, Physica D: Nonlinear Phenomena 297, 33 (2015).
  * Lucarini and Chekroun [2023] V. Lucarini and M. D. Chekroun, Nature Reviews Physics 5, 744 (2023), doi:10.1038/s42254-023-00650-8.
  * Chekroun  _et al._ [2017] M. D. Chekroun, H. Liu, and J. C. McWilliams, Computers & Fluids 151, 3 (2017).
  * Froyland and Padberg [2009b] G. Froyland and K. Padberg, Physica D 238, 1507 (2009b).
  * Lucarini [2008] V. Lucarini, Journal of Statistical Physics 131, 543 (2008).
  * Gutiérrez and Lucarini [2022] M. S. Gutiérrez and V. Lucarini, Journal of Physics A: Mathematical and Theoretical 55, 425002 (2022).
  * Schweitzer [1968] P. J. Schweitzer, Journal of Applied Probability 5, 401 (1968).