# Rapid and facile reconstruction of time-resolved fluorescence data with
exponentially modified Gaussians

Darien J. Morrow<EMAIL_ADDRESS>Center for Nanoscale Materials,
Argonne National Laboratory, Lemont, Illinois 60439, United States  Xuedan Ma
<EMAIL_ADDRESS>Center for Nanoscale Materials, Argonne National Laboratory,
Lemont, Illinois 60439, United States  Consortium for Advanced Science and
Engineering, University of Chicago, Chicago, Illinois 60637, United States
Northwestern-Argonne Institute of Science and Engineering, 2205 Tech Drive,
Evanston, IL 60208, USA

###### Abstract

Analyte response is convoluted with instrument response in time resolved
fluorescence data. Decoding the desired analyte information from the
measurement usually requires iterative numerical convolutions. Here in, we
show that time resolved data can be completely, analytically reconstructed
without numerical convolutions. Our strategy relies on a summation of
exponentially modified Gaussians which encode all convolutions within easily
evaluated complementary error functions. Compared to a numerical convolution
strategy implemented with Python, this new method is computationally cheaper
and scales less steeply with the number of temporal points in the experimental
dataset.

##  I Introduction

Time resolved fluorescence is a ubiquitous measurement used in fields as
diverse as materials physics,[1] physical chemistry,[2] structural biology,[3,
4] robotic surgery,[5] and art conservation.[6] In these uses, the metric of
interest is the set of fluorescence lifetimes, $\\{\tau_{i}\\}$, which is
defined by the physical identity and structure of the sample. In the Platonic
ideal of the experiment, a sample is excited with an infinitely fast pulse of
light, the sample then undergoes spontaneous emission with first order
kinetics, so the emitted intensity of light decays exponentially with time.
The emitted light is finally measured with an infinitely fast detector. In
actuality, excitation pulses have finite width and detectors have finite
response times (Figure 1). These non-idealities must be accounted for in data
analysis routines to extract useful and representative fluorescence lifetimes.

If $S$ is the ideal signal, the measured response, $M$, is the convolution of
the ideal signal with some instrument response function, $R$, which accounts
for the non-idealities of reality,

$\displaystyle M$ $\displaystyle=S*R.$ (1)

Deconvolution of $S$ from $M$ can be accomplished if $R$ is exactly known. If
$R$ is not known with certainty, then deconvolution is a non-covex (usually
unfeasible) optimization problem.[7] Instead most commercial software packages
and individual researchers rely on a reconvolution strategy in which a decay
model, $S_{\text{model}}$, is convolved with a measured (or supposed)
instrument response function, $R_{\text{measured}}$,

$\displaystyle M_{\text{model}}$
$\displaystyle=S_{\text{model}}*R_{\text{measured}}.$ (2)

The resultant signal, $M_{\text{model}}$, is compared to $M$. The parameters
which define $S_{\text{model}}$ are then iteratively adjusted to minimize the
difference between $M_{\text{model}}$ and $M$ in order to extract a
representation of $S$. In the case of methods like Fluorescence Lifetime
Imaging (FLIM), the extraction of $S$ must then be accomplished for thousands
of time traces which each represents a pixel of a spatial map.[8]

Figure 1:  Actual time-resolved flourescence data, $M$, versus the ideal
exponential decay predicted from first order kinetics, $S$.
$\tau=1/\lambda=1.5\text{ ns}$.

The reconvolution strategy is usually robust and generates an accurate view of
the dynamics encoded in fluorescence data. In practice, it can be difficult to
correctly zero-pad and interpolate $S_{\text{model}}$ and $R$. Moreover, if
one does not use an easily scriptable analysis software, then one must iterate
the parameters defining $S_{\text{model}}$ manually, which can be time
consuming. Closed source commercial software packages (e.g. SPCImage NG from
Becker & Hickl and EasyTau 2 from Picoquant) and open source packages (e.g.
DecayFit and FLIMfit,[9]) exist to help solve the problems associated with
fitting fluorescence data to extract lifetimes via numerical convolutions.
Without using these packages, in many cases, researchers choose to assume the
mid to late time dynamics encoded in their datasets are representative of $S$
and therefore just fit $S_{\text{model}}$ to
$M\left(t>t_{\text{cutoff}}\right)$, which is a good assumption in some
limits. It is likely that this choice is oftentimes made not for suitability
of assumption, but instead because it can be hard to correctly implement
numerical convolutions.

Herein we introduce an analytic solution to this numerical convolution
problem. We derive an integral equation which accounts for any number of
exponential decay processes in both the sample response and instrument
response. All numerically difficult convolutions are accounted for
analytically. Our result is constructed so that analysis of time-resolved
fluorescence data exclusively uses the robust algorithms which are packaged in
most software packages to calculate the values of exponential and error
functions. We demonstrate the utility of our result for reconstructing the
wavelength dependent response of a commercial single-photon avalanche diode
(SPAD) and also for retrieving the decay lifetimes of a model system. Finally
we show that our method is computationally cheaper than more established
methods.

##  II Derivation of primary result

In this section we derive our primary result (Equation 27).

###  II.1 Exponentially modified Gaussian: h(x)

We start by assuming that fluorescence decay processes can be described using
monoexponential decays

$\displaystyle f(x;\lambda)$ $\displaystyle\equiv\lambda\exp{\left[-\lambda
x\right]}\Theta\left[x\right],$ (3) $\displaystyle\Theta\left[x\right]$
$\displaystyle\equiv\begin{cases}1,x>0\\\ 0,x\leq 0\end{cases}.$ (4)

in which $\Theta$ is the Heaviside step function, which accounts for the
sample not emitting until after it is excited. The lifetime of a component is
given by $\tau_{i}=\frac{1}{\lambda_{i}}$. The ideal instrument response
function (IRF) is considered to be a Gaussian

$\displaystyle g(x;\mu,\sigma)$
$\displaystyle\equiv\frac{1}{\sigma\sqrt{2\pi}}\exp{\left[-\left(\frac{x-\mu}{\sqrt{2}\sigma}\right)^{2}\right]},$
(5)

with $\sigma$ characterizing the temporal width and $\mu$ defining the “time
zero” of the instrument. Later we will account for non-Gaussian IRFs. Note
that, for clarity, $f(x)$ and $g(x)$ are both written as normalized functions
so that $\int_{-\infty}^{\infty}f(x)\textrm{d}x=1$—this choice does not affect
the form of our final result.

We now define a composite function which will be the mainstay of our
derivation,

$\displaystyle h(x;\mu,\sigma,\lambda)$ $\displaystyle\equiv
f(x;\lambda)*g(x;\mu,\sigma),$ (6) $\displaystyle(f*g)(t)$
$\displaystyle\equiv\int_{-\infty}^{\infty}f(\tau)g(t-\tau)\textrm{d}\tau.$
(7)

Here, $*$ is the convolution operator which, importantly, is commutative,
associative, and distributive. In Appendix A we show that $h$ is a function
known as the _exponentially modified Gaussian_ (EMG) with form,

$\displaystyle\begin{split}h(x;\mu,\sigma,\lambda)=\frac{\lambda}{2}\exp{\left[\frac{\lambda}{2}\left(2\mu+\lambda\sigma^{2}-2x\right)\right]}\\\
\times\operatorname{erfc}{\left[\frac{\mu+\lambda\sigma^{2}-x}{\sqrt{2}\sigma}\right]},\end{split}$
(8)

where

$\displaystyle\operatorname{erfc}{(t)}\equiv
1-\operatorname{erf}{(t)}=\frac{2}{\sqrt{\pi}}\int_{t}^{\infty}\exp{\left[-\tau^{2}\right]}\textrm{d}\tau.$
(9)

In Equation 8, the convolution operation present in Equation 6 is rewritten as
a single _complementary error_ function. The EMG is commonly used in many
fields as a phenomenological quantification tool: in the chromatography and
mass spectrometry fields it is used to quantify tailed lineshapes,[10, 11, 12]
in psychophysiology it is used to quantify response times,[13, 14] and in
material science it has seen limited use fitting fluorescence spectra and
diffusion profiles from inhomogeneous semiconductors.[15, 16, 17, 18, 19]

Figure 2a graphs $f$, $g$, and $h$. At early times $h$ is similar to the
Gaussian IRF, $g$, but at late times, $h$ follows the exponential decay, $f$.
Figure 2 validates the intuition that if the dynamics of interest are much
longer than the width of the IRF, then one can extract the lifetimes of
interest by fitting the long-time tail of the measured response. Figure 2b
shows how larger widths of the IRF ($\sigma$) pushes the peak maximum of $h$
to later times.

Figure 2: Graphical form of the exponentially modified Gaussian. (a) Graph of
Equation 3, Equation 5, and Equation 6 for $\lambda=0.4$, $\sigma=1$ and
$\mu=0$. A scaled form of Equation 9 is also shown to highlight its sigmoid-
like behavior. (b) Equation 6 with $\lambda=0.5$ and
$\sigma\in\\{0.1,0.5,1\\}$ showing how larger values of $\sigma$ decreases the
angle of $h$ at negative $x$ and pushes the maximum to positive $x$.

###  II.2 Extension of h(x) with multiple decay functions

For multiple decay pathways we write a distribution of exponential decays,

$\displaystyle F^{(i)}(x)$
$\displaystyle\equiv\sum_{i}^{m}a_{i}\lambda_{i}\exp{\left[-\lambda_{i}x\right]}\Theta\left[x\right],$
(10) $\displaystyle=\sum_{i}^{m}a_{i}f_{i}(x),$ (11)

with $a_{i}$ being a weighting factor for each exponential decay with
$\alpha_{i}\equiv a_{i}\lambda_{i}$ being the unnormalized weight. Convolving
this distribution of decays with a Gaussian yields

$\displaystyle H^{(i)}(x)$ $\displaystyle\equiv g(x)*F^{(i)}(x),$ (12)
$\displaystyle=\sum_{i}^{m}g(x)*\left[a_{i}f_{i}(x)\right],$ (13)
$\displaystyle=\sum_{i}^{m}a_{i}h_{i}(x),$ (14)

in which we noted that convolution is a linear operation and therefore
distributive. Here we emphasize that $h_{i}(x)=h(x;\mu,\sigma,\lambda_{i})$
which means that there are single $\mu$ and $\sigma$ but multiple
$\lambda_{i}$. If a time-resolved fluorescence experiment has a Gaussian IRF,
then $H^{(i)}(x)$ describes the measured response of a sample with a set,
$\\{\lambda_{i}\\}$, of decay pathways.

###  II.3 Detectors with non-Gaussian, tailed IRFs

Oftentimes SPAD detectors do _not_ have Gaussian IRFs, but instead have a long
tail towards positive times. Phenomenologically these non-ideal IRFs may be
described with a distribution of weights and rates. In this case we interpret
$H^{(i)}(x)$ as the IRF. We then account for analyte response by convolving
$H^{(i)}(x)$ (our IRF) with a second $F^{(j)}(x)$ (the true sample decay
pathways) to form a representative measured signal, $\mathcal{M}$,

$\displaystyle\mathcal{M}(x)$ $\displaystyle=H^{(i)}(x)*F^{(j)}(x),$ (15)

which may be rewritten as

$\displaystyle\mathcal{M}(x)$ $\displaystyle=g(x)*F^{(i)}(x)*F^{(j)}(x),$ (16)
$\displaystyle=g(x)*\left\\{\sum_{i}^{m}f_{i}(x)\right\\}*\left\\{\sum_{j}^{n}f_{j}(x)\right\\},$
(17) $\displaystyle=g(x)*\left\\{\sum_{i,j}^{m,n}f_{i}(x)*f_{j}(x)\right\\},$
(18) $\displaystyle=g(x)*\left\\{\sum_{i,j}^{m,n}y_{i,j}(x)\right\\}.$ (19)

To solve the double convolution present in Equation 19 we first consider the
argument of the summand inside the braces

$\displaystyle y_{i,j}(t)$ $\displaystyle=f_{i}(t)*f_{j}(t),$ (20)
$\displaystyle=a_{i}\lambda_{i}a_{j}\lambda_{j}\int_{-\infty}^{\infty}e^{\left[-\lambda_{i}\tau\right]}\Theta\left[\tau\right]e^{\left[-\lambda_{j}(t-\tau)\right]}\Theta\left[t-\tau\right]\textrm{d}\tau.$
(21)

This integrand is nonzero only when $\tau>0$ and $t-\tau>0$ so

$\displaystyle y_{i,j}(t)$
$\displaystyle=a_{i}\lambda_{i}a_{j}\lambda_{j}\exp{\left[-\lambda_{j}\tau\right]}\Theta\left[t\right]\int_{0}^{t}\exp{\left[\tau(\lambda_{j}-\lambda_{i})\right]}\textrm{d}\tau,$
(22)
$\displaystyle=\left(\frac{a_{j}\lambda_{j}}{\lambda_{j}-\lambda_{i}}\right)f_{i}(t)+\left(\frac{a_{i}\lambda_{i}}{\lambda_{i}-\lambda_{j}}\right)f_{j}(t),$
(23)

in which we noted the standard integral

$\displaystyle\int_{0}^{t}\exp{\left[\tau(\lambda_{j}-\lambda_{i})\right]}\textrm{d}\tau=\frac{1-\exp{\left[(\lambda_{j}-\lambda_{i})t\right]}}{\lambda_{j}-\lambda_{i}}\text{,
for }\lambda_{j}\neq\lambda_{i}.$ (24)

Equation 23 is only valid when the sample’s decay rates are not the same as
the instrument response decay rates. Practically, this requirement is
satisfied when the sample has slower dynamics than the detector.

Substitution of $y_{i,j}$ from Equation 23 into Equation 19 yields our
measured response

$\displaystyle\mathcal{M}(x)$
$\displaystyle=g(x)*\sum_{i,j}^{m,n}\left\\{\left(\frac{a_{j}\lambda_{j}f_{i}(x)}{\lambda_{j}-\lambda_{i}}\right)+\left(\frac{a_{i}\lambda_{i}f_{j}(x)}{\lambda_{i}-\lambda_{j}}\right)\right\\},$
(25)
$\displaystyle=\sum_{i,j}^{m,n}\left\\{\left(\frac{a_{j}\lambda_{j}}{\lambda_{j}-\lambda_{i}}\right)h_{i}(x)+\left(\frac{a_{i}\lambda_{i}}{\lambda_{i}-\lambda_{j}}\right)h_{j}(x)\right\\},$
(26)
$\displaystyle=\sum_{i}^{m}h_{i}(x)\sum_{j}^{n}\frac{a_{j}\lambda_{j}}{\lambda_{j}-\lambda_{i}}+\sum_{j}^{n}h_{j}(x)\sum_{i}^{m}\frac{a_{i}\lambda_{i}}{\lambda_{i}-\lambda_{j}}.$
(27)

Equation 27 is our primary analytic result. It encodes the dynamics of an
experiment whose IRF has an arbitrarily large number of exponential tails, and
whose sample response also has an arbitrarily large number of monoexponential
decays. Importantly, there are no explicit convolutions present in Equation
27. Instead, Equation 27 requires evaluation of complementary error functions;
software libraries as diverse as Microsoft Excel and the Scientific Python
stack provide this capability using efficient, accurate algorithms.[20, 21]
Appendix B shows two specific cases of Equation 27 with the nested summations
worked out.

##  III Comparison to experimental results

We demonstrate the utility and suitability of our strategy of exponentially
modified Gaussians for time-resolved florescence in two case studies. The
figure of merit (model error, $\mathcal{E}$) to minimize while fitting is the
weighted square residuals,[22, 7]

$\displaystyle\mathcal{E}=\sum_{i}^{n}w_{i}\left(\mathcal{F}_{\text{data}}[t_{i}]-\mathcal{F}_{\text{model}}[t_{i}]\right)^{2},$
(28)

where $\mathcal{F}[t_{i}]$ is the measured fluorescence counts or model
prediction at the time point $t_{i}$. The weighting factor, $w_{i}$, is
related to the variance, $\sigma_{i}^{2}$, at each data point

$\displaystyle w_{i}$
$\displaystyle=\frac{1}{\sigma_{i}^{2}}\left[\frac{1}{n}\sum_{i}^{n}\frac{1}{\sigma_{i}^{2}}\right]^{-1},$
(29)
$\displaystyle=\frac{1}{\mathcal{F}_{\text{data}}[t_{i}]}\left[\frac{1}{n}\sum_{i}^{n}\frac{1}{\mathcal{F}_{\text{data}}[t_{i}]}\right]^{-1},$
(30)

in which we have used the fact that photon-counting data without systematic
error is well described by a Poisson distribution so the variance and count
number, $\mathcal{F}_{\text{data}}$, are equal for a specific time interval.
To visually access the goodness of fit, we can plot the weighted residuals,

$\displaystyle\text{weighted
residuals}=\sqrt{w_{i}}\left(\mathcal{F}_{\text{data}}[t_{i}]-\mathcal{F}_{\text{model}}[t_{i}]\right).$
(31)

###  III.1 Fitting wavelength dependent IRF

In the first study, we characterize the instrument response of a commercial
SPAD (PDM from Micro Photon Devices) and counting electronics (PicoQuant
HydraHarp 400) by excitation with what is effectively a Dirac delta impulse.
We couple the attenuated output of an ultrafast tuneable Ti:sapphire laser
(Coherent Chameleon Discovery with Harmonics Generator, $\sim 100$ fs pulse
width, 80 MHz repetition rate) into the SPAD. By tuning the ultrafast laser
from 500 to 1075 nm (Figure 3c) and keeping SPAD count rates constant, we
recover the wavelength dependent response of the SPAD and associated
electronics.

We fit the IRFs with $H^{(i)}(x)$ (Equation 14) and find that three
exponential decays are required for recapitulation of the measurement. Figure
3b shows the weighted residuals and their autocorrelation. At delay times
before 0.5 ns, the residuals oscillate about zero. This indicates that the
rising edge of the IRF is not perfectly described by a Gaussian and so
Equation 14 does not exactly capture the instrument response.

The decay profiles shown in Figure 3a indicates that as one excites a SPAD
with higher wavelength light, the response becomes longer.[23] The intensity
weighted average decay time (Figure 3a inset),[24]

$\displaystyle\tau_{\text{avg}}=\frac{\sum_{i}\frac{\alpha_{i}}{\lambda_{i}^{2}}}{\sum_{i}\frac{\alpha_{i}}{\lambda_{i}}},$
(32)

codifies this observation of longer instrument response times for longer
wavelength excitation.

Figure 3: Wavelength dependent SPAD IRF. (a) Temporal response from 500 to
1075 nm with thin, solid lines being the model, $H^{(i)}(x)$. Pulse profile of
excitation laser is shown as gray (note how it is effectively a delta function
on the nanosecond scale). Inset shows the intensity weighted average (Equation
32) of the three decay rates fit to the IRFs. (b) Weighted residuals of fits
and their zoomed-in autocorrelation. (c) Color coded spectra of laser used to
excite SPAD.

###  III.2 Fitting wavelength dependent POPOP fluorescence

In the second study, we characterize the response of powdered POPOP;
1,4-Bis(5-phenyl-2-oxazolyl)benzene; whose chemical structure is shown in
Figure 4a. POPOP is a scintillator with an absorption onset in the UV-A and
bright fluorescence across the visible. With an objective we focus the output
of a 375 nm pulsed picosecond diode laser (Picoquant, 20 MHz repetition rate)
onto powdered POPOP. Fluorescence is collected with the same objective,
filtered with a 425, 525, or 610 nm bandpass filter, and focused into a photon
counting module (PerkinElmer, PicoQuant HydraHarp 400 electronics). Figure 4c
shows a florescence spectrum of POPOP and the throughput of the three bandpass
filters used.

Figure 4: Fluorescence of POPOP excited with a 375 nm pulsed didode laser. (a)
Measured and modeled time-resolved florescence. The measured IRF of our SPAD
at 375 nm is shown in gray. We recover the following decay rates: for 425 nm,
0.09, 0.14, and 0.16 ns-1; for 525 nm, 0.34, 0.39, and 0.61 ns-1; and finally
for 610 nm, 0.59, 0.90, and 1.34 ns-1. A fit of the 610 nm dataset with
Grinvald and Steinberg [22]’s iterative method is shown in blue. (b) Weighted
residuals of fits and their zoomed-in autocorrelation. (c) Fluorescence
intensity of POPOP measured through bandpass filters color-coded to traces in
(a) along with emission spectra of POPOP which has not been for the wavelength
dependent efficiency of our camera and collection optics.

We first measure the response function of the detector at 375 nm (Figure 4a).
This IRF is fit in the same way as the IRFs in Figure 3a. Next, we measure the
florescence of POPOP over three wavelength bands. Finally, we use our primary
result, Equation 27, to fit each time trace (tri-exponential) while holding
the IRF parameters constant. We recover average lifetimes of 1.1, 2.2, and 3.6
ns for the 425, 525, and 610 nm bandpass filters, respectively. The model fits
all aspects of the time-resolved fluorescence exceedingly well. The weighted
residuals (Figure 3b) show slight systematic deviation at the rising edge of
the flouresence near $t=0$.

To compare our model to other methods, we fit the longest wavelength data
using the iterative method of Grinvald and Steinberg [22] using the measured
IRF. The result of this fit is plotted in blue in Figure 3a. There is striking
agreement between our model and Grinvald’s; however, their method better
captures the rising edge of the instrument response and does not show a
systemic deviation in the early time weighted residuals.

##  IV Comparison of computation time to other methods

In this Section we show how the computational evaluation time of our EMG
formalism compares to other established methods. Specifically, we look at a
quadruple exponential decay convoluted with a Gaussian IRF and compare the
evaluation times of (1) a calculation using Equation 14, (2) step-wise
(iterative) construction of an IRF convoluted decay following the procedure of
Grinvald and Steinberg [22], and (3) numerical convolution of Equation 11 with
Equation 5. We show the code used to compare these strategies in Algorithms 1
and 2 in the Appendices. Base functions are defined in Algorithm 1 using the
Scientific Python ecosystem.[25, 26, 27, 28] Next, these base functions are
used in Algorithm 2 to construct singular functions which each take the _same_
input parameters and output a single convoluted decay profile. Note, in these
constructions it is imperative that the numerical convolution model and
Grinvald and Steinberg [22]’s iterative model each be evaluated with an evenly
spaced independent (x, delay time) array because the numerical convolution
relies on a fast-Fourier-transform and the iterative model uses a trapezoidal
(step-sized dependent) sequential integration technique. Our EMG model does
not have this restriction, and can be used on datasets with e.g.
logarithmically spaced time points.

Figure 5a shows the decay trace from each method. Note, when the traces are
not shifted from one-another, they are effectively identical. The edge affects
inherent in numerical convolutions have not been clipped out in order to
highlight their existence—our method and Grinvald and Steinberg [22]’s method
do not have these edge affects.

Figure 5: Comparison of Equation 14, numerical convolution of Equation 11 with
Equation 5, and the stepwise method of Grinvald and Steinberg [22]. (a) shows
that all functions are identical other than edge affects caused by the finite
width numerical convolution. (b) shows the average computational evaluation
time (Windows 10 running on an Intel Quad Core i7-6700 3.4 GHz processor) as
the number of temporal points is increased. The vertical purple bar indicates
the number of points (65536 histogram bins) recorded by our HydraHarp 400
time-correlated single photon counting electronics.

We use python’s built in timeit module to measure the _average_ execution time
of each function as the number of time points in the evaluation are
logarithmically increased. We also logarithmically _decrease_ the number of
function evaluations from 30000 to 30 as the number of points are increased.
Figure 5b shows the results of this study. In all cases, the EMG method is
faster than the other two methods. On average, the EMG method is $30\times$
faster than Grinvald and Steinberg [22]’s method. This increase in speed
likely originates from the fact that our computer program can calculate the
exponentially modified Gaussian array in one contiguous batch using an
optimized low-level function while the iterative method requires repeated
array calls and assignments.[26, 25]

The EMG method also scales less steeply with number of points ($x^{1.3}$) than
the numerical convolution method ($x^{2}$) which will be an important aspect
when handling large datasets. For datasets with the same number of bins,
65536, as our commercial counting electronics, the EMG method is two orders-
of-magnitude faster than direct numerical convolution.

##  V Discussion and Conclusions

In this work, we derived expressions for fitting the complete transient of
time-resolved fluorescence measurements without requiring the use of numerical
convolutions. We showed that our expressions well reproduce the tailed
instrument response of a SPAD. We then showed that the SPAD IRF could be held
constant while a second expression was fit to sample fluorescence. Our results
suggest a clear workflow:

  1. 1.

Measure instrument response and fit to

$H^{(i)}(x;\mu,\sigma,\\{a_{i}\\},\\{\lambda_{i}\\})=\sum_{i}^{m}a_{i}h_{i}(x)$

using the minimum number of exponential terms to recapitulate the measured
response.

  2. 2.

Measure sample response and fit to

$\mathcal{M}(x;\\{a_{j}\\},\\{\lambda_{j}\\})=\sum_{i,j}^{m,n}\left(\frac{a_{j}\lambda_{j}h_{i}(t)}{\lambda_{j}-\lambda_{i}}+\frac{a_{i}\lambda_{i}h_{j}(t)}{\lambda_{i}-\lambda_{j}}\right)$

using the minimum number of exponential terms to recapitulate the measured
response while holding the IRF parameters constant.

While this workflow seems to be robust, there are specific cases where this
method should not be applied. Firstly, if the material has dynamics which last
longer than the period between excitation pulses then there will be some
steady-state population which is not accounted for in our model. Secondly, if
the material has kinetics which are not collectively well described by first
order rate laws, then the assumption of exponential decay is poor. In both of
these cases, a more holistic treatment of the fluorescence is needed, indeed
there are software packages which can numerically solve the set of pulsed rate
equation necessary to correctly describe the system.[29, 30] Finally, some
SPADs are known to exhibit non-ideal behaviors like “afterpulsing” which can
lead to IRFs which have a train of humps on the tail (note that in Figure 3a
slight afterpulsing is present around 0.2 ns).[31, 32] If a detector’s IRF is
not well described by a collection of exponentially modified Gaussians which
share the same width and center, then our treatment is not adequate. Instead,
the linearity of convolutions allows one to have multiple weighted, and
shifted ($\\{\mu_{k}\\}$) versions of the IRF which results in a third
summation over $\mu_{k}$. However, this level of nested summations is likely
excessive and the researcher should instead accomplish convolutions with their
measured, IRF.

The method presented herein offers a single functional form for the entire
rise and decay of time resolved fluorescence data. This closed-form may be
useful to researchers who do not have access to bespoke software packages for
fitting time-resolved fluorescence data. Perhaps more importantly, because
this method has a faster computation time compared to other fitting methods,
analytical techniques like fluorescence lifetime imaging microscopy which rely
on fitting thousands of decay traces to build up an image will be able to use
exponentially modified Gaussians to increase the processing rate of their
large streams of data.

###### Acknowledgements.

This work was performed at the Center for Nanoscale Materials, a U.S.
Department of Energy Office of Science User Facility, and supported by the
U.S. Department of Energy, Office of Science, under Contract No. DE-
AC02-06CH11357.

##  Appendix A Showing $h(x)$ to be an EMG

To rewrite $h$ in a form which does not require one to perform numerical
convolutions, we must cast the convolution into a complementary error
function. We first write out the convolution integral in full

$\displaystyle h(x;\mu,\sigma,\lambda)$ $\displaystyle\equiv
f(x;\lambda)*g(x;\mu,\sigma),$ (33)
$\displaystyle=\left\\{\lambda\exp{\left[-\lambda
x\right]}\Theta\left[x\right]\right\\}*\left\\{\frac{1}{\sigma\sqrt{2\pi}}\exp{\left[-\left(\frac{x-\mu}{\sqrt{2}\sigma}\right)^{2}\right]}\right\\},$
(34)
$\displaystyle=\frac{\lambda}{\sigma\sqrt{2\pi}}\int_{-\infty}^{\infty}\exp{\left[-\lambda
t\right]}\Theta\left[t\right]\exp{\left[-\left(\frac{x-t-\mu}{\sqrt{2}\sigma}\right)^{2}\right]}\textrm{d}t.$
(35)

Because of the step function, the integrand is non-zero only when $t>0$ so we
can rewrite the bounds of integration

$\displaystyle h(x;\mu,\sigma,\lambda)$
$\displaystyle=\frac{\lambda}{\sigma\sqrt{2\pi}}\int_{0}^{\infty}\exp{\left[-\lambda
t\right]}\exp{\left[-\left(\frac{x-t-\mu}{\sqrt{2}\sigma}\right)^{2}\right]}\textrm{d}t.$
(36)

We will perform this integration by substitution. We define a new variable,
$u$,

$\displaystyle u$
$\displaystyle\equiv\frac{\mu+\lambda\sigma^{2}-x+t}{\sqrt{2}\sigma}$ (37)
$\displaystyle\Longrightarrow\textrm{d}u$
$\displaystyle=\frac{\textrm{d}t}{\sqrt{2}\sigma}$ (38)
$\displaystyle\Longrightarrow t$ $\displaystyle=\sqrt{2}\sigma
u-\mu-\lambda\sigma^{2}+x$ (39)

Proper substitution of $u$ into Equation 36 yields

$\displaystyle h(x;\mu,\sigma,\lambda)$
$\displaystyle=\frac{\sqrt{2}\sigma\lambda}{\sigma\sqrt{2\pi}}\int_{0}^{\infty}\exp{\left[-\lambda\left(\sqrt{2}\sigma
u-\mu-\lambda\sigma^{2}+x\right)\right]}\exp{\left[-\left(\frac{x-\left(\sqrt{2}\sigma
u-\mu-\lambda\sigma^{2}+x\right)-\mu}{\sqrt{2}\sigma}\right)^{2}\right]}\textrm{d}t,$
(40)
$\displaystyle=\frac{\lambda}{\sqrt{\pi}}\int_{\frac{\mu+\lambda\sigma^{2}-x}{\sqrt{2}\sigma}}^{\infty}\exp{\left[\lambda\mu-\lambda
x+\frac{\lambda^{2}\sigma^{2}}{2}\right]}\exp{\left[-u^{2}\right]}\textrm{d}u.$
(41) $\displaystyle=\frac{\lambda}{\sqrt{\pi}}\exp{\left[\lambda\mu-\lambda
x+\frac{\lambda^{2}\sigma^{2}}{2}\right]}\int_{\frac{\mu+\lambda\sigma^{2}-x}{\sqrt{2}\sigma}}^{\infty}\exp{\left[-u^{2}\right]}\textrm{d}u.$
(42)

The integral is now in the form of the complementary error function,
$\int_{t}^{\infty}\exp{\left[-u^{2}\right]}\textrm{d}u=\frac{\sqrt{\pi}}{2}\operatorname{erfc}{(t)}$,

$\displaystyle h(x;\mu,\sigma,\lambda)$
$\displaystyle=\frac{\lambda}{\sqrt{\pi}}\exp{\left[\lambda\mu-\lambda
x+\frac{\lambda^{2}\sigma^{2}}{2}\right]}\frac{\sqrt{\pi}}{2}\operatorname{erfc}{(t)},t=\frac{\mu+\lambda\sigma^{2}-x}{\sqrt{2}\sigma}$
(43)
$\displaystyle=\frac{\lambda}{2}\exp{\left[\frac{\lambda}{2}\left(2\mu+\lambda\sigma^{2}-2x\right)\right]}\operatorname{erfc}{\left[\frac{\mu+\lambda\sigma^{2}-x}{\sqrt{2}\sigma}\right]}.$
(44)

This is the desired result, an exponentially modified Gaussian.

##  Appendix B Common cases of Equation 27

In this Appendix we explicitly write out specific cases of Equation 27:

  1. 1.

An IRF with one exponential decay $i\in\\{1\\}$ and two exponential decays in
material response $j\in\\{\alpha,\beta\\}$

  2. 2.

An IRF with two exponential decays $i\in\\{1,2\\}$ and two exponential decays
in material response $j\in\\{\alpha,\beta\\}$

###  B.1 One exponential decay in IRF and two decays in material response

Consider the case of an IRF which has one exponential decay, $i\in\\{1\\}$ and
a material response with two exponential decays, $j\in\\{\alpha,\beta\\}$.

$\displaystyle\mathcal{M}_{\\{1\\},\\{\alpha,\beta\\}}(x)={}$ $\displaystyle
h_{1}(x)\left(\frac{a_{\alpha}\lambda_{\alpha}}{\lambda_{\alpha}-\lambda_{1}}+\frac{a_{\beta}\lambda_{\beta}}{\lambda_{\beta}-\lambda_{1}}\right)+a_{1}\lambda_{1}\left\\{\frac{h_{\alpha}(x)}{\lambda_{1}-\lambda_{\alpha}}+\frac{h_{\beta}(x)}{\lambda_{1}-\lambda_{\beta}}\right\\}$
(45) $\displaystyle h_{1}(x)={}$
$\displaystyle\frac{\lambda_{1}}{2}\exp{\left[\frac{\lambda_{1}}{2}\left(2\mu+\lambda_{1}\sigma^{2}-2x\right)\right]}\operatorname{erfc}{\left[\frac{\mu+\lambda_{1}\sigma^{2}-x}{\sqrt{2}\sigma}\right]}$
(46) $\displaystyle h_{\alpha}(x)={}$
$\displaystyle\frac{\lambda_{\alpha}}{2}\exp{\left[\frac{\lambda_{\alpha}}{2}\left(2\mu+\lambda_{\alpha}\sigma^{2}-2x\right)\right]}\operatorname{erfc}{\left[\frac{\mu+\lambda_{\alpha}\sigma^{2}-x}{\sqrt{2}\sigma}\right]}$
(47) $\displaystyle h_{\beta}(x)={}$
$\displaystyle\frac{\lambda_{\beta}}{2}\exp{\left[\frac{\lambda_{\beta}}{2}\left(2\mu+\lambda_{\beta}\sigma^{2}-2x\right)\right]}\operatorname{erfc}{\left[\frac{\mu+\lambda_{\beta}\sigma^{2}-x}{\sqrt{2}\sigma}\right]}$
(48)

First one would fit the measured IRF using the parameters
$\\{\sigma,\mu,a_{1},\lambda_{1}\\}$. Then these parameters would be held
constant and the full response would be fit using the parameters
$\\{a_{\alpha},a_{\beta},\lambda_{\alpha},\lambda_{\beta}\\}$.

###  B.2 Two exponential decays in both IRF and material response

Consider the case of an IRF which has two exponential decays, $i\in\\{1,2\\}$
and a material response also with two exponential decays,
$j\in\\{\alpha,\beta\\}$.

$\displaystyle\begin{split}\mathcal{M}_{\\{1,2\\},\\{\alpha,\beta\\}}(x)={}&+h_{1}(x)\left(\frac{a_{\alpha}\lambda_{\alpha}}{\lambda_{\alpha}-\lambda_{1}}+\frac{a_{\beta}\lambda_{\beta}}{\lambda_{\beta}-\lambda_{1}}\right)+h_{2}(x)\left(\frac{a_{\alpha}\lambda_{\alpha}}{\lambda_{\alpha}-\lambda_{2}}+\frac{a_{\beta}\lambda_{\beta}}{\lambda_{\beta}-\lambda_{2}}\right)\\\
&+h_{\alpha}(x)\left(\frac{a_{1}\lambda_{1}}{\lambda_{1}-\lambda_{\alpha}}+\frac{a_{2}\lambda_{2}}{\lambda_{2}-\lambda_{\alpha}}\right)+h_{\beta}(x)\left(\frac{a_{1}\lambda_{1}}{\lambda_{1}-\lambda_{\beta}}+\frac{a_{2}\lambda_{2}}{\lambda_{2}-\lambda_{\beta}}\right)\end{split}$
(49) $\displaystyle h_{1}(x)={}$
$\displaystyle\frac{\lambda_{1}}{2}\exp{\left[\frac{\lambda_{1}}{2}\left(2\mu+\lambda_{1}\sigma^{2}-2x\right)\right]}\operatorname{erfc}{\left[\frac{\mu+\lambda_{1}\sigma^{2}-x}{\sqrt{2}\sigma}\right]}$
(50) $\displaystyle h_{2}(x)={}$
$\displaystyle\frac{\lambda_{2}}{2}\exp{\left[\frac{\lambda_{2}}{2}\left(2\mu+\lambda_{2}\sigma^{2}-2x\right)\right]}\operatorname{erfc}{\left[\frac{\mu+\lambda_{2}\sigma^{2}-x}{\sqrt{2}\sigma}\right]}$
(51) $\displaystyle h_{\alpha}(x)={}$
$\displaystyle\frac{\lambda_{\alpha}}{2}\exp{\left[\frac{\lambda_{\alpha}}{2}\left(2\mu+\lambda_{\alpha}\sigma^{2}-2x\right)\right]}\operatorname{erfc}{\left[\frac{\mu+\lambda_{\alpha}\sigma^{2}-x}{\sqrt{2}\sigma}\right]}$
(52) $\displaystyle h_{\beta}(x)={}$
$\displaystyle\frac{\lambda_{\beta}}{2}\exp{\left[\frac{\lambda_{\beta}}{2}\left(2\mu+\lambda_{\beta}\sigma^{2}-2x\right)\right]}\operatorname{erfc}{\left[\frac{\mu+\lambda_{\beta}\sigma^{2}-x}{\sqrt{2}\sigma}\right]}$
(53)

First one would fit the measured IRF using the parameters
$\\{\sigma,\mu,a_{1},a_{2},\lambda_{1},\lambda_{2}\\}$. Then these parameters
would be held constant and the full response would be fit using the parameters
$\\{a_{\alpha},a_{\beta},\lambda_{\alpha},\lambda_{\beta}\\}$.

##  Appendix C Function evaluation time

In this Appendix we show the code used to determine the computational
evaluation time of our exponentially modified Gaussian formalism versus other
established methods. Specifically, we look at a quadruple exponential decay
convoluted with a Gaussian IRF. We compare the evaluation times of:

  1. 1.

Analytical calculation using Equation 14.

  2. 2.

Numerical convolution of Equation 11 with Equation 5.

  3. 3.

Stepwise (iterative) construction of an IRF convoluted decay following the
procedure of Grinvald and Steinberg [22].

Algorithm 1 python 3.8 construction of exponential, Gaussian, and
exponentially modified Gaussian functions.

    import numpy as np
    import timeit
    from scipy import special
    Ψ
    def f(x, l, a=1):
        # exponential decay with heaviside
        out = a * l * np.exp(-l * x)
        out[x < 0] = 0
        return out
    Ψ
    def g(x, mu, s):
        # Gaussian
        arg = -1 * ((x - mu) / (np.sqrt(2) * s)) ** 2
        return np.exp(arg) / (s * np.sqrt(2 * np.pi))
    Ψ
    def h(x, mu, s, l, a=1):
        # exponentially modified Gaussian (EMG)
        erfc = special.erfc
        arg1 = l / 2 * (2 * mu + l * s**2 - 2 * x)
        arg2 = (mu + l * s**2 - x) / (np.sqrt(2) * s)
        return a * l / 2 * np.exp(arg1) * erfc(arg2)
    Ψ
    def F(x, li, ai):
        # summation of exponentials with Heaviside
        # assumes x, li, and ai are all 1D. li and ai have same size
        out = ai[None, :] * li[None, :] * np.exp(-li[None, :] * x[:, None])
        out = np.sum(out, axis=-1)
        out[x < 0] = 0
        return out
    Ψ
    def H(x, mu, s, li, ai):
        # summation of EMGs with only one Gaussian convolved.
        # assumes x, li, and ai are all 1D. li and ai have same size
        out = h(x[:, None], mu, s, li[None, :], ai[None, :])
        return np.sum(out, axis=-1)

Algorithm 2 python 3.8 construction of quad-exponential functions and
evaluation time assessment.

    def quadexp_unpack(p):
        # unpack parameter array
        mu, s = p[0], p[1]
        li = p[2:6]
        ai = p[6:10]
        off = p[10]
        return mu, s, li, ai, off
    Ψ
    def quadexp_EMG_model(p, x):
        # exponentially modified Gaussian
        mu, s, li, ai, off = quadexp_unpack(p)
        return H(x, mu, s, li, ai) + off
    Ψ
    def quadexp_convolution_model(p, x):
        # direct numerical convolution of Gaussian with exponentials
        mu, s, li, ai, off = quadexp_unpack(p)
        IRF = g(x - x.mean(), mu, s)
        IRF /= IRF.sum()  # area normalization
        decay = F(x, li, ai) + off
        return np.convolve(decay, IRF, mode="same")
    Ψ
    def quadexp_iterative_model(p, x):
        # Iterative model of Grinvald and Steinberg 1974
        # See equation A-5 of DOI: 10.1016/0003-2697(74)90312-1
        outs = np.zeros((x.size, 4))
        epsilon = x[1] - x[0]
        mu, s, li, ai, off = quadexp_unpack(p)
        IRF = g(x, mu, s)
        exp = np.exp(-epsilon * li)
        epsilon_a = epsilon * ai * li
        for i in range(x.size - 1):
            outs[i + 1] = (outs[i] + 0.5 * epsilon_a * IRF[i]) * exp + 0.5 * epsilon_a * IRF[i + 1]
        return np.sum(outs, axis=-1) + off
    Ψ
    def time_function_evaluation(xsize, numevals):
        # returns array with average evaluation time of 3 different models
        # [’EMG’, ’numerical convolution’, ’iterative’]
        x = np.linspace(-20, 300, xsize)
        out = np.zeros(3)
        p = np.array([0, 1, 1 / 3, 1 / 10, 1 / 30, 1 / 100, 0.3, 0.4, 0.5, 0.6, 0.01])
        out[0] = timeit.timeit(lambda: quadexp_EMG_model(p, x), number=numevals)
        out[1] = timeit.timeit(lambda: quadexp_convolution_model(p, x), number=numevals)
        out[2] = timeit.timeit(lambda: quadexp_iterative_model(p, x), number=numevals)
        return out / numevals
    Ψ
    # average evaluation time in seconds
    #[’EMG’, ’numerical convolution’, ’iterative’]
    print(time_function_evaluation(xsize=65536, numevals=10))
    >>> [0.01594715 0.50203172 0.39133182]

## References

  * Fu  _et al._ [2017] Y. Fu, M. T. Rea, J. Chen, D. J. Morrow, M. P. Hautzinger, Y. Zhao, D. Pan, L. H. Manger, J. C. Wright, R. H. Goldsmith, and S. Jin, Selective stabilization and photophysical properties of metastable perovskite polymorphs of CsPbI3 in thin films, Chem. Mater. 29, 8385 (2017).
  * Goldsmith and Moerner [2010] R. H. Goldsmith and W. E. Moerner, Watching conformational- and photodynamics of single fluorescent proteins in solution, Nat. Chem. 2, 179 (2010).
  * Elson  _et al._ [2004] D. Elson, J. Requejo-Isidro, I. Munro, F. Reavell, J. Siegel, K. Suhling, P. Tadrous, R. Benninger, P. Lanigan, J. McGinty, C. Talbot, B. Treanor, S. Webb, A. Sandison, A. Wallace, D. Davis, J. Lever, M. Neil, D. Phillips, G. Stamp, and P. French, Time-domain fluorescence lifetime imaging applied to biological tissue, Photochem. Photobiol. Sci. 3, 795 (2004).
  * Bastiaens and Squire [1999] P. I. Bastiaens and A. Squire, Fluorescence lifetime imaging microscopy: spatial resolution of biochemical processes in the cell, Trends Cell Biol. 9, 48 (1999).
  * Gorpas  _et al._ [2019] D. Gorpas, J. Phipps, J. Bec, D. Ma, S. Dochow, D. Yankelevich, J. Sorger, J. Popp, A. Bewley, R. Gandour-Edwards, L. Marcu, and D. G. Farwell, Autofluorescence lifetime augmented reality as a means for real-time robotic surgery guidance in human patients, Sci. Rep. 9, 1187 (2019).
  * Comelli  _et al._ [2004] D. Comelli, C. D’Andrea, G. Valentini, R. Cubeddu, C. Colombo, and L. Toniolo, Fluorescence lifetime imaging and spectroscopy as tools for nondestructive analysis of works of art, Appl. Opt. 43, 2175 (2004).
  * Knight and Selinger [1971] A. Knight and B. Selinger, The deconvolution of fluorescence decay curves, Spectrochim. Acta - A: Mol. Biomol. 27, 1223 (1971).
  * Becker [2012] W. Becker, Fluorescence lifetime imaging - techniques and applications, J. Microsc. 247, 119 (2012).
  * Warren  _et al._ [2013] S. C. Warren, A. Margineanu, D. Alibhai, D. J. Kelly, C. Talbot, Y. Alexandrov, I. Munro, M. Katan, C. Dunsby, and P. M. W. French, Rapid global fitting of large fluorescence lifetime imaging microscopy datasets, PLoS ONE 8, e70687 (2013).
  * Jeansonne and Foley [1991] M. S. Jeansonne and J. P. Foley, Review of the exponentially modified gaussian (EMG) function since 1983, J. Chromatogr. Sci. 29, 258 (1991).
  * Kalambet  _et al._ [2011] Y. Kalambet, Y. Kozmin, K. Mikhailova, I. Nagaev, and P. Tikhonov, Reconstruction of chromatographic peaks using the exponentially modified gaussian function, J. Chemometrics 25, 352 (2011).
  * Purushothaman  _et al._ [2017] S. Purushothaman, S. A. S. Andrés, J. Bergmann, T. Dickel, J. Ebert, H. Geissel, C. Hornung, W. Plaß, C. Rappold, C. Scheidenberger, Y. Tanaka, and M. Yavor, Hyper-EMG: A new probability distribution function composed of exponentially modified gaussian distributions to analyze asymmetric peak shapes in high-resolution time-of-flight mass spectrometry, Int. J. Mass. Spectrom. 421, 245 (2017).
  * Golubev [2017] A. Golubev, Exponentially modified peak functions in biomedical sciences and related disciplines, Comput. Math. Methods Med. 2017, 7925106 (2017).
  * Matzke and Wagenmakers [2009] D. Matzke and E.-J. Wagenmakers, Psychological interpretation of the ex-gaussian and shifted wald parameters: A diffusion model analysis, Psychon. Bull. Rev. 16, 798 (2009).
  * Pan  _et al._ [2020] D. Pan, Y. Fu, N. Spitha, Y. Zhao, C. R. Roy, D. J. Morrow, D. D. Kohler, J. C. Wright, and S. Jin, Deterministic fabrication of arbitrary vertical heterostructures of two-dimensional ruddlesden–popper halide perovskites, Nat. Nanotechnol. 16, 159 (2020).
  * Elbaz  _et al._ [2017] G. A. Elbaz, D. B. Straus, O. E. Semonin, T. D. Hull, D. W. Paley, P. Kim, J. S. Owen, C. R. Kagan, and X. Roy, Unbalanced hole and electron diffusion in lead bromide perovskites, Nano Lett. 17, 1727 (2017).
  * Lockwood and Wasilewski [2004] D. J. Lockwood and Z. R. Wasilewski, Optical phonons in AlxGa1-xAs: Raman spectroscopy, Phys. Rev. B 70, 155202 (2004).
  * Ardekani  _et al._ [2019] H. Ardekani, R. Younts, Y. Yu, L. Cao, and K. Gundogdu, Reversible photoluminescence tuning by defect passivation via laser irradiation on aged monolayer MoS2, ACS Appl. Mater. Interfaces 11, 38240 (2019).
  * Yoo  _et al._ [2015] W. S. Yoo, K. Kang, G. Murai, and M. Yoshimoto, Temperature dependence of photoluminescence spectra from crystalline silicon, ECS J. Solid State Sci. Technol. 4, P456 (2015).
  * Zaghloul and Ali [2011] M. R. Zaghloul and A. N. Ali, Algorithm 916: Computing the faddeyeva and voigt functions, ACM Trans. Math. Softw. 38, 1 (2011).
  * Abrarov and Quine [2018] S. M. Abrarov and B. M. Quine, A rational approximation of the dawson’s integral for efficient computation of the complex error function, Appl. Math 321, 526 (2018).
  * Grinvald and Steinberg [1974] A. Grinvald and I. Z. Steinberg, On the analysis of fluorescence decay kinetics by the method of least-squares, Anal. Biochem. 59, 583 (1974).
  * Van Den Zegel  _et al._ [1986] M. Van Den Zegel, N. Boens, D. Daems, and F. De Schryver, Possibilities and limitations of the time-correlated single photon counting technique: a comparative study of correction methods for the wavelength dependence of the instrument response function, Chem. Phys. 101, 311 (1986).
  * Li  _et al._ [2020] Y. Li, S. Natakorn, Y. Chen, M. Safar, M. Cunningham, J. Tian, and D. D.-U. Li, Investigations on average fluorescence lifetimes for visualizing multi-exponential decays, Front. Phys. 8, 447 (2020).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, and P. van Mulbregt, SciPy 1.0: fundamental algorithms for scientific computing in python, Nat. Methods , 261 (2020).
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Array programming with NumPy, Nature 585, 357 (2020).
  * van der Walt  _et al._ [2011] S. van der Walt, S. C. Colbert, and G. Varoquaux, The NumPy array: A structure for efficient numerical computation, Comput. Sci. Eng. 13, 22 (2011).
  * van Rossum  _et al._ [01 ] G. van Rossum _et al._ , Python (2001–), [Online; accessed 2022-09-15].
  * Spitha  _et al._ [2020] N. Spitha, D. D. Kohler, M. P. Hautzinger, J. Li, S. Jin, and J. C. Wright, Discerning between exciton and free-carrier behaviors in ruddlesden–popper perovskite quantum wells through kinetic modeling of photoluminescence dynamics, J. Phys. Chem. C 124, 17430 (2020).
  * Manger  _et al._ [2017] L. H. Manger, M. B. Rowley, Y. Fu, A. K. Foote, M. T. Rea, S. L. Wood, S. Jin, J. C. Wright, and R. H. Goldsmith, Global analysis of perovskite photophysics reveals importance of geminate pathways, J. Phys. Chem. C 121, 1062 (2017).
  * Brown  _et al._ [1986] R. G. W. Brown, K. D. Ridley, and J. G. Rarity, Characterization of silicon avalanche photodiodes for photon correlation measurements 1: Passive quenching, Appl. Opt. 25, 4122 (1986).
  * Ziarkash  _et al._ [2018] A. W. Ziarkash, S. K. Joshi, M. Stipčević, and R. Ursin, Comparative study of afterpulsing behavior and models in single photon counting avalanche photo diode detectors, Sci. Rep. 8, 5076 (2018).


††thanks: Correspondence should be addressed to X.R.D., P.S., and Y.T.

Also at ]Huawei Technologies Co. Ltd., China.

# Koopman neural operator as a mesh-free solver of non-linear partial
differential equations

Wei Xiong  Department of Earth System Science, Tsinghua University, Beijing,
100084, China.  Xiaomeng Huang  Department of Earth System Science, Tsinghua
University, Beijing, 100084, China.  Ziyang Zhang  Huawei Technologies Co.
Ltd., China.  Ruixuan Deng<EMAIL_ADDRESS>Department of
Biomedical Engineering, Tsinghua University, Beijing, 100084, China.  Pei Sun
<EMAIL_ADDRESS>Department of Psychology & Tsinghua Laboratory of Brain
and Intelligence, Tsinghua University, Beijing, 100084, China.  Yang Tian
<EMAIL_ADDRESS>Department of Psychology & Tsinghua Laboratory
of Brain and Intelligence, Tsinghua University, Beijing, 100084, China.  [

###### Abstract

The lacking of analytic solutions of diverse partial differential equations
(PDEs) gives birth to series of computational techniques for numerical
solutions. In machine learning, numerous latest advances of solver designs are
accomplished in developing neural operators, a kind of mesh-free approximators
of the infinite-dimensional operators that map between different
parameterization spaces of equation solutions. Although neural operators
exhibit generalization capacities for learning an entire PDE family
simultaneously, they become less accurate and explainable while learning long-
term behaviours of non-linear PDE families. In this paper, we propose Koopman
neural operator (KNO), a new neural operator, to overcome these challenges.
With the same objective of learning an infinite-dimensional mapping between
Banach spaces that serves as the solution operator of target PDE family, our
approach differs from existing models by formulating a non-linear dynamic
system of equation solution. By approximating the Koopman operator, an
infinite-dimensional linear operator governing all possible observations of
the dynamic system, to act on the flow mapping of dynamic system, we can
equivalently learn the solution of an entire non-linear PDE family by solving
simple linear prediction problems. In zero-shot prediction and long-term
prediction experiments on representative PDEs (e.g., the Navier-Stokes
equation), KNO exhibits notable advantages in breaking the tradeoff between
accuracy and efficiency (e.g., model size) while previous state-of-the-art
models are limited. These results suggest that more efficient PDE solvers can
be developed by the joint efforts from physics and machine learning.

††preprint: APS/123-QED

In science and engineering, partial differential equation (PDE) is a
fundamental tool to characterize various problems [1]. However, numerous
important PDEs, such as the Navier-Stokes equation, still lack analytic
solutions [2, 3]. Consequently, the urgent needs of solving complicated PDEs
have given birth to diverse techniques for computationally approximating PDE
solutions [4].

Given $\Phi=\Phi\left(D;\mathbb{R}^{d_{\phi}}\right)$, a Banach space of
inputs, and $\Gamma=\Gamma\left(D;\mathbb{R}^{d_{\gamma}}\right)$, a Banach
space of solutions, defined on a bounded open set $D\subset\mathbb{R}^{d}$,
these PDE solvers are expected to approximate a solution operator
$\mathcal{O}$ that relates $\Phi$ with $\Gamma$ for a typically time-dependent
PDE family

$\displaystyle\partial_{t}\gamma\left(x_{t}\right)$
$\displaystyle=\left(\mathcal{L}_{\phi}\gamma\right)\left(x_{t}\right)+\eta\left(x_{t}\right),\;x_{t}\in
D\times T,$ (1) $\displaystyle\gamma\left(x_{t}\right)$
$\displaystyle=\gamma_{B},\;x_{t}\in\partial D\times T,$ (2)
$\displaystyle\gamma\left(x_{0}\right)$ $\displaystyle=\gamma_{I},\;x_{0}\in
D\times\\{0\\}.$ (3)

In Eq. (1-3), notions $\gamma_{B}$ and $\gamma_{I}$ denote the boundary and
initial conditions, set $T=\left[0,\infty\right)$ is the time domain,
differential operator $\mathcal{L}_{\phi}$ is characterized depending on
$\phi$, fixed function $\eta\left(\cdot\right)$ lives in a function space
determined by $\mathcal{L}_{\phi}$, and $\gamma\left(\cdot\right)$ is the
solution of the PDE family. This formulation gives rise to the solution
operator $\mathcal{O}:\left(\phi,\gamma_{B},\gamma_{I}\right)\mapsto\gamma$.
For convenience, we always consider fixed boundary and initial conditions in
our derivations. In application, researchers always consider a parametric
counterpart $\mathcal{O}_{\theta}\simeq\mathcal{O}$ parameterized by $\theta$
[5].

To date, diverse types of PDE solvers have been developed. In Sec. I of our
supplementary materials (SM), we present a review of classic numerical and
neural-network-based solvers. Compared with classic numerical solvers, neural-
network-based solvers, especially neural operators, are more efficient in
dealing with science and engineering questions where PDEs are increasingly
complicated [5, 6]. Therefore, our research primarily focus on neural
operator, a kind of mesh-free PDE solver [5, 6].

## Challenges.

Despite the progress in theoretical foundations [7], approximator designs [6,
8], and applications [9], there still remain numerous challenges in existing
neural operators, among which, a critical one lies in the limited capacity of
existing models to learn the long-term dynamics of complicated PDEs.

Let us rethink about the iterative update strategy of neural operators for
solving Eqs. (1-(3) [5]. As suggested in Sec. II of supplementary materials
(SM), the iteration is similar to a a dynamic system perspective where we
study the evolution mapping $\zeta:\mathbb{R}^{d_{\gamma}}\times
T\rightarrow\mathbb{R}^{d_{\gamma}}$ of an infinite-dimensional non-linear
dynamic system $\gamma_{t}=\gamma\left(D\times\\{t\\}\right)$ (notion
$\gamma\left(D\times\\{t\\}\right)$ represents that function $\gamma$ acts on
all elements in $D\times\\{t\\}$)

$\displaystyle\gamma_{t+\varepsilon}=\gamma_{t}+\int_{t}^{t+\varepsilon}\zeta\left(\gamma_{\tau},\tau\right)\mathsf{d}\tau,\;\forall
t\in T.$ (4)

For non-linear PDE families, the evolution mapping $\zeta$ is even more
intricate than the equation solution $\gamma$ itself. The infinite-dimensional
non-linear flow mapping $\zeta$ makes the accurate prediction of long-term
dynamics of $\gamma_{t}$ nearly impossible. To maintain accuracy, existing
models are required to enlarge model size or complexity, inevitably facing the
trade-off between accuracy and efficiency. Therefore, overcoming the
limitation in learning long-term PDE behaviours is necessary for optimizing
PDE solvers, which is the core objective of our research.

## Our contribution.

In this paper, we attempt to overcome the challenge of predicting long-term
PDE dynamics by proposing a new neural operator named as Koopman neural
operator (KNO). Besides learning the solution operator of an entire target PDE
family, we formalize a non-linear dynamic system of equation solution
described by Eq. (4) in the meanwhile. This characterization supports to
optimize the iterative update strategy of neural operators using dynamic
system theory. To capture intricate long-term dynamics, our model is designed
to learn the Koopman operator, an infinite-dimensional linear operator
governing all observations of a dynamic system, to act on the evolution
mapping $\zeta$ of the dynamic system of equation solution. By doing so, we
can transform the original task into an equivalent but simpler linear
prediction problem. We compare KNO with existing state-of-the-art models (e.g.
the Fourier neural operator [6]) in zero-shot and long-term prediction
experiments with representative PDEs. While previous methods suffer from the
tradeoff between accuracy and efficiency (e.g., model size), KNO is shown to
achive higher accuracy with a smaller model size.

Figure 1:  Visualization of the neural network architecture of KNO and the
instances of its prediction. (a) A KNO architecture is illustrated, where
$r=\frac{t^{\prime}-t}{\varepsilon}$ is prediction length. Note that the
layout of every part is slightly reorganized to offer a clear version. (b) Two
cases of the prediction results on the 2-dimensional Navier-Stokes equation
with a viscosity coefficient $\nu=10^{-4}$ are visualized, where KNO realizes
a 10-second forward prediction of the equation solution given initial values.
High consistency can be observed between predicted results and ground turthes.

## Computational framework.

We begin with formalizing the non-linear dynamic system of equation solution

$\displaystyle\partial_{t}\gamma_{t}=\zeta\left(\gamma_{t},t\right),\;\forall\gamma_{t}\in\mathbb{R}^{d_{\gamma}}\times
T,$ (5)

which can be either non-autonomous (i.e., the flow mapping $\theta$ associated
with $\zeta\left(\cdot,\cdot\right)$ is time-dependent) or autonomous (i.e.,
the flow mapping $\theta$ is time-independent such that
$\partial_{t}\zeta\left(\cdot,t\right)\equiv 0$) in different PDE families.
The cocycle property of the flow mapping $\theta$ [10]

$\displaystyle\theta_{t}^{t^{\prime}}=\theta_{t+\tau}^{t^{\prime}}\circ\theta_{t}^{t+\tau},\forall
t\leq t+\tau\leq t^{\prime}\in T$ (6)

defines how $\gamma$, the equation solution, evolves across adjoining time
intervals, where notion $\circ$ denotes the composition of mappings. The non-
linear, and potentially non-autonomous, dynamic system in Eqs. (5-6) makes the
long-term prediction of $\gamma_{t}$ a daunting challenge.

Our idea originates from dynamic system theory. Given with an appropriate
space defined by a set of observation function
$\mathcal{G}\left(\mathbb{R}^{d_{\gamma}}\times
T\right)=\\{\mathbf{g}|\mathbf{g}:\mathbb{R}^{d_{\gamma}}\times
T\rightarrow\mathbb{C}\\}$, we may search a family of Koopman operators
$\mathcal{K}^{t+\varepsilon}_{t}:\mathcal{G}\left(\mathbb{R}^{d_{\gamma}}\times
T\right)\rightarrow\mathcal{G}\left(\mathbb{R}^{d_{\gamma}}\times T\right)$
parameterized by time difference $\varepsilon$ to linearize the non-linear
dynamics of $\gamma_{t}$ [10]

$\displaystyle\mathcal{K}^{t+\varepsilon}_{t}\mathbf{g}\left(\gamma_{t}\right)=\mathbf{g}\left(\theta^{t+\varepsilon}_{t}\left(\gamma_{t}\right)\right)=\mathbf{g}\left(\gamma_{t+\varepsilon}\right),\;\forall
t\leq t+\varepsilon\in T.$ (7)

Therefore, after proposing an optimal approximation of $\mathbf{g}$, we can
derive a linear evolution of $\gamma_{t}$ in a new space

$\displaystyle\partial_{t}\mathbf{g}\left(\gamma_{t}\right)=\lim_{\varepsilon\rightarrow
0}\frac{\mathcal{K}^{t+\varepsilon}_{t}\mathbf{g}\left(\gamma_{t}\right)-\mathbf{g}\left(\gamma_{t}\right)}{\varepsilon},$
(8)

which is beneficial to solving Eq. (5). Eq. (7) is generally applicable to
both non-autonomous and autonomous dynamic systems. An autonomous case of Eq.
(7) is presented in Sec. II in our supplementary materials (SM).

In mathematics, the adjoint of the Koopman operator defined by Eq. (8) is the
Perron-Frobenius operator of dynamic systems [11] while the adjoint of the
associated Lie operator (see Sec. II in our supplementary materials (SM)) is
the Liouville operator of Hamiltonian dynamics [12, 13]. These properties
relate our idea with well-known theories about linear representation of
dynamic systems in statistical physics and quantum mechanics [14].

The challenges underlying Eq. (7) arise from the non-triviality of
approximating observation function $\mathbf{g}$ and estimating the Koopman
operator $\mathcal{K}^{t+\varepsilon}_{t}$. Different from existing machine-
learning-based Koopman operator models that either are limited to autonomous
dynamic systems (e.g., the case described by Eq. (7)) [15, 16, 17, 18] or
require _a priori_ knowledge about the eigenvalue spectrum (e.g, the numbers
of real and complex eigenvalues) of Koopman operator [19], we pursue to
develop an approximation applicable to both non-autonomous and autonomous
dynamic systems without strong requirement of _a priori_ knowledge [20]. In
Sec. III of supplementary materials (SM), we propose a computational framework
to learn $\mathbf{g}$ and $\mathcal{K}^{t+\varepsilon}_{t}$ based on the
Krylov subspace method [21] and the Galerkin projection. This framework is
valid as long as the target Koopman operator has a discrete spectrum [10] to
support an invariant sub-space approximation. Because the on-line learning of
a time-dependent operator $\mathcal{K}^{t+\varepsilon}_{t}$ may unavailable in
application due to its high requirement on data size and computation resource,
we further suggest $\overline{\mathcal{K}}_{\varepsilon}$, an expectation of
the Koopman operator controlled by time resolution $\varepsilon$. This
constant operator only requires economic off-line training and serves as a
reasonable counterpart of $\mathcal{K}^{t+\varepsilon}_{t}$ under ergodic
condition [22, 23]. Operator $\overline{\mathcal{K}}_{\varepsilon}$ is used as
a default setting in our framework, yet one can freely modify it to a time-
dependent operator if on-line training is accessible.

## Neural network.

Our neural network model, named as Koopman neural operator (KNO), not only
implements the mathematics defined above but also utilizes various numerical
techniques to optimize its efficiency. Below, we summarize the main
architecture of KNO.

  * (1)

_Part 1: Observation._ An input $\phi_{t}=\phi\left(D\times\\{t\\}\right)$ of
the PDE is transformed into $\mathbf{g}\left(\widehat{\gamma_{t}}\right)$ by
an observation function $\mathbf{g}\left(\cdot\right)$ parameterized by an
encoder.

  * (2)

_Part 2: Fourier transform._ To accelerate the iterative update strategy [6]
and to extract the low-frequency information (i.e., main component and
background) of the represented equation solution
$\mathbf{g}\left(\widehat{\gamma_{t}}\right)$, we apply the Fourier transform
to map $\mathbf{g}\left(\widehat{\gamma_{t}}\right)$ as
$\mathbf{g}_{\mathcal{F}}\left(\widehat{\gamma_{t}}\right)=\mathcal{F}\circ\mathbf{g}\left(\widehat{\gamma_{t}}\right)$
and truncate the Fourier series at $\omega$, a maximum number of frequency
modes. Compared with the truncated high-frequency information (i.e.,
perturbations), the remaining low-frequency information with more clear
evolution patterns may be easier to learn.

  * (3)

_Part 3: Koopman operator approximator._ Given low-frequency information
$\mathbf{g}_{\mathcal{F}}\left(\widehat{\gamma_{t}}\right)$, we define a
Hankel matrix $\widehat{\mathcal{H}}_{m\times n}$ to sample the Krylov
sequence of $\mathbf{g}_{\mathcal{F}}\left(\widehat{\gamma_{t}}\right)$ in our
computational framework (note that $n$ equals the number of accessible samples
and $m$ is a dimension of delay-embedding). Then, we train a $o\times o$
linear layer to learn a potential Koopman operator
$\overline{\mathcal{K}}_{\varepsilon}:\mathcal{G}\left(\mathbb{R}^{d_{\widehat{\gamma}}}\times
T\right)\rightarrow\widehat{\mathbb{K}}$ to predict the future state as
$\mathbf{g}_{\mathcal{F}}\left(\widehat{\gamma}_{t+r\varepsilon}\right)=\overline{\mathcal{K}}_{\varepsilon}^{r}\mathbf{g}_{\mathcal{F}}\left(\widehat{\gamma_{t}}\right)$.

  * (4)

_Part 4: Inverse Fourier transform._ Given each predicted state
$\mathbf{g}_{\mathcal{F}}\left(\widehat{\gamma}_{t+r\varepsilon}\right)$ in
_Part 3_ , we transform it from the Fourier space to
$\mathcal{G}\left(\mathbb{R}^{d_{\widehat{\gamma}}}\times T\right)$ by an
inverse Fourier transform, i.e,
$\mathbf{g}\left(\widehat{\gamma}_{t+r\varepsilon}\right)=\mathcal{F}^{-1}\circ\mathbf{g}_{\mathcal{F}}\left(\widehat{\gamma}_{t+r\varepsilon}\right)$.

  * (5)

_Part 5: High-frequency information complement._ Parallel to _Parts 2-4_ , we
do not omit high-frequency information given its non-linear effects on
$\widehat{\gamma}$. Because convolutional layers can amplify high-frequency
components [24], we train a convolutional layer $\mathcal{C}$ on the outputs
of _Part 1_ to realize a prediction of high-frequency information, i.e,
$\left[\mathbf{g}_{\mathcal{C}}\left(\widehat{\gamma}_{t+\left(r-m\right)\varepsilon}\right),\ldots,\mathbf{g}_{\mathcal{C}}\left(\widehat{\gamma}_{t+r\varepsilon}\right)\right]^{\mathsf{T}}=\mathcal{C}\left[\mathbf{g}\left(\widehat{\gamma}_{t-m\varepsilon}\right),\ldots,\mathbf{g}\left(\widehat{\gamma}_{t}\right)\right]^{\mathsf{T}}$,
where $\mathsf{T}$ is transpose.

  * (6)

_Part 6: Inverse observation._ Given two future states,
$\mathbf{g}\left(\widehat{\gamma}_{t+r\varepsilon}\right)$ and
$\mathbf{g}_{\mathcal{C}}\left(\widehat{\gamma}_{t+r\varepsilon}\right)$,
predicted by _Parts 2-4_ and _Part 5_ , we unify them by
$\mathbf{g}_{\mathcal{U}}\left(\widehat{\gamma}_{t+r\varepsilon}\right)=\mathbf{g}\left(\widehat{\gamma}_{t+r\varepsilon}\right)+\mathbf{g}_{\mathcal{C}}\left(\widehat{\gamma}_{t+r\varepsilon}\right)$.
We train an non-linear decoder to represent the inverse of observation
function
$\mathbf{g}^{-1}\left(\cdot\right)\simeq\mathbf{g}_{\mathcal{U}}^{-1}\left(\cdot\right)$
and transform
$\mathbf{g}_{\mathcal{U}}\left(\widehat{\gamma}_{t+r\varepsilon}\right)$ to
$\widehat{\gamma}_{t+r\varepsilon}$, the target state of equation solution.

Figure 2:  Experiment results of KNO. (a) Results of the mesh-independent experiment. (b) Results of the long-term prediction experiment on the 2-dimensional Navier-Stokes equation with $\nu=10^{-3}$. (c) Zero-shot prediction experiment concerning prediction length. (d) Visualization of prediction results on the 2-dimensional Navier-Stokes equation with $\nu=10^{-4}$.  Models | Parameter number |  Accuracy at $\nu=10^{-3}$ |  Accuracy at $\nu=10^{-4}$
---|---|---|---
FNO (default settings, one-unit) | 233897 |  1.78$\times 10^{-3}$ |  2.61$\times 10^{-2}$
FNO (default settings, two-unit) | 464717 | $2.12\times 10^{-4}$ | 7.72$\times 10^{-3}$
ConvLSTM (default settings) | 10001 |  6.05$\times 10^{-2}$ |  1.69$\times 10^{-1}$
U-Net (default settings) | 24950491 |  1.77$\times 10^{-3}$ |  7.90$\times 10^{-2}$
KNO ($o=24$, $f=10$, $r=12$, one-unit) | 80850 |  3.00$\times 10^{-4}$ |  9.60$\times 10^{-3}$
KNO ($o=32$, $f=10$, $r=12$, one-unit) | 206538 | 2.37$\times 10^{-4}$ | $6.66\times 10^{-3}$
Table 1: Full results of long-term prediction experiment. Accuracy is measured
by MSE. The best model is marked in bold while the second best one is marked
by the underline.

In Fig. 1(a), we illustrate the one-unit architecture of KNO. Similar to
Fourier neural operator [6], a $x$-unit KNO architecture can be produced by
cascading the copy of _Parts 2-5_ $x$ times. The full details of KNO
architecture can be seen in Sec. IV of supplementary materials (SM).

## Data sets.

We implement our experiments on the 1-dimensional Bateman–Burgers equation
[25] and the 2-dimensional Navier-Stokes equation [26]. Please see Sec. V of
supplementary materials (SM) for data description.

## Experiment designs.

We conduct five experiments to validate our model: mesh-independent
experiment, long-term prediction experiment, zero-shot experiment on
discretization granularity, zero-shot experiment on prediction length, and
ablation experiment. See Sec. IV of supplementary materials (SM) for full
descriptions.

## Details of implemented models.

Besides KNO, we implement Fourier neural operator (FNO) [6], U-shaped neural
operator (UNO) [27], convolutional LSTM (ConvLSTM) [28], and U-Net [29] for
comparison. See Sec. VI of supplementary materials (SM) for details.

## Mesh-independent experiment.

We implement experiment on the data of 1-dimensional Bateman–Burgers equation
generated under different resolutions, where FNO serves as a baseline. We
implement multiple versions of KNO and FNO with different hyper-parameters
(e.g., operator size $o$, frequency mode number $f$, the iteration number of
the Koopman operator $r=\frac{t^{\prime}-t}{\varepsilon}$ in KNO, and the
width $w$ of FNO). As illustrated in Fig. 2(a), KNO achieves almost constant
prediction error on every resolution. Compared with FNO models, the prediction
error of KNO models maintains more stable across different conditions.
Notably, a one-unit KNO only requires about $5\times 10^{3}$ parameters to
outperform a one-unit FNO with more than $2\times 10^{7}$ parameters,
suggesting the potential of KNO to break the trade-off of accuracy and
efficiency (i.e., KNO achieves higher accuracy with fewer parameters). See
Sec. VI of supplementary materials (SM) for experiment details.

## Long-term prediction experiment.

Our long-term prediction experiment is implemented on the 2-dimensional
Navier-Stokes equation, where we consider the low viscosity cases, i.e.,
$\nu\in\\{10^{-3},10^{-4}\\}$. In the cases where $\nu=10^{-3}$ and
$\nu=10^{-4}$, we consider a 40-time-interval and a 10-time-interval
prediction tasks, respectively. Table 1 reports our full results, Fig. 2(b)
visualizes the performance of experiment with $\nu=10^{-3}$, and Fig. 2(d)
visualizes instances of predicted results with $\nu=10^{-4}$. These results
suggest that KNO is optimal in characterizing the long-term evolution of
equation solution. See Sec. VI of supplementary materials (SM) for more
details.

Grid number | $2^{8}$ | $2^{9}$ | $2^{10}$ | $2^{11}$ | $2^{12}$ | $2^{13}$
---|---|---|---|---|---|---
Setting 1 |  3.820924$\times 10^{-5}$ |  3.820929$\times 10^{-5}$ |  3.820926$\times 10^{-5}$ |  3.820927$\times 10^{-5}$ |  3.820931$\times 10^{-5}$ |  3.820932$\times 10^{-5}$
Setting 2 |  1.521201$\times 10^{-5}$ |  1.521201$\times 10^{-5}$ |  1.521200$\times 10^{-5}$ |  1.521202$\times 10^{-5}$ |  1.521202$\times 10^{-5}$ |  1.521200$\times 10^{-5}$
Setting 3 |  5.403627$\times 10^{-6}$ |  5.403630$\times 10^{-6}$ |  5.403628$\times 10^{-6}$ |  5.403630$\times 10^{-6}$ |  5.403632$\times 10^{-6}$ |  5.403615$\times 10^{-6}$
Setting 4 |  3.553927$\times 10^{-6}$ |  3.553980$\times 10^{-6}$ |  3.553943$\times 10^{-6}$ |  3.553957$\times 10^{-6}$ |  3.553942$\times 10^{-6}$ |  3.553938$\times 10^{-6}$
Table 2: Zero-shot experiment (discretization granularity) on the 1-dimensional Bateman–Burgers equation. KNO settings $\left(o,f,r\right)\in\\{\left(8,10,10\right),\left(16,10,10\right),\left(32,10,16\right),\left(128,10,16\right)\\}$ are referred to as Settings 1-4. Accuracy is measured by MSE. | $\nu=10^{-3}$ | $\nu=10^{-4}$
---|---|---
| 32 grids | 64 grids | 32 grids | 64 grids
FNO (default settings, one-unit) |  1.78$\times 10^{-3}$ |  1.79$\times 10^{-3}$ |  1.41$\times 10^{-2}$ |  1.54$\times 10^{-2}$
KNO ($o=24$, $f=10$, $r=12$, one-unit) | $3.22\times 10^{-4}$ | $3.31\times 10^{-4}$ | $9.63\times 10^{-3}$ | $9.60\times 10^{-3}$
Table 3: Zero-shot experiment (discretization granularity) on the
2-dimensional Navier-Stokes equation. Accuracy is measured by MSE.

## Zero-shot experiment on discretization granularity.

As suggested by Ref. [6], a mesh-independent neural operator model can be
trained only on the data with lower resolution to predict the data with higher
resolution (referred to as zero-shot super-resolution [6]). In our research,
we validate this property by implementing a zero-shot experiment. Same as our
long-term prediction experiment, there are two data sets of the 2-dimensional
Navier-Stokes equation distinguished according to the viscosity
$\nu\in\\{10^{-3},10^{-4}\\}$ and one data set of the 1-dimensional
Bateman–Burgers equation. The settings and objective of prediction tasks on
each data set is exactly same as the long-term prediction experiment expect
that all models are trained on a lower resolution and tested on a higher
resolution. Tables 2-3 suggest the optimality of KNO in this task. See Sec. VI
of supplementary materials (SM) for experiment settings.

## Zero-shot prediction experiment concerning prediction length.

Because KNO is proposed to model the intricate dynamics of PDE solutions, it
is natural to wonder if KNO can achieve optimal performance in zero-shot
prediction with an untrained prediction length, i.e., train the model on a
prediction length $t^{\prime}-t=r\varepsilon$ and test it on a larger
prediction length $r^{\prime}\varepsilon>r\varepsilon$. This objective can be
simply realized by increasing $r=\frac{t^{\prime}-t}{\varepsilon}$, the
iteration number of the Koopman operator, during prediction. To validate the
optimality of KNO in such a task, we implement an experiment on the
2-dimensional Navier-Stokes equation with a viscosity coefficient
$\nu=10^{-3}$ and $2^{12}$ grids. All models are trained on a time interval of
$\left[0,10\right)$ to predict the interval of $\left(10,30\right]$
(supervised learning). Apart from the trained supervised learning task, models
are also required to predict interval of $\left(30,50\right]$, an untrained
time interval. As illustrated by Fig. 2(c), KNO outperforms other models in
both supervised prediction and zero-shot prediction tasks. See Sec. VI of
supplementary materials (SM) for more information.

## Ablation experiment.

As shown in Sec. VI of supplementary materials (SM), excluding or modifying
the strict combination of the observation function, the Koopman operator, and
the inverse of observation function from our model (e.g., let the whole model
directly implement forward prediction) reduces the accuracy significantly.
Therefore, the proposed architecture design of KNO is not redundant.

## Conclusion.

In summary, we develop the Koopman neural operator (KNO), a mesh-independent
neural-network-based solver of partial differential equations. The basic code
of KNO is provided in Sec. IV of supplementary materials (SM) and the official
toolbox can be seen in Ref. [30]. Compared with the existing state-of-the-art
approaches, such as the Fourier neural operator, KNO exhibits a higher
capacity to capture the long-term evolution of equation solution and maintain
accuracy across different mesh resolutions. This property suggests the
potential of KNO to be considered as a basic unit for constructing large-scale
frameworks of solving complex physics equations with non-linear dynamics.
Moreover, the robust mesh-independence enables KNO to be trained on the low
resolution data before being applied on the high resolution data. This
property may create more possibilities to reduce the computational cost of
partial differential equation solving.

## Acknowledgements.

This project is supported by the Artificial and General Intelligence Research
Program of Guo Qiang Research Institute at Tsinghua University (2020GQG1017)
as well as the Tsinghua University Initiative Scientific Research Program.

## References

  * Debnath and Debnath [2005] L. Debnath and L. Debnath,  _Nonlinear partial differential equations for scientists and engineers_ (Springer, 2005).
  * Gockenbach [2005] M. S. Gockenbach,  _Partial differential equations: analytical and numerical methods_ , Vol. 122 (Siam, 2005).
  * Tanabe [2017] H. Tanabe,  _Functional analytic methods for partial differential equations_ (CRC Press, 2017).
  * Mattheij  _et al._ [2005] R. M. Mattheij, S. W. Rienstra, and J. T. T. Boonkkamp,  _Partial differential equations: modeling, analysis, computation_ (SIAM, 2005).
  * Li  _et al._ [2020a] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, Neural operator: Graph kernel network for partial differential equations, arXiv preprint arXiv:2003.03485 (2020a).
  * Li  _et al._ [2020b] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, Fourier neural operator for parametric partial differential equations, arXiv preprint arXiv:2010.08895 (2020b).
  * Kovachki  _et al._ [2021] N. Kovachki, S. Lanthaler, and S. Mishra, On universal approximation and error bounds for fourier neural operators, Journal of Machine Learning Research 22, Art (2021).
  * Li  _et al._ [2022] Z. Li, D. Z. Huang, B. Liu, and A. Anandkumar, Fourier neural operator with learned deformations for pdes on general geometries, arXiv preprint arXiv:2207.05209 (2022).
  * Guibas  _et al._ [2021] J. Guibas, M. Mardani, Z. Li, A. Tao, A. Anandkumar, and B. Catanzaro, Efficient token mixing for transformers via adaptive fourier neural operators, in  _International Conference on Learning Representations_ (2021).
  * Brunton  _et al._ [2022] S. L. Brunton, M. Budisic, E. Kaiser, and J. N. Kutz, Modern koopman theory for dynamical systems, SIAM Review 64, 229 (2022).
  * Lasota and Mackey [1985] A. Lasota and M. C. Mackey,  _Probabilistic properties of deterministic systems_ (Cambridge university press, 1985).
  * Gaspard  _et al._ [1995] P. Gaspard, G. Nicolis, A. Provata, and S. Tasaki, Spectral signature of the pitchfork bifurcation: Liouville equation approach, Physical Review E 51, 74 (1995).
  * Gaspard [2005] P. Gaspard, Chaos, scattering and statistical mechanics, Chaos (2005).
  * Lasota and Mackey [1998] A. Lasota and M. C. Mackey,  _Chaos, fractals, and noise: stochastic aspects of dynamics_ , Vol. 97 (Springer Science & Business Media, 1998).
  * Takeishi  _et al._ [2017] N. Takeishi, Y. Kawahara, and T. Yairi, Learning koopman invariant subspaces for dynamic mode decomposition, Advances in Neural Information Processing Systems 30 (2017).
  * Azencot  _et al._ [2020] O. Azencot, N. B. Erichson, V. Lin, and M. Mahoney, Forecasting sequential data using consistent koopman autoencoders, in  _International Conference on Machine Learning_ (PMLR, 2020) pp. 475–485.
  * Otto and Rowley [2019] S. E. Otto and C. W. Rowley, Linearly recurrent autoencoder networks for learning dynamics, SIAM Journal on Applied Dynamical Systems 18, 558 (2019).
  * Alford-Lago  _et al._ [2022] D. J. Alford-Lago, C. W. Curtis, A. T. Ihler, and O. Issan, Deep learning enhanced dynamic mode decomposition, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 033116 (2022).
  * Lusch  _et al._ [2018] B. Lusch, J. N. Kutz, and S. L. Brunton, Deep learning for universal linear embeddings of nonlinear dynamics, Nature communications 9, 1 (2018).
  * Macesic  _et al._ [2018] S. Macesic, N. Crnjaric-Zic, and I. Mezic, Koopman operator family spectrum for nonautonomous systems, SIAM Journal on Applied Dynamical Systems 17, 2478 (2018).
  * Saad [2011] Y. Saad,  _Numerical methods for large eigenvalue problems: revised edition_ (SIAM, 2011).
  * Arbabi and Mezic [2017] H. Arbabi and I. Mezic, Ergodic theory, dynamic mode decomposition, and computation of spectral properties of the koopman operator, SIAM Journal on Applied Dynamical Systems 16, 2096 (2017).
  * Cornfeld  _et al._ [2012] I. P. Cornfeld, S. V. Fomin, and Y. G. Sinai,  _Ergodic theory_ , Vol. 245 (Springer Science & Business Media, 2012).
  * Park and Kim [2022] N. Park and S. Kim, How do vision transformers work?, arXiv preprint arXiv:2202.06709 (2022).
  * Benton and Platzman [1972] E. R. Benton and G. W. Platzman, A table of solutions of the one-dimensional burgers equation, Quarterly of Applied Mathematics 30, 195 (1972).
  * Wang [1991] C. Wang, Exact solutions of the steady-state navier-stokes equations, Annual Review of Fluid Mechanics 23, 159 (1991).
  * Rahman  _et al._ [2022] M. A. Rahman, Z. E. Ross, and K. Azizzadenesheli, U-no: U-shaped neural operators, arXiv preprint arXiv:2204.11127 (2022).
  * Wang  _et al._ [2020] R. Wang, K. Kashinath, M. Mustafa, A. Albert, and R. Yu, Towards physics-informed deep learning for turbulent flow prediction, in  _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_ (2020) pp. 1457–1466.
  * Ronneberger  _et al._ [2015] O. Ronneberger, P. Fischer, and T. Brox, U-net: Convolutional networks for biomedical image segmentation, in  _International Conference on Medical image computing and computer-assisted intervention_ (Springer, 2015) pp. 234–241.
  * Xiong  _et al._ [2023] W. Xiong, M. Ma, P. Sun, and Y. Tian, Koopmanlab: A module of koopman neural operator family for solving partial differential equations, https://github.com/Koopman-Laboratory/KoopmanLab (2023).