# Multifractality of ab initio wave functions in doped semiconductors

Edoardo G. Carnio  Physikalisches Institut, Albert-Ludwigs-Universität
Freiburg, Hermann-Herder-Straße 3, 79104 Freiburg, Germany  Nicholas D.M. Hine
Department of Physics, University of Warwick, Coventry CV4 7AL, United Kingdom
Rudolf A. Römer

###### Abstract

In Refs. [1, 2] we have shown how a combination of modern linear-scaling DFT,
together with a subsequent use of large, effective tight-binding Hamiltonians,
allows to compute multifractal wave functions yielding the critical properties
of the Anderson metal-insulator transition (MIT) in doped semiconductors. This
combination allowed us to construct large and atomistically realistic samples
of sulfur-doped silicon (Si:S). The critical properties of such systems and
the existence of the MIT are well known, but experimentally determined values
of the critical exponent $\nu$ close to the transition have remained different
from those obtained by the standard tight-binding Anderson model. In Ref. [1],
we found that this “exponent puzzle” can be resolved when using our novel _ab
initio_ approach based on scaling of multifractal exponents in the realistic
impurity band for Si:S. Here, after a short review of multifractality, we give
details of the multifractal analysis as used in [1] and show the obtained
_critical_ multifractal spectrum at the MIT for Si:S.

###### keywords:

Multifractality , Anderson localization , Semiconductors

††journal: Physica E

##  1 Introduction

The Anderson metal-insulator transition (MIT) [3] is one of the fundamental
manifestations of the quantum mechanical nature of disordered materials [4, 5,
6]. In his 1958 publication [3], Anderson studies the localization of
electrons in doped semiconductors. The existence of the MIT in these materials
was later confirmed by measuring the scaling of the conductance when
increasing the dopant concentration beyond a critical value [7, 8, 9, 10]. The
critical properties of the MIT, such as the exponent $\nu$ of the
conductivity, should be universal quantities [11, 12]. For classical waves
[13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23] and cold atom systems [24, 25,
26, 27, 28, 29, 12] results agree well with many of the non-interacting
Anderson model estimates [30, 31]. However, in experiments with
semiconductors, $\nu$ is found to vary with sample-specific properties, namely
the dopant concentration around the transition point, the homogeneity of the
doping, and the purity of the sample itself [32]. The term “exponent puzzle”
[8, 9] was hence coined to describe this inability to characterize the
Anderson transition in terms of a single, universal value for $\nu$.

Recently, we presented a study [1] that moves beyond the paradigmatic, highly-
simplified, tight-binding Anderson model, and employs atomistically correct ab
initio simulations [33, 34] of the doped semiconductor Si:S [35]. With this
approach we observe how the impurity band (IB) forms and eventually merges
with the conduction band upon increasing the dopant concentration. We then
exploit the multifractal nature of the (near-)critical electronic wave
functions as a basis for a finite-size scaling analysis that aims to retrieve
the critical properties of the transition in the thermodynamic limit. In order
to reach sufficiently large system sizes, we devised a hybrid approach:
linear-scaling DFT calculations using the ONETEP code [36], on prototype
systems of $8\times 8\times 8$ diamond-cubic unit cells ($4096$ atoms) to
construct catalogs of local Hamiltonian blocks to describe the
no/single/double dopant situation. For each concentration of impurities and
disorder realization, we then built much larger, effective tight-binding
Hamiltonians $H$ and overlap matrices $O$ from these catalogs for system size
$L$, and solved the large generalized eigenvalue problem
$H\psi_{j}=\epsilon_{j}O\psi_{j}$ for eigenenergies $\epsilon_{j}$ and
eigenvectors $\psi_{j}(\mathbf{r})$, with $\mathbf{r}=(x,y,z)$ for coordinates
$x,y,z=1,\ldots,L$.

In this paper, we show in detail how to perform the multifractal analysis of
the $\psi_{j}$ and present the resulting _critical_ multifractal spectrum. Our
results suggest that it is different from the spectrum in the Anderson model
[2]. In addition, we give further details in the finite-size scaling analysis
needed to ascertain the existence and the properties of the MIT when
$L\rightarrow\infty$.

##  2 Some of the basics of multifractals

For the Anderson transition from a metal to an insulator upon increasing the
disorder [3], the absence of length scales at criticality means that the wave
function intensity $|\psi(\boldsymbol{r})|^{2}$ at the critical disorder is
_self-similar_ [37, 4, 6]. It needs to have a “filamentary” structure [38],
i.e. it needs to be extended throughout the volume, a property of the metal
phase, but also to occupy only an infinitesimal fraction of it, a property of
the localised phase. This structure allows the critical phase to be
continuously connected to both the extended and localised phases. In
conjunction with its self-similar property, the critical wave function
qualifies as a _fractal_ [39], at least as long as we can disregard the lower
limit imposed by the lattice spacing $a$. Castellani _et al._ [40] realised,
based on the earlier work of Wegner [41], that the critical wave function is
not a simple fractal, but rather an “interwoven family” of fractals, each with
its own dimension and distribution. Such an object is a _multifractal_ [42,
43].

###  2.1 Self-similarity

Let us consider a system occupying a finite region of space
$\Sigma\subset\mathbb{R}^{D}$ with a local density $\rho(\boldsymbol{r})$.
Following Ref. [44], we define the pair correlation function,
$g(\boldsymbol{r})=\langle\rho(\boldsymbol{r}+\boldsymbol{r}^{\prime})\rho(\boldsymbol{r}^{\prime})\rangle_{\boldsymbol{r}^{\prime}}$,
which gives the probability that two points separated by $\boldsymbol{r}$ both
belong to the region $\Sigma$. For simplicity, we now assume that the
correlation function is isotropic, $g(\boldsymbol{r})=g(r)$. In the absence of
length scales, $g$ obeys _homogeneity laws_ (or _scale-invariance_) with
respect to a _resolution_ or _coarse-graining_ $\lambda$. More specifically,
if we rescale lengths as $r\rightarrow r^{\prime}=\lambda r$ we have that
$g(r^{\prime})=\lambda^{\kappa}g(r)$, where $\kappa$ is a _homogeneity
exponent_. The solution to this equation is given by a _power-law behaviour_ ,
$g(r)\propto r^{\kappa}$. If we then fix $r$ as the reference length scale,
and $g(r)=1$, we see that

$g(\lambda)=\lambda^{\kappa}$ (1)

translates the concept of self-similarity into a mathematical relation. We can
hence describe fractal objects as self-similar structures whose observed
spatial extent (e.g. volume) depends, with a power-law behaviour, on the
resolution at which we look at it. For fractals originating from a
mathematical relation, the dependence on the resolution can extend over an
infinite range. For fractals appearing in physical systems, instead, the range
of $\lambda$ is usually limited by macro- and/or microscopic scales. A very
comprehensive list of examples can be found, e.g., in [45].

###  2.2 Measures, fractals and multifractals

Let $\psi(\boldsymbol{r})$ be the wave function of an electron in a $L\times
L\times L$ volume. The modulus square $|\psi(\boldsymbol{r})|^{2}$ defines a
normalised measure on this volume, which we partition in boxes of linear size
$l=\lambda L$. The number of boxes will then be $\lambda^{-d}$, where $d=3$ is
the Euclidean dimension of the support of the system. The probability of
finding the electron in box $\mathcal{B}_{i}$ is the _box-probability_
$\mu_{i}=\int_{\mathcal{B}_{i}}|\psi(\boldsymbol{r})|^{2}\operatorname{d}\boldsymbol{r}$.
We can then compute the _fractal dimension_ $D$ of the system by counting the
number of boxes where the box-probability does not vanish:111It is customary
to use $\sim$ to indicate that the proportionality constant is independent of
the resolution and can thus be ignored. This constant might appear, for
instance, when the boxes, whichever their shape, do not perfectly cover the
system. Since in our analysis we are covering boxes with boxes, most of the
relations in this section are actually equalities.
$N(\lambda)\sim\lambda^{-D}$. Because the electron can access any portion of
the volume, i.e. there is no region of space with vanishing probability, we
trivially conclude that the fractal dimension is $D=d=3$.

Compared to the fractal dimension, more insightful is actually the study of
the powers of the box-probability $\mu^{q}_{i}$, which is the idea behind
_multifractal_ analysis. _If_ the wave function is a multifractal, we expect
to see the power-law behaviour of (1):

$\langle\mu^{q}\rangle_{L}\sim\lambda^{D+\tau_{q}}\,,$ (2)

where $\langle\ldots\rangle_{L}$ denotes the average over all boxes in the
volume. Equivalently, we can introduce the _partition sum_ $R_{q}(\lambda)$
(also the _generalised inverse participation ratio_) as
$\langle\mu^{q}\rangle_{L}=\lambda^{D}R_{q}(\lambda)$ and write

$R_{q}(\lambda)=\sum_{i}\mu_{i}^{q}\sim\lambda^{\tau_{q}}\,.$ (3)

The _mass exponents_ $\tau_{q}$ describe the scaling behaviour of the moments
and do not depend on $\lambda$.

Let us stress again that multifractality holds if, in the power-law relation
of Eq. (2), we have $\tau_{q}\neq 0$ for a finite range of $\lambda$: the box
size $l$ should be smaller than the system size, but also larger than the
microscopic scale $a$, usually equal to the lattice spacing in the Anderson
case. At the same time, for critical states at the Anderson transition, the
system size is much smaller than the correlation length $\xi$, such that

$a\ll l<L\ll\xi$ (4)

Additionally, the wave function is truly critical (and hence multifractal)
only in the thermodynamic limit, where both $L$ and $\xi$ diverge. Thus
$\tau_{q}$ is defined in the limit $\lambda\rightarrow 0$. For finite systems,
instead, we choose states and coarse-grainings that satisfy (4). In this case,
we can _estimate_ $\tau_{q}$ by fitting the slope of $\log R_{q}(\lambda)$
versus $\log\lambda$. We are assuming here that multifractality survives in
finite systems [46], and postpone the discussion of this non-trivial
assumption to section 5.

From (2) and the normalisation of the wave function, it is possible to show
that $\tau_{0}=-D$ and $\tau_{1}=0$. This implies that we can generalise the
definition of the fractal dimension to the _anomalous dimensions_ $D_{q}$ such
that $D_{0}=D$ and _mass strengths_ $\tau_{q}=D_{q}\,(q-1)$. In the case of a
simple fractal $D_{q}\equiv D$, while for a multifractal $D_{q}$ has a non-
trivial dependence on $q$. The deviation from the simple-fractal case is
captured by the _anomalous scaling exponent_
$\Delta_{q}=(D_{q}-d)(q-1)=\tau_{q}-d(q-1)$.

Figure 1:  Wave function of the impurity-band state closest to
$\varepsilon_{\text{F}}$ for an exemplary systems of $4067$ Si atoms and $29$
S impurities. (a) Shows the state calculated with ONETEP, while in (b)
calculated in the effective tight-binding model. We have represented the top
90% wave function values $\psi_{i}$ with spheres of volume proportional to
$|\psi_{i}|^{2}$. Opacity and colour are proportional to
$-\log_{L}|\psi_{i}|^{2}$, with $L=16$ here, so that lower (higher) values are
in red transparent (violet solid). Orange dots indicate the position of the
impurities. We use $a$ to denote the Si lattice parameter.

###  2.3 The multifractal spectrum

The scaling of the moments $R_{q}$, yielding $\tau_{q}$, is enough to fully
characterise the multifractal nature of the wave function. Now we present an
equivalent description of the multifractal that will be useful, in the
sections that follow, to validate our results and compare them to the 3D
Anderson model. This description is founded on a _multifractal measure_ [47],
a distribution such that, around each box $\mathcal{B}_{i}$,
$\mu_{i}=\lambda^{\alpha_{i}}$. The set of boxes with
$\alpha_{i}\in[\alpha,\alpha+\operatorname{d}\alpha]$, then, constitutes a
simple fractal with dimension $f(\alpha)$, such that the number of said boxes
is

$N_{\lambda}(\alpha)\sim\lambda^{-f(\alpha)}\qquad\text{and}\qquad\alpha_{i}=\frac{\log\mu_{i}}{\log\lambda}\,.$
(5)

This is the formalisation of the idea of Castellani [40] that the multifractal
is composed of different simple fractals. We re-express the partition sum of
Eq. (3) as

$R_{q}(\lambda)=\sum_{i}\mu^{q}_{i}=\sum_{i}\lambda^{q\alpha_{i}}=\int
N(\alpha)\lambda^{q\alpha}\operatorname{d}\alpha\sim\int\lambda^{q\alpha-f(\alpha)}\operatorname{d}\alpha\,.$
(6)

For small $\lambda$, we can use the saddle point approximation and find that
the biggest contribution in the integral (6) comes from the value of $\alpha$
that maximises (since $\lambda<1$) the argument of the exponential, i.e. the
$\alpha_{q}$ such that $f^{\prime}(\alpha_{q})=q$. We can then write, from
(3), $\tau_{q}=q\alpha_{q}-f(\alpha_{q})$. If we identify
$f_{q}=f(\alpha_{q})$ we can see that $(q,\tau_{q})$ and $(\alpha_{q},f_{q})$
are related by a Legendre transformation

$f_{q}=q\,\alpha_{q}-\tau_{q}\qquad\text{and}\qquad\alpha_{q}=\frac{\operatorname{d}\tau_{q}}{\operatorname{d}q}\,.$
(7)

It can be proven, e.g. in [48], that $\tau_{q}$ is a monotonically increasing
function in $q$, which implies that $\alpha_{q}>0,\forall q$.

Figure 2: Multifractal spectrum $f(\alpha)$ of the highest-occupied molecular
orbital wave function of the Onetep prototype, computed for $q$ from $-2$ to
$5$ in steps of $0.1$ (increasing from right to left). Dashed lines indicate
the functions $f_{0}\equiv D$ and $f_{1}(\alpha)=\alpha$.

We can combine _singularity strengths_ $\alpha_{q}$ and the _singularity
spectrum_ $f_{q}$ to obtain the _multifractal spectrum_ $f(\alpha)$. This
function is equivalent to the generalised dimensions $D_{q}$ in characterising
the multifractal, and in the case of a simple fractal analogously reduces to
the point $(D,D)$ in a $(\alpha,f(\alpha))$ plot. As shown in the example of
Fig. 2, $f(\alpha)$ is a convex function reaching its maximum at $\alpha_{0}$
with a value $f_{0}=\tau_{0}=D$. From (7) we further notice that
$f_{1}=\alpha_{1}$, since $\tau_{1}=0$. The spectrum is therefore tangential
to the functions $f_{0}(\alpha)\equiv D$ and $f_{1}(\alpha)=\alpha$.

###  2.4 Symmetry of the multifractal spectrum

Using the nonlinear $\sigma$ model, Mirlin _et al._ [49] have analytically
proven that at criticality the multifractal exponents (7) satisfy the exact
symmetry relation

$\alpha_{q}+\alpha_{1-q}=2d\qquad f_{1-q}=f_{q}+d-\alpha_{q}\;.$ (8)

Assuming the universality of the critical properties at the Anderson
transition, this result is expected to generally hold for the Wigner-Dyson
symmetry classes [6]. Indeed, Eq. (8) was confirmed numerically for different
systems, including the 3D Anderson model [50, 51] and in experiments [20]

##  3 Multifractal analysis of the wave function

We are mainly interested in the singularity strengths $\alpha_{q}$, which,
together with $\tau_{q}$ and $\Delta_{q}$ are called multifractal exponents
(MFE). In this section we recast the exponents derived in Sec. 2 in a form
that is more convenient for numerical calculations, mostly by reducing the
loss of precision. We then extend our definitions to include a disorder
ensemble average.

Before moving to the determination of these exponents, let us briefly discuss
their expected values [51, 52]. In the case of very low disorder, the wave
function intensities $|\psi(\mathbf{r})|^{2}$ will be nearly plane-wave like,
i.e. extended throughout the whole of the volume $L^{3}$. Hence, as shown in
section 2.2, $\alpha_{q}=D=d=3$ for all $q$. This implies that
$\tau_{q}=d(q-1)$, $\Delta_{q}=0$ and $f_{q}=d$. Hence $f(\alpha)$ is seen to
contract, and eventually converge, to one point $f(d)=d$. On the other hand,
for very large disorder, well into the insulating/localized regime,
$|\psi(\mathbf{r})|^{2}$ can be seen as localized on a few (single) sites
only. In this case $\alpha_{q}\rightarrow\infty$ for $q\leq 0$, and
$\alpha_{q}\rightarrow 0$ for $q>0$. Similarly, $\tau_{q}\rightarrow-\infty$,
$\Delta_{q}\rightarrow-\infty$ for $q\leq 0$, and $\tau_{q}=0$,
$\Delta_{q}=d(1-q)$ for $q>0$. The $f(\alpha)$ spectrum broadens and in the
limit of strong localization, converges to the points $f(0)=0$ and
$f(\infty)=d$. At criticality, the behaviour is even richer [52], with a non-
parabolic $f(\alpha)$ [53] and $\alpha_{0}=4.048(4.045,4.050)$,
$\alpha_{1}=1.958(1.953,1.953)$ the current best estimates for the non-
interacting 3D Anderson model [52]. As a rough guide for the following
sections, the $\alpha_{q}$ tend towards $d=3$ for extended states as weak
disorder, while $\alpha_{0}$ increases without bounds in the insulating
regime.

###  3.1 Numerical calculation

Following [54], it is convenient to define, from (3) and (7), the auxiliary
quantity

$S_{q}(\lambda)=\frac{\operatorname{d}R_{q}(\lambda)}{\operatorname{d}q}=\sum_{i}\mu_{i}^{q}\log\mu_{i}\,.$
(9)

This ratio can be interpreted as an average with respect to the measure
defined by $\mu^{q}$. The latter is also called _$q$ -microscope_, because it
increases the large (small) box-probabilities for $q>0$ ($q<0$). A
computationally-friendly formulation of the MFE reads

$\tau_{q}=\lim_{\lambda\rightarrow 0}\frac{\log
R_{q}(\lambda)}{\log\lambda}\,,\qquad\Delta_{q}=\tau_{q}-q(d-1)\,,\qquad\alpha_{q}=\lim_{\lambda\rightarrow
0}\frac{S_{q}(\lambda)}{R_{q}(\lambda)\log\lambda}\,.$ (10)

To comply with (4), we choose $\lambda\leq 1/2$, namely we consider boxes of
linear size up to $l\leq L/2$. We coarse-grain the wave function using the
partitioning scheme proposed by [55]. Here, the box size $l$ can take any
integer value (up to $L/2$), so that $\lambda^{-1}=L/l$ can take non-integer
values. This is achieved by first periodically replicating the original
system, such that it can be exactly covered by an integer number of boxes, and
then by averaging over the possible equivalent box origins. The increased
number of available box sizes translates, in the linear fits, in reduced
uncertainties in the estimated slopes.

###  3.2 Ensemble averaging

So far we have computed the multifractal properties of a single wave function.
The multifractal analysis of the Anderson transition is usually performed by
taking an average over the disorder realisations. The definitions of the MFE
can be extended by defining the _ensemble average_ of the partition sum as
$\langle R_{q}(\lambda)\rangle\sim\lambda^{\tau_{q}^{\text{ens}}}$, such that
$\tau_{q}^{\text{ens}}=\lim_{\lambda\rightarrow 0}{\log\langle
R_{q}(\lambda)\rangle}/{\log\lambda}$. We then proceed to take the Legendre
transform and define

$f_{q}^{\text{ens}}=q\alpha_{q}^{\text{ens}}-\tau_{q}^{\text{ens}}\qquad\text{and}\qquad\alpha_{q}^{\text{ens}}=\frac{\operatorname{d}\tau_{q}^{\text{ens}}}{\operatorname{d}q}=\lim_{\lambda\rightarrow
0}\frac{\langle S_{q}(\lambda)\rangle}{\langle
R_{q}(\lambda)\rangle\log\lambda}\,.$ (11)

Notice that, in the ensemble average of $\alpha_{q}$, the $q$-microscope
$\mu_{i}^{q}$ is normalised by $\langle R_{q}\rangle$, namely the averaged
partition sum of the wave function. If we normalised the $\mu_{i}^{q}$ terms
for every wave function we would obtain the _typical_ average
$\alpha_{q}^{\text{typ}}=\langle S_{q}/R_{q}\rangle/\log\lambda$
($\lambda\rightarrow 0$). While in the ensemble average all wave functions,
including rare events, are equally weighted, the typical average is dominated
by the behaviour of “typical” wave functions. The presence of rare events
translates in the appearance of negative fractal dimensions (see Sec. 4.1), a
feature of the $f(\alpha)$ that is best captured by ensemble averaging [6].

##  4 Results for eigenstates of the effective Hamiltonians

For finite systems, the multifractal exponents are computed, as explained in
Sec. 2.2, by estimating the slope of a $\log\langle R_{q}(\lambda)\rangle$ vs.
$\log\lambda$ plot, in the case of $\tau_{q}^{\text{ens}}$. Accordingly, the
statistical uncertainties at fixed $\lambda$ have to be multiplied by a factor
$\log\lambda$.222The standard deviation $\sigma_{\alpha_{q}}$ associated to
$\alpha_{q}^{\text{ens}}$ (at a fixed $\lambda$) is related to the standard
deviations $\sigma_{S_{q}}$ and $\sigma_{R_{q}}$ and the covariance
$\operatorname{cov}(S_{q},R_{q})$ via propagation of the variance [52]
$\sigma_{\alpha_{q}}=\alpha_{q}^{\text{ens}}\sqrt{\frac{\sigma^{2}_{S_{q}}}{\langle
S_{q}\rangle^{2}}+\frac{\sigma^{2}_{R_{q}}}{\langle
R_{q}\rangle^{2}}-2\frac{\operatorname{cov}(S_{q},R_{q})}{\langle
S_{q}\rangle\langle R_{q}\rangle}}\quad\text{ and
}\quad\sigma_{\tau_{q}}=\frac{\sigma_{R_{q}}}{\langle
R_{q}\rangle\log\lambda}\,,$ (12) so that
$\sigma_{f_{q}}^{2}=\sigma_{\alpha_{q}}^{2}+\sigma_{\tau_{q}}^{2}$. The
standard error of the mean is obtained by dividing every standard deviation by
$\sqrt{\mathcal{N}}$, where $\mathcal{N}$ is the number of available
realisations. In Fig. 3 we show the average singularity spectrum for the
ensemble of $10\,648$ atoms with $140$ impurities, a system that is critical
at energy close to $-0.249\text{\,}\mathrm{e}\mathrm{V}$ [1]. The increase in
the ensemble size does not change the shape of the spectrum significantly, but
has the effect of reducing the error bars on the data points. In particular,
it is the calculation of the average $\langle R_{q}\rangle$ and $\langle
S_{q}\rangle$ in Eq. (11) that benefits from larger ensembles, since smaller
error bars in the data used for the linear fits results in smaller
uncertainties on the fit parameters, see Fig. 3(b). We quantitatively report
the quality of said fits in the lower panel of Fig. 3(a), where we show the
linear correlation coefficient $r^{2}$ and the $p$ value. As noted in [50],
while $r^{2}\approx 1$ indicates a good linear behaviour, small $p$ values
suggest that the uncertainties on the data point are too small to support the
deviation from the linear behaviour we are fitting. This is likely due to the
limited number of realisations available for the ensemble averaging. For
comparison, at the end of the two branches, i.e. for large $|q|$ values, error
bars are larger and hence the quality-of-fit increases again.

###  4.1 Negative fractal dimensions

Error bars increase on the two ends of the spectrum, for negative $q$ (right)
and positive (left). For $q<0$, the $q$-microscope increases the weight of
small values of the wave function, which are more sensitive to numerical
fluctuations from the diagonalization. The other end of the spectrum ($q>0$,
left) describes instead the presence of _rare_ critical wave functions with
small values of $\alpha$ and hence large $|\psi|^{2}\sim L^{-\alpha}$ [53].
The set of these values scales with a negative fractal dimension $f(\alpha)$,
which means that their occurrence frequency vanishes in the
$L\rightarrow\infty$ limit. This is a known effect arising from ensemble
averaging [56] and known since the pioneering work of Mandelbrot [39]. This
finite-size effect is further observed and commented in [50], where larger
systems are accessible and studied. In comparison, the single state used to
produce Fig. 2 does not show said rare boxes with large probability
amplitudes, as indicated by $f_{q}>0$.

###  4.2 Width of the multifractal spectrum

Let us comment on the width of the distribution in Fig. 3, as compared to the
Anderson model studied in [50]. A narrow $f(\alpha)$ spectrum implies that
extreme values (either large or small) occur less frequently. This means that,
in our case, the average state near criticality in our model looks more
homogeneous or extended than in the Anderson model. In Fig. 4, instead, we
show the singularity spectrum for a system of $10\,648$ atoms, which, for
$230$ impurities, is close to criticality at the Fermi energy
$\varepsilon_{\text{F}}=0$ and deeper in the impurity band at
$-0.320\text{\,}\mathrm{e}\mathrm{V}$ (estimated from [1]). While both spectra
are narrower than the Anderson model, the critical wave function at
$\varepsilon_{\text{F}}$ appears on average more extended than deeper in the
impurity band.

Our results are reminiscent of those found by [57] for the power-law random
banded matrix (PRBM) model, which describes a 1D chain with random long-range
hopping decaying as $r^{-\alpha}$ over distances larger than a band width $b$.
For the critical value $\alpha=1$, the model undergoes an Anderson transition
for any value of $b$, which parametrises a family of critical models that can
be studied from the weak- ($b\gg 1$) to the strong-coupling ($b\ll 1$) regime.
For $b\gg 1$ the model shows a “quasi-metallic” behaviour, where the critical
wave functions has statistical properties similar to the delocalised phase.
The singularity spectrum becomes correspondingly narrower with a parabolic
shape, a regime called _weak multifractality_. In this case the multifractal
spectrum follows the parabolic _approximation_ [48]:

$f(\alpha)\simeq
d-\frac{(\alpha-\alpha_{0})^{2}}{4(\alpha_{0}-d)}\qquad\text{and}\qquad\alpha_{0}=d+\gamma\,.$
(13)

In Fig. 3 we fit our full-ensemble data to (13) to find the estimate value
$\alpha_{0}\approx 3.55$. We only report an approximate value without
uncertainty because the parabolic behaviour is an approximation and does not
necessarily hold for the whole spectrum, since $f(\alpha)$ is defined only for
positive $\alpha$. In fact, for $q_{c}=(d+\gamma)/2\gamma$ [6] we have
$\alpha_{q_{c}}=0$ and $f(0)$ is finite (a _termination point_), whereas our
results in Fig. 3 and those in [50] seems to suggest that
$f(\alpha)\rightarrow-\infty$ in the limit $\alpha\rightarrow 0$ (hence no
termination point).

###  4.3 Symmetry of the multifractal spectrum

In Fig. 3 we also show the symmetrised spectrum obtained by computing and
plotting $\alpha_{1-q}$ and $f_{1-q}$ from Eq. (8). As expected from the
previous paragraphs, the uncertainty on the data points increase at the
extremities. Within these error bars, the spectra are in good agreement with
each other. We verify the same symmetry relation also for the spectra in Fig.
4. While at $-0.320\text{\,}\mathrm{e}\mathrm{V}$ there is excellent agreement
between the spectra, at the Fermi energy there is a slightly higher
discrepancy, especially at extreme values of $q$. Since this discrepancy would
be resolved by taking two standard deviations as confidence intervals, instead
of one, we cannot attribute this discrepancy to any specific underlying
physical factor or systematic error.

Figure 3:  (a) top panel: singularity spectrum for $L=22$ and
$N_{\text{S}}=140$, sampled for values of $q$ from $-2$ to $5$ in steps of
$1/4$ (increasing from right to left) at energy
$-0.249\text{\,}\mathrm{e}\mathrm{V}$. Blue diamonds show the results for the
ensemble of the first 100 disorder realisations, while red circles indicate
the results from all available realisations (597). Simple error bars, without
data point, indicate the symmetrised spectrum to the full ensemble. Dashed
lines indicate the functions $f_{0}\equiv D$ and $f_{1}(\alpha)=\alpha$. The
dotted line indicates the spectrum for the Anderson model at criticality,
reproduced from [50], while the dot-dashed line indicates the fit to the
parabolic approximation (13). (a) bottom panel: linear correlation coefficient
$r^{2}$ (red) and quality of fit $p$ (black) for the linear fits used to
extrapolate the thermodynamic limit of $\alpha_{q}$ (diamonds) and $f_{q}$
(circles), shown in Fig. 3. (b) Linear fits used to produce the data plotted
in Fig. 3. The slopes of the lines yield $\alpha_{q}$ (panel above) and
$f_{q}$ (below). For clarity we only show data for integer values of $q$ from
$-2$ (red) to $5$ (grey), with data for $q=0$ highlighted with a full symbol.
Figure 4: Singularity spectrum for $L=22$ and $N_{\text{S}}=230$, sampled for
values of $q$ from $-2$ to $5$ in steps of $1/4$ (increasing from right to
left). The ensemble contains $500$ realisations. Green circles indicate the
average over the ensemble of states near the Fermi energy
$\varepsilon_{\text{F}}=0$, while red squares indicate the ensemble over
states closest to $-0.320\text{\,}\mathrm{e}\mathrm{V}$. The corresponding
symmetrised spectra are indicated with the same colours by the error bars
only.

##  5 Validity of the scaling assumption

The question that arises when dealing with finite systems is whether the wave
function is still a multifractal. Formally speaking, the wave function is a
true multifractal only at the critical point. For a finite $L$, however, the
effective critical point shifts away as $L^{-1/\nu}$ from its thermodynamic
limit [11]. Luckily this is not a problem, since, as shown by [46], states on
the two sides of the transition still show multifractal features
characteristic of a critical wave function.

Now that we can construct a multifractal measure from the wave functions away
from the critical point, we can actually check the most important assumption
we have taken so far, namely that a localisation-delocalisation transition
occurs in our model. The histogram distribution $N_{\lambda}(\alpha)$ of the
measure $\alpha$ depends, at the critical point, only on the coarse-graining
$\lambda=l/L$, rather than separately on the system size $L$ and the box size
$l$. At the critical point, then, $N_{\lambda}(\alpha)$ has the same shape for
any $L$, provided that the wave functions are coarse-grained with the matching
$l$ box size. The dependence of $N_{\lambda}(\alpha)$ on $L$ gradually
reappears away from the critical point, where in the strong (weak) disorder
regime, larger systems become more localized (delocalized). This is shown in
Fig. 5, where we plot the ensemble histogram
$\operatorname{PDF}(\alpha)=N_{\lambda}(\alpha)\lambda^{d}/\mathcal{N}$ at
$\lambda=1/2$ for three values of $n$: the lowest in the localized regime, the
intermediate close to the critical point and the highest in the delocalised
regime. The ensemble $\operatorname{PDF}(\alpha)$ is built by filling one
single histogram with the $\mathcal{N}$ available wave functions.

Figure 5: Ensemble PDF of the multifractal measure $\alpha$ at coarse-graining
$\lambda=1/2$ and energy
$\varepsilon-\varepsilon_{\text{F}}=$-0.249\text{\,}\mathrm{e}\mathrm{V}$$, as
a function of the concentration $n$ (in units of
${10}^{20}\text{\,}\mathrm{c}\mathrm{m}^{-3}$), for two system sizes
$L^{3}=4096$ (blue dots) and $10648$ (red crosses). For clarity we show the
histogram for three concentrations: before the transition
($n=$4.6\text{\times}{10}^{20}\text{\,}\mathrm{c}\mathrm{m}^{-3}$$), near the
critical point
($6.8\text{\times}{10}^{20}\text{\,}\mathrm{c}\mathrm{m}^{-3}$), and after
($8.8\text{\times}{10}^{20}\text{\,}\mathrm{c}\mathrm{m}^{-3}$). The critical
point
($n_{\text{c}}=$6.7\text{\times}{10}^{20}\text{\,}\mathrm{c}\mathrm{m}^{-3}$$)
is indicated by a black dashed line [1]. On the bottom plane we show the
position of the average $\alpha_{0}$ also for the intermediate concentrations,
connected by lines to guide the eye. We use again blue dots with a solid line
for $L^{3}=$4096$$ and red crosses with a dashed line for $10648$. Reproduced
from the Supplemental Materials to [1].

Because of the limited spread in system sizes and the large common coarse-
graining, the difference in the PDFs outside the critical point is not very
well pronounced in Fig. 5. An alternative check consists in fixing the system
size and study how the PDF renormalises with the coarse graining [58]. The
“scaling” variable $\xi/L$ [59] scales like $\lambda^{-1}$, which implies
that, with increasing $\lambda$, $\xi/L$ becomes smaller. Physically this
means that, upon coarse-graining, localised (delocalised) states become more
localised (delocalised), or, equivalently, that the renormalisation flow
rescales the disorder away from its critical value,if a phase transition, and
hence a critical point, exists [2]. We verify this in Fig. 6 where the PDF’s
move in opposite directions upon increasing the box size $l$ in a system of
$L^{3}=$4096$$ atoms.

Figure 6: Ensemble PDF of the multifractal measure $\alpha$ at energy
$\varepsilon-\varepsilon_{\text{F}}=$-0.249\text{\,}\mathrm{e}\mathrm{V}$$ for
$L^{3}=4096$ atoms. In blue we show the PDF for $N_{\text{S}}=20$, in orange
for $N_{\text{S}}=100$. Filled symbols indicate a coarse-graining with box
size $l=1$, while empty symbols connected by a dotted line indicate $l=4$. The
thin vertical lines indicate the average values of $\alpha_{0}$ for each case.

##  6 Conclusions

After reviewing the foundations of multifractal analysis in the study of a
disordered system at criticality, we have shown the multifractal nature,
captured by the $f(\alpha)$ spectrum, of the Kohn-Sham wave function. This
result is in line with previous experimental [60], theoretical [61], and
numerical [62, 63, 64] studies, where the critical fluctuations of the wave
function at criticality are expected to survive in the presence of the Coulomb
interaction.

We have shown that multifractality persists also in the wave functions near
the critical concentrations computed in the effective tight-binding model
presented in Refs. [1, 2]. The multifractal spectrum follows the symmetry
relations derived from field-theoretical models [49] and the ensemble average
over hundreds of realizations results in negative fractal dimensions ascribed
to rare events. These are known features of (finite) multifractal critical
states that have been studied also in the non-interacting 3D Anderson model
[31]. A difference with this case, however, lies in the nearly-parabolic shape
of the $f(\alpha)$ spectrum observed here. This behaviour, referred to as weak
multifractality, is a common trait with the Anderson transition in
$2+\epsilon$ deminsions, with $\epsilon\ll 1$, and with the power-law random
banded matrix model with $b\gg 1$ [6].

Finally, we compute the $f(\alpha)$ spectrum near the two mobility edges,
respectively, far from and at the Fermi energy. In the latter case we notice
that the spectrum is narrower, indicating that the critical wave function at
this energy is, on average, more delocalized than deeper in the band. In Refs.
[1, 2] we propose the origin of this quasi-metallic behaviour to be the
hybridization of the impurity states near the Fermi energy with states from
the conduction band. The presence of this second extended band might alter the
physics of the metal-insulator transition close to the Fermi energy, leading
to a varying critical exponent across the impurity band.

A similar phenomenon has very recently been reported in a Hartree-Fock study
of the Anderson transition in the presence of tunable Coulomb interaction
[64]. There the authors show that a second mobility edge appears near the
Fermi energy, where the Coulomb gap forms. At this energy, the critical state
shows a narrower multifractal spectrum compared to the mobility edge at higher
energy. While in this model quasi-metallic behaviour appears where the Coulomb
gap forms in the centre of the band of the Anderson model, it is present in
our model when the impurity band forms and then merges with the conduction
band.

To conclude, the numerical analysis of the Anderson transition in the presence
of interactions, and, in particular, in real materials, is still very
challenging. Progress relies on a large amount of resources needed for the
simulations, with recent first-principles or self-consistent calculations
reaching system sizes of the order of $10^{2}$–$10^{4}$ sites, and of hundreds
of realisations [1, 64, 65] — while studies of the non-interacting Anderson
model can nowadays easily achieve more then $10^{6}$ sites and $10^{4}$
samples [52, 58]. Still, ab initio studies of real materials have observed new
phenomena and traced these back to the electronic interaction. Future
investigations with increasing system and sample sizes will undoubtedly
clarify the universal properties of the transitions and critical points.

## Acknowledgements

RAR thanks Presidency University, Kolkata, where part of this work was
written, for their warm hospitality. We thank EPSRC for support via the ARCHER
RAP project e420 and the MidPlus Regional HPC Centre (EP/K000128/1). UK
research data statement: data is available at
http://wrap.warwick.ac.uk/id/eprint/92910.

## References

  * [1] E. G. Carnio, N. D. M. Hine, and R. A. Römer, Physical Review B 99, 081201(R) (2019).
  * [2] E. G. Carnio, Ph.D. thesis, University of Warwick, 2018.
  * [3] P. W. Anderson, Physical Review 109, 1492 (1958).
  * [4] B. Kramer and A. MacKinnon, Reports on Progress in Physics 56, 1469 (1993).
  * [5] D. Belitz and T. R. Kirkpatrick, Reviews of Modern Physics 66, 261 (1994).
  * [6] F. Evers and A. D. Mirlin, Reviews of Modern Physics 80, 1355 (2008).
  * [7] G. A. Thomas, M. Paalanen, and T. F. Rosenbaum, Physical Review B 27, 3897 (1983).
  * [8] G. A. Thomas, Philosophical Magazine Part B 52, 479 (1985).
  * [9] H. Stupp, M. Hornung, M. Lakner, O. Madel, and H. v. Löhneysen, Physical Review Letters 71, 2634 (1993).
  * [10] S. Waffenschmidt, C. Pfleiderer, and H. v. Löhneysen, Physical Review Letters 83, 3005 (1999).
  * [11] J. L. Cardy, Scaling and renormalization in statistical physics (Cambridge University Press, Cambridge, 1996), p. 238.
  * [12] M. Lopez, J.-F. Clément, P. Szriftgiser, J. C. Garreau, and D. Delande, Physical Review Letters 108, 095701 (2012).
  * [13] D. S. Wiersma, P. Bartolini, A. Lagendijk, and R. Righini, Nature 390, 671 (1997).
  * [14] F. Scheffold, R. Lenke, R. Tweer, and G. Maret, Nature 398, 206 (1999).
  * [15] M. Störzer, P. Gross, C. M. Aegerter, and G. Maret, Physical Review Letters 96, 063904 (2006).
  * [16] T. Sperling, L. Schertel, M. Ackermann, G. J. Aubry, C. M. Aegerter, and G. Maret, New Journal of Physics 18, 013039 (2016).
  * [17] S. E. Skipetrov and J. H. Page, New Journal of Physics 18, 021001 (2016).
  * [18] T. Schwartz, G. Bartal, S. Fishman, and M. Segev, Nature 446, 52 (2007).
  * [19] H. Hu, A. Strybulevych, J. H. Page, S. E. Skipetrov, and B. A. van Tiggelen, Nature Physics 4, 945 (2008).
  * [20] S. Faez, A. Strybulevych, J. H. Page, A. Lagendijk, and B. A. van Tiggelen, Physical Review Letters 103, 155703 (2009).
  * [21] A. Aubry, L. A. Cobus, S. E. Skipetrov, B. A. Van Tiggelen, A. Derode, and J. H. Page, Physical Review Letters 112, 1 (2014).
  * [22] W. K. Hildebrand, A. Strybulevych, S. E. Skipetrov, B. A. Van Tiggelen, and J. H. Page, Physical Review Letters 112, 1 (2014).
  * [23] L. Cobus, S. Skipetrov, A. Aubry, B. van Tiggelen, A. Derode, and J. Page, Physical Review Letters 116, 193901 (2016).
  * [24] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, D. Clément, L. Sanchez-Palencia, P. Bouyer, and A. Aspect, Nature 453, 891 (2008).
  * [25] G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort, M. Zaccanti, G. Modugno, M. Modugno, and M. Inguscio, Nature 453, 895 (2008).
  * [26] D. Clément, A. F. Varón, J. A. Retter, L. Sanchez-Palencia, A. Aspect, and P. Bouyer, New Journal of Physics 8, (2006).
  * [27] S. S. Kondov, W. R. McGehee, J. J. Zirbel, and B. DeMarco, Science 334, 66 (2011).
  * [28] F. Jendrzejewski, A. Bernard, K. Müller, P. Cheinet, V. Josse, M. Piraud, L. Pezzé, L. Sanchez-Palencia, A. Aspect, and P. Bouyer, Nature Physics 8, 398 (2012).
  * [29] G. Semeghini, M. Landini, P. Castilho, S. Roy, G. Spagnolli, A. Trenkwalder, M. Fattori, M. Inguscio, and G. Modugno, Nature Physics 11, 554 (2015).
  * [30] K. Slevin and T. Ohtsuki, Physical Review Letters 82, 382 (1999).
  * [31] A. Rodriguez, L. J. Vasquez, K. Slevin, and R. A. Römer, Physical Review Letters 105, 046403 (2010).
  * [32] K. M. Itoh, M. Watanabe, Y. Ootuka, E. E. Haller, and T. Ohtsuki, Journal of the Physical Society of Japan 73, 173 (2004).
  * [33] J. Neugebauer and T. Hickel, Wiley Interdisciplinary Reviews: Computational Molecular Science 3, 438 (2013).
  * [34] A. Jain, Y. Shin, and K. A. Persson, Nature Reviews Materials 1, 15004 (2016).
  * [35] M. T. Winkler, D. Recht, M.-J. Sher, A. J. Said, E. Mazur, and M. J. Aziz, Physical Review Letters 106, 178701 (2011).
  * [36] C.-K. Skylaris, P. D. Haynes, A. A. Mostofi, and M. C. Payne, The Journal of Chemical Physics 122, 084119 (2005).
  * [37] H. Aoki, Physica A: Statistical Mechanics and its Applications 114, 538 (1982).
  * [38] H. Aoki, Journal of Physics C: Solid State Physics 16, L205 (1983).
  * [39] B. B. Mandelbrot, Journal of Statistical Physics 34, 895 (1984).
  * [40] C. Castellani and L. Peliti, Journal of Physics A: Mathematical and General 19, L429 (1986).
  * [41] F. Wegner, Zeitschrift für Physik B Condensed Matter and Quanta 36, 209 (1980).
  * [42] B. B. Mandelbrot, Possible refinement of the lognormal hypothesis concerning the distribution of energy dissipation in intermittent turbulence (Springer, Berlin, Heidelberg, 1972), pp. 333–351.
  * [43] B. B. Mandelbrot, Journal of Fluid Mechanics 62, 331 (1974).
  * [44] L. Pietronero, The Superworld III (Springer US, Boston, MA, 1990), pp. 175–183.
  * [45] O. Malcai, D. A. Lidar, O. Biham, and D. Avnir, Physical Review E 56, 2817 (1997).
  * [46] E. Cuevas and V. E. Kravtsov, Physical Review B 76, 235119 (2007).
  * [47] U. Frisch and G. Parisi, in Turbulence and predictability in geophysical fluid dynamics and climate dynamics, Proceedings of the International School of Physics ‘E. Fermi’ (Varenna, Italy, 1983), edited by M. Ghil, R. Benzi, and G. Parisi (North-Holland Publishing Company, Amsterdam, 1985), pp. 84–87.
  * [48] M. Janssen, International Journal of Modern Physics B 08, 943 (1994).
  * [49] A. D. Mirlin, Y. V. Fyodorov, A. Mildenberger, and F. Evers, Physical Review Letters 97, 1 (2006).
  * [50] A. Rodriguez, L. J. Vasquez, and R. A. Römer, Physical Review B 78, 195107 (2008).
  * [51] L. J. Vasquez, A. Rodriguez, and R. A. Römer, Physical Review B 78, 195106 (2008).
  * [52] A. Rodriguez, L. J. Vasquez, K. Slevin, and R. A. Römer, Physical Review B 84, 134209 (2011).
  * [53] A. Rodriguez, L. J. Vasquez, and R. A. Römer, Physical Review Letters 102, 106406 (2009).
  * [54] A. B. Chhabra, C. Meneveau, R. V. Jensen, and K. R. Sreenivasan, Physical Review A: Atomic, Molecular, and Optical Physics 40, 5284 (1989).
  * [55] S. Thiem and M. Schreiber, The European Physical Journal B 86, 48 (2013).
  * [56] A. B. Chhabra and K. R. Sreenivasan, Physical Review A 43, 1114 (1991).
  * [57] A. D. Mirlin and F. Evers, Physical Review B 62, 7920 (2000).
  * [58] J. Lindinger and A. Rodríguez, Physical Review B 96, 134202 (2017).
  * [59] A. MacKinnon and B. Kramer, Physical Review Letters 47, 1546 (1981).
  * [60] A. Richardella, P. Roushan, S. Mack, B. Zhou, D. A. Huse, D. D. Awschalom, and A. Yazdani, Science 327, 665 (2010).
  * [61] I. S. Burmistrov, I. V. Gornyi, and A. D. Mirlin, Physical Review Letters 111, 066601 (2013).
  * [62] M. Amini, V. E. Kravtsov, and M. Müller, New Journal of Physics 16, 015022 (2014).
  * [63] Y. Harashima and K. Slevin, Physical Review B 89, 205108 (2014).
  * [64] H.-J. Lee and K.-S. Kim, Physical Review B 97, 155105 (2018).
  * [65] Y. Zhang, R. Nelson, K.-M. Tam, W. Ku, U. Yu, N. S. Vidhyadhiraja, H. Terletska, J. Moreno, M. Jarrell, and T. Berlijn, Physical Review B 98, 174204 (2018).


# Artificial to Spiking Neural Networks Conversion for Scientific Machine
Learning ††thanks: Citation: Authors. Title. Pages…. DOI:000000/11111.

Qian Zhang, Chenxi Wu, Adar Kahana, George Em Karniadakis
Division of Applied Mathematics
Brown University
Providence
{qian_zhang1, chenxi_wu, adar_kahana<EMAIL_ADDRESS>
&Youngeun Kim, Yuhang Li, Priyadarshini Panda
Department of Electrical Engineering
Yale University
New Haven
{youngeun.kim, yuhang.li<EMAIL_ADDRESS>

###### Abstract

We introduce a method to convert Physics-Informed Neural Networks (PINNs),
commonly used in scientific machine learning, to Spiking Neural Networks
(SNNs), which are expected to have higher energy efficiency compared to
traditional Artificial Neural Networks (ANNs). We first extend the calibration
technique of SNNs to arbitrary activation functions beyond ReLU, making it
more versatile, and we prove a theorem that ensures the effectiveness of the
calibration. We successfully convert PINNs to SNNs, enabling computational
efficiency for diverse regression tasks in solving multiple differential
equations, including the unsteady Navier-Stokes equations. We demonstrate
great gains in terms of overall efficiency, including Separable PINNs
(SPINNs), which accelerate the training process. Overall, this is the first
work of this kind and the proposed method achieves relatively good accuracy
with low spike rates.

_K_ eywords Spiking Neural Networks $\cdot$ Conversion $\cdot$ Nonlinear
activation

##  1 Introduction

The use of machine learning techniques in the scientific community has been
spreading widely, reaching many fields such as physics [1, 2, 3, 4], chemistry
[5, 6], biology [7, 8, 9], geophysics [10, 11], epidemiology[12, 13] and many
more. The advances in computation capabilities have enabled many researchers
to reformulate diverse problems as data-driven problems, by combining prior
knowledge of the problem with fitting a model for the available data. A
prominent drawback of Scientific Machine Learning (SciML) techniques is that
they are usually expensive in terms of computational cost. They require either
knowledge of the governing equations that determine the process (approximating
them is a costly procedure), or a large amount of data to fit (expensive as
well). The SciML community is striving for a more efficient method for
training and inferring neural networks. Neuromorphic chips are one edge
computing component that SciML applications could benefit from. In this work,
we explore methods for enabling this.

An important breakthrough in the field of SciML was the invention of the
Physics-Informed Neural Networks (PINNs) [14, 15]. PINNs incorporate the
knowledge of the physical experiment and governing equations into the network
training step, making it a hybrid (physics and data) training. PINNs and its
extensions [16, 17, 18, 19] have achieved great success in many fields and
applications [20, 1, 21, 22, 23, 24, 7, 25, 26] . A disadvantage of PINNs is
that like other deep neural networks, they are prone to long training times.
In addition, when changing the problem conditions (initial conditions,
boundary conditions, domain properties, etc.), the PINN has to be trained from
scratch. Therefore, for real-time applications, a more efficient solution is
sought after. For training a PINN, one usually uses smooth activation
functions (such as the Tanh or Sine activation functions), where in most ANNs
the ReLU activation is dominant. Using smooth activation function in a SNN is
a new challenge we address in this paper with theoretical justification.

Spiking Neural Networks (SNNs) have been gaining traction in the machine
learning community for the past few years. The main reason is their expected
efficiency [27, 28, 29, 30], in terms of energy consumption, compared to their
Artificial Neural Network (ANN) counterparts that are commonly used for many
applications. In addition, the advances in neuromorphic hardware (such as
Intel’s Loihi 2 chip [31, 32]), call for innovative algorithms and software
that can utilize the chips and produce lighter and faster machine learning
models. However, developing an SNN is a challenging task, especially for
regression [33, 34]. Studies have been conducted for translating components
from the popular ANNs into a spiking framework [35], but many components are
not yet available in the spiking regime. In this paper we focus on that
specific aspect.

There are three popular approaches for training SNNs. The first involves using
mathematical formulations of the components of the brain, such as the membrane
[36, 37, 38], the synapse [39], etc. In this case, one uses a Hebbian learning
rule [40] to find the weights of the synapses (the trainable parameters) using
forward propagation (without backward propagation [41, 42]). The second method
involves building surrogate models for the elements in the SNN that block the
back-propagation, such as the non-differentiable activation functions used in
SNNs. The third method, which is discussed in this paper, addresses converting
a trained ANN into a SNN.

The main contributions of this paper are as follows:

  1. 1.

We propose a method to convert PINNs, a type of neural network commonly used
for regression tasks, to Spiking Neural Networks (SNNs). The conversion allows
for utilizing the advantages of SNNs in the inference stage, such as
computational efficiency, in regression tasks.

  2. 2.

We extend the calibration techniques used in previous studies to arbitrary
activation functions, which significantly increases the applicability of the
conversion method. Furthermore, we provide a convergence theorem to guarantee
the effectiveness of the calibration.

  3. 3.

We apply the conversion to separable PINNs (SPINNs), which accelerates the
training process of the PINNs.

Overall, the proposed method extends the application of SNNs in regression
tasks and provides a systematic and efficient approach to convert existing
neural networks for diverse regression tasks to SNNs.

##  2 Related Work

#### Physics-informed neural networks (PINNs):

An innovative framework that combines neural networks with physical laws to
learn complex physical phenomena. In PINNs, the physical equations are
integrated into the loss function, which allows the network to learn from both
the given data and the underlying physics. This approach significantly
improves the network’s ability to handle incomplete or noisy data and performs
well with limited training data. PINNs have been successfully applied to a
range of problems in fluid dynamics, solid mechanics, and more [22, 1, 21, 23,
7, 24].

#### Separable PINNs (SPINNs):

Cho et al. [43] proposed a novel neural network architecture called SPINNs,
which aims to reduce the computational demands of PINNs and alleviate the
curse of dimensionality. Unlike vanilla PINNs that use point-wise processing,
SPINN works on a per-axis basis, thereby reducing the number of required
network forward passes. SPINN utilizes factorized coordinates and separated
sub-networks, where each sub-network takes an independent one-dimensional
coordinate as input, and the final output is generated through an outer
product and element-wise summation. Because SPINN eliminates the need to query
every multidimensional coordinate input pair, it is less affected by the
exponential growth of computational and memory costs associated with grid
resolution in standard PINNs. Furthermore, SPINNs operate on a per-axis basis,
which allows for parallelization with multiple GPUs.

#### Spiking Neural Networks (SNNs):

A type of Artificial Neural Network (ANN), that differs in the implementation
of the core components. The purpose is to create more biologically-plausible
training and inference procedures. Unlike traditional ANNs, which process
information through numerical values, SNNs process information through spikes,
which occur in response to stimulation (much like the human brain). SNNs are
becoming increasingly popular as they can mimic the temporal nature of
biological neurons. Additionally, SNNs are computationally efficient and have
the potential for efficient hardware implementation, making them well-suited
for real-time applications. Combining SNN implementation with edge computing,
training could be significantly faster, and inference as well [27, 28, 29].
Recent results have shown that SNNs can achieve high accuracy on image
classification tasks, with up to 99.44% on the MNIST dataset [44] and up to
79.21% on ImageNet [45].

#### SNN conversion

A technique to transform a trained ANN into a SNN. SNN conversion usually
involves mapping the weights and activation functions of the ANN to the
synaptic strengths and spike rates of SNNs. It is considered the most
efficient way to train deep SNNs, as it avoids the challenges of direct SNN
training, such as gradient estimation and spike generation [46, 47, 48]. The
algorithm of SNN conversion can be divided into two steps: offline conversion
and online inference. In the offline conversion step, the trained ANN model is
converted into an equivalent SNN model by adjusting the network parameters. In
the online inference step, the converted SNN model is used for the inference.
In the online step, the SNN is intended to be deployed on neuromorphic
hardware, to unlock its full potential and energy efficiency.

##  3 Method

A complementary process to the conversion technique. The SNN calibration [49]
is a method that minimizes the loss of accuracy and efficiency when converting
an ANN into a SNN. SNN calibration leverages the knowledge of the pre-trained
ANN and corrects the conversion error layer by layer. The algorithm consists
of two steps: replacing the ReLU activation function with a spike activation
function, and applying the calibration to adjust only the biases (light) or
weights and biases (advanced) of each layer. The calibration method is based
on the theoretical analysis of the conversion error and its propagation
through the network. SNN calibration can achieve comparable or even better
performance than the original ANN on various datasets and network
architectures. This paper is a generalization of this work. We propose an
extension to the SNN conversion, which is apropriate for regression tasks.

###  3.1 SNN Conversion setup

We consider a dataset $D=(X,Y)$, and an ANN $\mathcal{A}$ with $n$ hidden
layers, trained to fit it. Let $\mathbf{x}^{(n)}=\mathcal{A}(X)$ be the output
of $\mathcal{A}$. The goal of SNN conversion is to find an SNN $\mathcal{S}$,
whose (averaged) output is $\bar{\mathbf{s}}^{(n)}=\mathcal{S}(X)$, such that
$\mathbf{x}^{(n)}$ is close to $\bar{\mathbf{s}}^{(n)}$. In other words, we
want to minimize the norm of the error
$\mathbf{e}^{(n)}\stackrel{{\scriptstyle
d}}{{=}}\mathbf{x}^{(n)}-\bar{\mathbf{s}}^{(n)}$ for a given $\mathcal{A}$ and
$D$. The same network structure as $\mathcal{S}$ is used for $\mathcal{A}$ and
the activation function of $\mathcal{A}$ is replaced with IF. Then we can
analyze the factors that influence the total conversion error
$\mathbf{e}^{(n)}$.

###  3.2 SNN Conversion with Calibration

In this section, we will briefly explain the SNN conversion with calibration
proposed in [49]. Consider an MLP model with $n$ layers, the first $n-1$
layers use ReLU activation and the last layer has no activation. We denote
$\mathbf{W}^{(l)}$ as the weights and bias for layer $l$.

The naive way of SNN conversion is simply replacing the ReLU activation layers
with IF activation layers. In this case, we define $\mathbf{x}^{(l)}$ as the
output of layer $l$ recursively, which is
$\mathbf{x}^{(l)}=\operatorname{ReLU}(\mathbf{W}^{(l)}\mathbf{x}^{(l-1)})$ and
$\mathbf{x}^{(0)}$ is the input. Similarly, we define
$\bar{\mathbf{s}}^{(l)}$, the output of the layer $l$ of converted SNN, as
$\bar{\mathbf{s}}^{(l)}=\operatorname{IF}(\mathbf{W}^{(l)}\bar{\mathbf{s}}^{(l-1)})$
and $\bar{\mathbf{s}}^{(0)}=\mathbf{x}^{(0)}$ is the input. In fact, we can
compute the expected output spikes as
$\bar{\mathbf{s}}^{(l)}=\operatorname{ClipFloor}(\mathbf{W}^{(l)}\bar{\mathbf{s}}^{(l-1)})$.
The $\operatorname{ClipFloor}$ function is an approximation to
$\operatorname{ReLU}$ and is illustrated in Figure 1.

Figure 1: Illustration of $\operatorname{ClipFloor}$ function and its
relationship to $\operatorname{ReLU}$ (left) and $\tanh$ (right).

Then we can define the conversion error of layer $l$ as
$\mathbf{e}^{(l)}=\mathbf{x}^{(l)}-\bar{\mathbf{s}}^{(l)}$ and decompose it as

$\displaystyle\mathbf{e}^{(l)}$
$\displaystyle=\mathbf{x}^{(l)}-\bar{\mathbf{s}}^{(l)}$ (1)
$\displaystyle=\operatorname{ReLU}(\mathbf{W}^{(l)}\mathbf{x}^{(l-1)})-\operatorname{ClipFloor}(\mathbf{W}^{(l)}\bar{\mathbf{s}}^{(l-1)})$
$\displaystyle=\operatorname{ReLU}(\mathbf{W}^{(l)}\mathbf{x}^{(l-1)})-\operatorname{ReLU}(\mathbf{W}^{(l)}\bar{\mathbf{s}}^{(l-1)})+\operatorname{ReLU}(\mathbf{W}^{(l)}\bar{\mathbf{s}}^{(l-1)})-\operatorname{ClipFloor}(\mathbf{W}^{(l)}\bar{\mathbf{s}}^{(l-1)})$
$\displaystyle=\mathbf{e}^{(l)}_{r}+\mathbf{e}^{(l)}_{c}$

Here $\mathbf{e}^{(l)}_{r}$ represents the error caused by approximating the
continuous input by the spiking input, and $\mathbf{e}^{(l)}_{c}$ represents
the local conversion error caused by changing the smooth activation function
to the spiking activation (IF). A major result in [49] is that we can bound
the conversion error $\mathbf{e}^{(n)}$, with the weighted sum of the local
conversion errors $\mathbf{e}^{(l)}_{c}$. This allows us to minimize the
conversion error via optimizing each $\mathbf{e}^{(l)}_{c}$.

###  3.3 Results for general activation functions

In fact, the results can be generalized to activation functions other than
$\operatorname{ReLU}$ by similar techniques in [49]. Since the generalized
activation functions may have negative values, we introduce the idea of
negative threshold, a concept in SNNs that allows neurons to fire both
positive and negative spikes, depending on their membrane potential [50]. A
positive spike occurs when the membrane potential exceeds the positive
threshold, and a negative spike occurs when it is below the negative
threshold. This mimics the biological behavior of neurons that do not fire a
spike when the membrane potential does not reach the threshold. Negative
threshold can be applied to different types of SNNs and learning functions
depending on the problem domain and the data characteristics. It is very
helpful when the dataset contains negative values.

To formulate the conversion with calibration for generalized activations, we
consider the conversion error decomposition

$\displaystyle\mathbf{e}^{(l)}$
$\displaystyle=\mathbf{x}^{(l)}-\bar{\mathbf{s}}^{(l)}$ (2)
$\displaystyle=f(\mathbf{W}^{(l)}\mathbf{x}^{(l-1)})-\operatorname{ClipFloor}(\mathbf{W}^{(l)}\bar{\mathbf{s}}^{(l-1)})$
$\displaystyle=f(\mathbf{W}^{(l)}\mathbf{x}^{(l-1)})-f(\mathbf{W}^{(l)}\bar{\mathbf{s}}^{(l-1)})+f(\mathbf{W}^{(l)}\bar{\mathbf{s}}^{(l-1)})-\operatorname{ClipFloor}(\mathbf{W}^{(l)}\bar{\mathbf{s}}^{(l-1)})$
$\displaystyle=\mathbf{e}^{(l)}_{r}+\mathbf{e}^{(l)}_{c}$

where $\mathbf{e}^{(l)}_{r}$ and $\mathbf{e}^{(l)}_{c}$ have the same meaning
as in Eq 1. Then we can also use the weighted sum of the local conversion
errors $\mathbf{e}^{(l)}_{c}$ to bound the total conversion error
$\mathbf{e}^{(n)}$ by the following theorem.

###### Theorem 1

For any activation function whose function values and first-order derivatives
can be uniformly approximated by piecewise linear functions and up to second-
order derivatives are bounded, then the conversion error in the final network
output space can be bounded by a weighted sum of local conversion errors,
given by

$\mathbf{e}^{(n),\top}\mathbf{H}^{(n)}e^{(n)}\leq\sum_{l=1}^{n}2^{n-l+1}\mathbf{e}_{c}^{(l),\top}(\mathbf{H}^{(l)}+K_{L}^{(l)}\sqrt{L}\mathbf{I})\mathbf{e}_{c}^{(l)}$
(3)

where $L$ is the training loss.

We present the detailed proof in the Appendix. The main technique is to find
piecewise linear functions that approximate the smooth activation function.

Therefore, the total conversion error is bounded from above by the local
conversion errors. To achieve more accurate conversion performance, we can
minimize the total conversion error by minimizing the local conversion error
layerwise, which can be easily implemented as in [49]. We enlarge the last
layer threshold to preserve the maximum value of the output. The details are
discussed in the Appendix.

##  4 Results

###  4.1 Function regression

We first present an example of function regression: using neural networks to
approximate the $\sin$ function. For the training dataset, the input is the
uniform mesg points on $[-\pi.\pi]$, and the output is the values of $\sin$ on
these points. The ANN model has two intermediate layers, each of them has 40
neurons. The activation function is $\tanh$ except for the last layer; there
is no activation function for the last layer. The network is trained with the
Adam optimizer until the training error is less then $1^{-7}$. Then, we
convert the ANN to SNN with advanced calibration with different numbers of
time steps. The results are shown in Figure 2.

(a) T=8

(b) T=32

(c) T=128

Figure 2: Results of converting an ANN, which is trained to approximate
$\sin$, to SNN. The number of time steps is $T=8,32,128$ from left to right.
The output of SNN is close to the ground truth, which is the output of ANN.
With increasing T, the conversion error becomes smaller. When $T=128$, the SNN
output curve is almost smooth and similar to the ground truth.

To further investigate the impact of $T$ on the conversion error, we train
networks with different intermediate layer numbers $L=2,3,4$ and neuron
numbers per layer $N=20,40,60,80,100$. All the other setups are the same. We
obtain the results shown in Figure 3.

(a) L=2

(b) L=3

(c) L=4

Figure 3: The number of layers are $L=2,3,4$ from left to right. The error
here refers to $\mathbf{e}^{(n)}$. When $T\leq 64$, the conversion error
follows $\mathbf{e}^{(n)}\approx 1/T$, and $\mathbf{e}^{(n)}$ will get stable
after $T$ is already very large.

We can find that the conversion error decreases with larger $T$ (more
specifically, conversion error $\sim 1/T$ where ) when $T<32$, but it becomes
stable or larger after $T\geq 32$. And for neural network with fixed depth,
larger layer width usually leads to smaller conversion error. However, for
fixed layer width, deeper neural networks will not bring significantly better
conversion performance.

To validate the Theorem 1, we need to compute
$\mathbf{e}^{(l),\top}\mathbf{H}^{(l)}\mathbf{e}^{(l)}$ for $l=1,2,\ldots,n$.
Since $\mathbf{H}^{(l)}$ are intractable, we can only use the identity matrix
to replace it and obtain some qualitative results. That is, compute
$\norm{\mathbf{e}^{(l)}}^{2}$ instead of
$\mathbf{e}^{(l),\top}\mathbf{H}^{(l)}\mathbf{e}^{(l)}$. So the computed RHS
is $\sum_{l=1}^{n}2^{n-l+1}\mathbf{e}^{(l),\top}\mathbf{e}^{(l)}$. We train
the ANN with layer numbers $L=2,3,4$ with $100$ neurons per layer. All the
other setups are the same. The results are shown in Figure 4.

(a) L=2

(b) L=3

(c) L=4

Figure 4: The number of layers are $L=2,3,4$ from left to right. The error
here refers to $\mathbf{e}^{(n)}$. RHS is defined as before. We can find that
the total conversion error is smaller than the computed RHS and decrease with
the RHS as well, which is an emprical validation of the Theorem 1.

Although the RHS term is not exactly the same as in 1, the trend agrees with
our statement, that conversion error decreases as the RHS term does.

###  4.2 PINNs

To show the power of SNN conversion in regression tasks, we train PINNs, a
MLP-based neural network which can solve PDEs, and convert them to SNNs. We
present results for the Poisson equation, the diffusion-reaction equation, the
wave equation, the Burgers equation, and the Navier-Stokes equations.

#### Poisson equation

The Poisson equation is often used to describe the potential field in
electromagnetic theory. Here we solve the following boundary value problem of
Poisson equation

$\displaystyle-\Delta u(x)$ $\displaystyle=1,\quad
x\in\Omega=[-1,1]\times[-1,1]$ (4) $\displaystyle u(x)$ $\displaystyle=0,\quad
x\in\partial\Omega$

with PINN. The network has $3$ intermediate layers, each of which has 100
neurons. The activation function is $\tanh$ except for the last layer. The
network is trained for 50,000 epochs. Then we convert it into SNN. The results
are shown in Figure 5.

(a) Reference solution

(b) PINN result

(c) Conversion w/o calibration

(d) Conversion w/ calibration

Figure 5: Poisson equation: The results of converting a PINN solving the
poisson equation (4). Figure 5(a) is the reference solution. Figure 5(b) is
the PINN result. Figure 5(c) is the result of the SNN converted from the PINN
without using calibration. Figure 5(c) is the result of the SNN converted from
the PINN using calibration. Here the L2 error and relative error are defined
as $\norm{\mathbf{x}^{(n)}-\bar{\mathbf{s}}^{(n)}}_{2}$ and
$\norm{\mathbf{x}^{(n)}-\bar{\mathbf{s}}^{(n)}}_{2}/\norm{\mathbf{x}^{(n)}}_{2}$,
where $\mathbf{x}^{(n)}$ is the reference solution and
$\bar{\mathbf{s}}^{(n)}$ is the neural network output. $\norm{\cdot}_{2}$ is
the $l^{2}$ norm, which is the root of mean square error.

We observe that the SNN converted with calibration can achieve the same
magnitude of error as the PINN evaluated as an ANN, while the SNN converted
without using calibration can just obtain a rough shape of the solution and
has much larger error.

#### Diffusion-reaction equation

This equation models reactive transport in physical and biologicak systems.
Here we use PINN to solve the diffusion-reaction equation with the following
initial condition:

$\displaystyle u_{t}-u_{xx}$ $\displaystyle=ku^{2},\quad x\in\Omega=[-1,1]$
(5) $\displaystyle u(x,0)$ $\displaystyle=\exp(-\frac{x^{2}}{2\sigma^{2}})$

where $k=1$, $\sigma=0.25$ up to time $T=0.01$. The network has $3$
intermediate layers, each of which has 100 neurons. The activation function is
$\tanh$ except for the last layer. The network is trained for 100,000 epochs.
Then we convert it into SNN. The results are shown in 6.

(a) Reference solution

(b) PINN result

(c) Conversion w/o calibration

(d) Conversion w/ calibration

Figure 6: Reaction-diffusion equation: The results of converting a PINN
solving the nonlinear heat equation (5). Figure 6(a) is the reference
solution. Figure 6(b) is the PINN result. Figure 6(c) is the result of the SNN
converted from the PINN without using calibration. Figure 6(c) is the result
of the SNN converted from the PINN using calibration.

#### Wave equation

Here we use PINN to solve the a wave equation with the following initial and
boundary conditions:

$\displaystyle u_{tt}-u_{xx}$ $\displaystyle=0,\quad x\in\Omega=[-1,1]$ (6)
$\displaystyle u(x,0)$ $\displaystyle=\left\\{\begin{aligned} &1\quad
x\in[-0.245,0.245]\\\ &0\quad x\in[-1,-0.6]\cap[0.6,1]\\\
&\mathrm{linear}\quad\mathrm{otherwise}\end{aligned}\right.$ $\displaystyle
u(-1,t)$ $\displaystyle=u(1,t)=0$

up to time $T=0.5$. The network has $3$ intermediate layers, each of which has
100 neurons. The activation function is $\tanh$ except for the last layer. The
network is trained for 100,000 epochs. Then we convert it into SNN. The
results are shown in 7.

(a) Reference solution

(b) PINN result

(c) Conversion w/o calibration

(d) Conversion w/ calibration

Figure 7: Wave equation: The results of converting a PINN solving the wave
equation (6). Figure 7(a) is the reference solution. Figure 7(b) is the PINN
result. Figure 7(c) is the result of the SNN converted from the PINN without
using calibration. Figure 7(c) is the result of the SNN converted from the
PINN using calibration.

#### Viscous Burgers equation

The Burgers equation is a prototype PDE representing nonlinear advection-
diffusion occurs fluid mechanics. Here we solve the following problem of
viscous Burgers equation:

$\displaystyle\frac{\partial u}{\partial t}-\frac{\partial}{\partial
x}(\frac{1}{2}u^{2})$ $\displaystyle=\nu\frac{\partial^{2}u}{\partial
x^{2}},\quad(x,t)\in[0,2\pi]\times[0,4]$ (7) $\displaystyle u(x,0)$
$\displaystyle=\sin(x)$

with PINN. The network has $6$ intermediate layers, each of which has $40$
neurons. The activation function is $\tanh$ except for the last layer. The
network is trained for 100,000 epochs. Then we convert it into SNN. The
results are shown in 8.

(a) Reference solution

(b) PINN result

(c) Conversion w/o calibration

(d) Conversion w/ calibration

Figure 8: Burgers equation: The results of converting a PINN solving the
viscous Burgers equation (7). Figure 8(a) is the reference solution. Figure
8(b) is the PINN result. Figure 8(c) is the result of the SNN converted from
the PINN without using calibration. Figure 8(c) is the result of the SNN
converted from the PINN using calibration.

Here we can find that conversion without calibration does not give correct
position of steep gradient as the position is moving left (see Figure 8(c)).
But the conversion with calibration keeps the steep gradient position correct,
which is important for physics.

Due to the discontinuous nature of SNN, the conversion results are not smooth.
However, we can still apply some filters to smooth the outputs. For example,
we apply FFT to the conversion results and remove the high frequencies, and
the results are shown in Figure

(a) Conversion w/ calibration

(b) Conversion w/ calibration

Figure 9: Figure 9(a) is the original conversion result. Figure 9(b) is the
smoothed conversion result.

After the smoothing, conversion error becomes lower, hence can be adopted as
the postprocessing procedure.

###  4.3 Accelerating training with separable PINNs (SPINNs)

A current limitation of PINN-SNN conversion is the computational resources it
requires, especially for training the ANN. Herein, we highlight the
effectiveness of using the separable physics-informed neural networks to
spiking neural networks (SPINN-SNN) conversion. The SPINN-SNN conversion
pipeline enhances the speed of operations and provides great computational
efficiency compared to the direct application of PINN. We implement the SPINN-
SNN conversion to address two PDEs: a two-dimensional viscous Burgers equation
and a three-dimensional unsteady Beltrami flow. We provide a comparative
analysis of the training time taken by standard PINN and the SPINN-SNN
conversion pipeline. Experimental results reveal that the application of
SPINN-SNN conversion greatly enhances the speed of the training, particularly
when solving high dimensional PDEs. Figure 10 displays a runtime comparison
between SPINN-SNN conversion and PINN, applied to two and three-dimensional
problems. When addressing the two-dimensional Burgers equation, the SPINN-SNN
conversion process is approximately 1.7 times faster than PINN. The
superiority of SPINN-SNN conversion becomes more apparent while solving the
three-dimensional Beltrami flow problem, which is over 60 times faster than
PINN. Notably, the time necessary for the SNN calibration remains relatively
constant, regardless of the problem dimensionality. This indicates that the
benefits of the SPINN-SNN conversion pipeline become increasingly prominent
with the rise in the dimensionality of the problem.

(a) 2D Burgers equation runtime

(b) 3D Beltrami flow runtime

Figure 10: The runtimes of the SPINN-SNN conversions solving the viscous
Burgers equation (7) and Beltrami flow (8). Figure 10(a) is the runtime for
the 2D viscous Burgers equation. Figure 10(b) is the runtime for the 3D
Beltrami flow.

#### Viscous Burgers equation

In order to conduct a fair comparison between the PINN-SNN and SPINN-SNN
conversions, the setup of the viscous Burgers equation is kept consistent with
that described in Equation 7, and the same set of hyperparameters is utilized.
Figure 11 presents the results of converting a SPINN to solve the Burgers’
equation. The SPINN has individual subnetworks for each independent variable,
$x$ and $t$. Each of these subnetworks comprises three intermediate layers,
each layer containing 40 neurons, and employs the tanh activation function,
except in the last layer. As depicted in Figure 11, the conversion from SPINN
provides an accuracy level comparable to that of PINN. An SNN converted with
calibration achieves a significantly smaller error than one converted without
calibration.

(a) Reference solution

(b) Separable PINN result

(c) Conversion w/o calibration

(d) Conversion w/ calibration

Figure 11: Burgers equation with Separable PINN (SPINN): The results of
converting a Separable PINN solving the viscous Burgers equation (7). Figure
11(a) is the reference solution. Figure 11(b) is the Separable PINN result.
Figure 11(c) is the result of the SNN converted from the Separable PINN
without using calibration. Figure 11(c) is the result of the SNN converted
from the Separable PINN using calibration.

#### Beltrami Flow

The Navier-Stokes equations are fundamental in fluid mechanics as they
mathematically represent the conservation of momentum in fluid systems, and
there has been significant advancement in solving Navier-Stokes flow problems
using scientific machine learning methods [51, 52, 53]. The Navier-Stokes
equations can be presented in two forms: the velocity-pressure (VP) form and
the vorticity-velocity (VV) form. The incompressible Navier-Stokes equations,
in their VP form, are as follows:

$\displaystyle\frac{\partial{\mathbf{u}}}{\partial
t}+(\mathbf{u}\cdot\nabla)\mathbf{u}$ $\displaystyle=-\nabla
p+\frac{1}{Re}\nabla^{2}\mathbf{u}$ (8) $\displaystyle\nabla\cdot\mathbf{u}$
$\displaystyle=0$

We chose the spatial domain to be $\Omega\in[-1,1]^{2}$ and time interval to
be $\Gamma\in[0,1]$. Here, $t$ is non-dimensional time,
$\mathbf{u}(\mathbf{x},t)=[u,v]^{T}$ the non-dimensional velocity in the
$(x,y)$-directions, $p$ the non-dimensional pressure, and the Reynolds number
$Re=\frac{U_{ref}D_{ref}}{v}$ is defined by characteristic length ($D_{ref}$),
reference velocity ($U_{ref}$), and kinematic viscosity ($v$). In this
example, we simulate a three-dimensional unsteady laminar Beltrami flow where
$Re=1$. The analytical solution of the Beltrami flow is [54]:

$\displaystyle u(x,y,t)$ $\displaystyle=-\cos x\sin y\text{ }e^{-2t}$ (9)
$\displaystyle v(x,y,t)$ $\displaystyle=\sin x\cos y\text{ }e^{-2t}$
$\displaystyle p(x,y,t)$ $\displaystyle=-\frac{1}{4}(\cos 2x+\cos 2y)\text{
}e^{-4t}$

The boundary and initial conditions are extracted from this exact solution.
The PINN network comprises 4 intermediate layers, with each one containing 128
neurons. The activation function is tanh applied taoi all layers except the
final layer. The network is trained for 20,000 epochs before being converted
into an SNN. The SPINN consists of separate subnetworks for each independent
variable, u, v, and p. Every subnetwork has 2 intermediate layers, with each
layer consisting of 50 neurons. They all utilize the tanh activation function,
apart from the last layer. Figure 12 illustrates that the SPINN conversion
offers similar accuracy compared to PINN. The error for pressure is slightly
more than the velocity errors in both the $x$ and $y$ axes. Nevertheless, the
SNN, once calibrated and converted using SPINN, attains good accuracy and
notable speed improvements in comparison to PINN.

Figure 12: Beltrami flow with Separable PINN (SPINN): The results of
converting a SPINN solving the Beltrami flow problem(7). Figure 12(a,e,i) are
the reference solutions for $\mathbf{u},\mathbf{v},\mathbf{p}$ respectively.
Figure 12(b,f,j) are the SPINN results. Figure 12(c,g,k) are the PINN results.
Figure 12(d,h,i) are the results of the SNN converted from the SPINN with
calibration.

###  4.4 Firing rates of the converted SNN

To demonstrate the potential efficiency of converting the ANN to the SNN, we
computed the spiking rates of different equations. The spiking rate is defined
as the ratio of non-zero values in the output of each layer. Prior works have
suggested that SNNs with lower spiking rates will translate to energy
efficient impelmentations on neurmorphic hardware. The results are shown in
Table 1. We observe $<0.5$ spiking rate in most cases demonstrating that SNNs
only expend $<50\%$ of their network computations.

Table 1: The spiking rate of SNNs for different equations. Equations | Number of parameters | Spiking rate
---|---|---
Poisson | 20601 | 0.3727
Diffusion-reaction | 20601 | 0.2879
Wave | 20601 | 0.5721
Burgers | 8361 | 0.7253
Burgers (Separable PINN) | 15500 | 0.2831
N-S (Beltrami flow, Separable PINN) | 30900 | 0.1754
Table 2: The $L^{2}$ and relative $L^{2}$ error of converted SNN for different equations. Equations |  $L^{2}$ error | Relative $L^{2}$ error
---|---|---
Poisson | $7.8508\times 10^{-3}$ | $4.7095\times 10^{-2}$
Diffusion-reaction | $2.8766\times 10^{-2}$ | $6.3267\times 10^{-2}$
Wave | $3.8965\times 10^{-2}$ | $6.5308\times 10^{-2}$
Burgers | $3.9884\times 10^{-2}$ | $6.9781\times 10^{-2}$
Burgers (Separable PINN) | $6.6495\times 10^{-2}$ | $1.1634\times 10^{-1}$
N-S (Beltrami flow, Separable PINN) | $8.2512\times 10^{-3}$ | $4.3273\times 10^{-2}$

##  5 Conclusion

We have successfully extended the SNN calibration method proposed in [49] to a
more generalized class of activation functions beyond ReLU. The original proof
relied on the specific property of ReLU’s second-order derivative being zero,
but our approach relaxed this constraint by incorporating the training loss as
an additional term in the bound. We demonstrated the effectiveness of our
method through various examples, including PINN and Separable PINN [43] (a
variant of PINN), where the activation functions are not ReLU. The results
demonstrated that our approach achieved good accuracy with low spike rates,
making it a promising and energy-efficient solution for scientific machine
learning. By enabling the conversion of a wider range of neural networks to
SNNs, this method opens up new possibilities for harnessing the temporal
processing capabilities and computational efficiency of SNNs in various
scientific applications. The proposed approach contributes to advancing the
field of spiking neural networks and their potential for practical real-world
implementations in edge computing. Future research can explore further
extensions of the calibration technique to other types of activation functions
and investigate its performance on more complex neural network architectures.

##  6 Acknowledgements

This work was supported in part by the DOE SEA-CROGS project (DE-SC0023191),
the ONR Vannevar Bush Faculty Fellowship (N00014-22-1-2795), CoCoSys- a
JUMP2.0 center sponsored by DARPA and SRC, and the National Science Foundation
CAREER Award.

## References

  * [1] Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis.  Physics-informed neural networks (PINNs) for fluid mechanics: A review.  Acta Mechanica Sinica, pages 1–12, 2022.
  * [2] Qianying Cao, Somdatta Goswami, George Em Karniadakis, and Souvik Chakraborty.  Deep neural operators can predict the real-time response of floating offshore structures under irregular waves.  arXiv preprint arXiv:2302.06667, 2023.
  * [3] Pengzhan Jin, Zhen Zhang, Aiqing Zhu, Yifa Tang, and George Em Karniadakis.  Sympnets: Intrinsic structure-preserving symplectic networks for identifying hamiltonian systems.  Neural Networks, 132:166–179, 2020.
  * [4] Enrui Zhang, Adar Kahana, Eli Turkel, Rishikesh Ranade, Jay Pathak, and George Em Karniadakis.  A hybrid iterative numerical transferable solver (hints) for pdes based on deep operator network and relaxation methods, 2022.
  * [5] Somdatta Goswami, Ameya D Jagtap, Hessam Babaee, Bryan T Susi, and George Em Karniadakis.  Learning stiff chemical kinetics using extended deep neural operators.  arXiv preprint arXiv:2302.12645, 2023.
  * [6] Qianying Cao, Somdatta Goswami, and George Em Karniadakis.  Lno: Laplace neural operator for solving differential equations.  arXiv preprint arXiv:2303.10528, 2023.
  * [7] Mitchell Daneker, Zhen Zhang, George Em Karniadakis, and Lu Lu.  Systems biology: Identifiability analysis and parameter identification via systems-biology informed neural networks.  arXiv preprint arXiv:2202.01723, 2022.
  * [8] Juan Diego Toscano, Christian Zuniga-Navarrete, Wilson David Jo Siu, Luis Javier Segura, and Hongyue Sun.  Teeth mold point cloud completion via data augmentation and hybrid rl-gan.  Journal of Computing and Information Science in Engineering, 23(4):041008, 2023.
  * [9] Qian Zhang, Konstantina Sampani, Mengjia Xu, Shengze Cai, Yixiang Deng, He Li, Jennifer K. Sun, and George Em Karniadakis.  AOSLO-net: A Deep Learning-Based Method for Automatic Segmentation of Retinal Microaneurysms From Adaptive Optics Scanning Laser Ophthalmoscopy Images.  Translational Vision Science & Technology, 11(8):7–7, 08 2022\.
  * [10] Zhongyi Jiang, Min Zhu, Dongzhuo Li, Qiuzi Li, Yanhua O Yuan, and Lu Lu.  Fourier-mionet: Fourier-enhanced multiple-input neural operators for multiphase modeling of geological carbon sequestration.  arXiv preprint arXiv:2303.04778, 2023.
  * [11] Min Zhu, Shihang Feng, Youzuo Lin, and Lu Lu.  Fourier-deeponet: Fourier-enhanced deep operator networks for full waveform inversion with improved accuracy, generalizability, and robustness.  arXiv preprint arXiv:2305.17289, 2023.
  * [12] Sheng Zhang, Joan Ponce, Zhen Zhang, Guang Lin, and George Karniadakis.  An integrated framework for building trustworthy data-driven epidemiological models: Application to the covid-19 outbreak in new york city.  PLOS Computational Biology, 17(9):1–29, 09 2021.
  * [13] Ehsan Kharazmi, Min Cai, Xiaoning Zheng, Zhen Zhang, Guang Lin, and George Em Karniadakis.  Identifiability and predictability of integer- and fractional-order epidemiological models using physics-informed neural networks.  Nature Computational Science, 1(11):744–753, Nov 2021.
  * [14] Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational physics, 378:686–707, 2019.
  * [15] Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Karniadakis.  Deepxde: A deep learning library for solving differential equations.  SIAM review, 63(1):208–228, 2021.
  * [16] Chenxi Wu, Min Zhu, Qinyang Tan, Yadhu Kartha, and Lu Lu.  A comprehensive study of non-adaptive and residual-based adaptive sampling for physics-informed neural networks.  Computer Methods in Applied Mechanics and Engineering, 403:115671, 2023.
  * [17] Sokratis J Anagnostopoulos, Juan Diego Toscano, Nikolaos Stergiopulos, and George Em Karniadakis.  Residual-based attention and connection to information bottleneck theory in pinns.  arXiv preprint arXiv:2307.00379, 2023.
  * [18] Jeremy Yu, Lu Lu, Xuhui Meng, and George Em Karniadakis.  Gradient-enhanced physics-informed neural networks for forward and inverse pde problems.  Computer Methods in Applied Mechanics and Engineering, 393:114823, 2022.
  * [19] Zongren Zou and George Em Karniadakis.  L-hydra: Multi-head physics-informed neural networks, 2023.
  * [20] Benjamin Fan, Edward Qiao, Anran Jiao, Zhouzhou Gu, Wenhao Li, and Lu Lu.  Deep learning for solving and estimating dynamic macro-finance models.  arXiv preprint arXiv:2305.09783, 2023.
  * [21] Enrui Zhang, Ming Dao, George Em Karniadakis, and Subra Suresh.  Analyses of internal structures and defects in materials using physics-informed neural networks.  Science Advances, 8(7):eabk0644, 2022.
  * [22] Maziar Raissi, Alireza Yazdani, and George Em Karniadakis.  Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations.  Science, 367(6481):1026–1030, 2020.
  * [23] Enrui Zhang, Minglang Yin, and George Em Karniadakis.  Physics-informed neural networks for nonhomogeneous material identification in elasticity imaging.  arXiv preprint arXiv:2009.04525, 2020.
  * [24] Yuyao Chen, Lu Lu, George Em Karniadakis, and Luca Dal Negro.  Physics-informed neural networks for inverse problems in nano-optics and metamaterials.  Optics Express, 28(8):11618–11633, 2020.
  * [25] Zongren Zou, Xuhui Meng, Apostolos F Psaros, and George Em Karniadakis.  Neuraluq: A comprehensive library for uncertainty quantification in neural differential equations and operators, 2022.
  * [26] Apostolos F. Psaros, Xuhui Meng, Zongren Zou, Ling Guo, and George Em Karniadakis.  Uncertainty quantification in scientific machine learning: Methods, metrics, and comparisons.  Journal of Computational Physics, 477:111902, 2023.
  * [27] Riccardo Massa, Alberto Marchisio, Maurizio Martina, and Muhammad Shafique.  An efficient spiking neural network for recognizing gestures with a dvs camera on the loihi neuromorphic processor.  In 2020 International Joint Conference on Neural Networks (IJCNN), pages 1–9. IEEE, 2020.
  * [28] Youngeun Kim, Hyoungseob Park, Abhishek Moitra, Abhiroop Bhattacharjee, Yeshwanth Venkatesha, and Priyadarshini Panda.  Rate coding or direct coding: Which one is better for accurate, robust, and energy-efficient spiking neural networks?  In ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 71–75. IEEE, 2022.
  * [29] Maxence Bouvier, Alexandre Valentian, Thomas Mesquida, Francois Rummens, Marina Reyboz, Elisa Vianello, and Edith Beigne.  Spiking neural networks hardware implementations and challenges: A survey.  ACM Journal on Emerging Technologies in Computing Systems (JETC), 15(2):1–35, 2019.
  * [30] Kaushik Roy, Akhilesh Jaiswal, and Priyadarshini Panda.  Towards spike-based machine intelligence with neuromorphic computing.  Nature, 575(7784):607–617, 2019.
  * [31] Mike Davies, Narayan Srinivasa, Tsung-Han Lin, Gautham Chinya, Yongqiang Cao, Sri Harsha Choday, Georgios Dimou, Prasad Joshi, Nabil Imam, Shweta Jain, et al.  Loihi: A neuromorphic manycore processor with on-chip learning.  Ieee Micro, 38(1):82–99, 2018.
  * [32] Garrick Orchard, E Paxon Frady, Daniel Ben Dayan Rubin, Sophia Sanborn, Sumit Bam Shrestha, Friedrich T Sommer, and Mike Davies.  Efficient neuromorphic signal processing with loihi 2.  In 2021 IEEE Workshop on Signal Processing Systems (SiPS), pages 254–259. IEEE, 2021.
  * [33] Adar Kahana, Qian Zhang, Leonard Gleyzer, and George Em Karniadakis.  Spiking neural operators for scientific machine learning, 2022.
  * [34] Qian Zhang, Adar Kahana, George Em Karniadakis, and Panos Stinis.  Sms: Spiking marching scheme for efficient long time integration of differential equations, 2022.
  * [35] Jayram Moorkanikara Nageswaran, Nikil Dutt, Yingxue Wang, and Tobi Delbrueck.  Computing spike-based convolutions on gpus.  In 2009 IEEE International Symposium on Circuits and Systems (ISCAS), pages 1917–1920. IEEE, 2009.
  * [36] Elie L Bienenstock, Leon N Cooper, and Paul W Munro.  Theory for the development of neuron selectivity: orientation specificity and binocular interaction in visual cortex.  Journal of Neuroscience, 2(1):32–48, 1982.
  * [37] Alan L Hodgkin, Andrew F Huxley, and Bernard Katz.  Measurement of current-voltage relations in the membrane of the giant axon of loligo.  The Journal of physiology, 116(4):424, 1952.
  * [38] Eugene M Izhikevich.  Dynamical systems in neuroscience.  MIT press, 2007.
  * [39] Wulfram Gerstner and Werner M Kistler.  Spiking neuron models: Single neurons, populations, plasticity.  Cambridge university press, 2002.
  * [40] Wulfram Gerstner, Raphael Ritz, and J Leo Van Hemmen.  Why spikes? hebbian learning and retrieval of time-resolved excitation patterns.  Biological cybernetics, 69(5-6):503–515, 1993.
  * [41] David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams.  Learning representations by back-propagating errors.  nature, 323(6088):533–536, 1986.
  * [42] Timothy P Lillicrap, Adam Santoro, Luke Marris, Colin J Akerman, and Geoffrey Hinton.  Backpropagation and the brain.  Nature Reviews Neuroscience, 21(6):335–346, 2020.
  * [43] Junwoo Cho, Seungtae Nam, Hyunmo Yang, Seok-Bae Yun, Youngjoon Hong, and Eunbyung Park.  Separable pinn: Mitigating the curse of dimensionality in physics-informed neural networks.  arXiv preprint arXiv:2211.08761, 2022.
  * [44] Bodo Rueckauer, Iulia-Alexandra Lungu, Yuhuang Hu, Michael Pfeiffer, and Shih-Chii Liu.  Conversion of continuous-valued deep networks to efficient event-driven networks for image classification.  Frontiers in neuroscience, 11:682, 2017.
  * [45] Yuhang Li, Shikuang Deng, Xin Dong, Ruihao Gong, and Shi Gu.  A free lunch from ann: Towards efficient, accurate spiking neural networks calibration.  In International Conference on Machine Learning, pages 6316–6325. PMLR, 2021.
  * [46] Yujie Liu, Zhe Wang, Peng Li, Yiran Chen, and Hai Li.  Optimal ann-snn conversion for fast and accurate inference in deep spiking neural networks.  arXiv preprint arXiv:2105.11654, 2021.
  * [47] Jihwan Kim, Jangho Lee, and Junmo Lee.  Reducing ann-snn conversion error through residual learning.  arXiv preprint arXiv:2302.02091, 2021.
  * [48] Bodo Rueckauer, Iulia-Alexandra Lungu, Yuhuang Hu, Michael Pfeiffer, and Shih-Chii Liu.  Spiking neural network conversion toolbox.  https://github.com/NeuromorphicProcessorProject/snn_toolbox, 2017\.
  * [49] Yuhang Li, Shikuang Deng, Xin Dong, and Shi Gu.  Converting artificial neural networks to spiking neural networks via parameter calibration, 2022.
  * [50] Zhen Liu, Zhihui Li, Zhihong Li, Yujie Wang, Zhen Li, Zhihui Li, Zhihong Li, and Yujie Wang.  Backeisnn: A deep spiking neural network with adaptive self-feedback and balanced excitatory–inhibitory neurons.  Neural Networks, 138:1–14, 2021.
  * [51] Xiaowei Jin, Shengze Cai, Hui Li, and George Em Karniadakis.  Nsfnets (navier-stokes flow nets): Physics-informed neural networks for the incompressible navier-stokes equations.  Journal of Computational Physics, 426:109951, 2021.
  * [52] Xin-Yang Liu, Hao Sun, Min Zhu, Lu Lu, and Jian-Xun Wang.  Predicting parametric spatiotemporal dynamics by multi-resolution pde structure-preserved deep learning.  arXiv preprint arXiv:2205.03990, 2022.
  * [53] Min Zhu, Handi Zhang, Anran Jiao, George Em Karniadakis, and Lu Lu.  Reliable extrapolation of deep neural operators informed by physics or sparse observations.  Computer Methods in Applied Mechanics and Engineering, 412:116064, 2023.
  * [54] John Kim and Parviz Moin.  Application of a fractional-step method to incompressible navier-stokes equations.  Journal of computational physics, 59(2):308–323, 1985.

## Appendix

###### Lemma 1

If the activation function is piecewise linear, then the conversion error in
the final network output space can be bounded by a weighted sum of local
conversion errors, given by

$\mathbf{e}^{(n),\top}\mathbf{H}^{(n)}\mathbf{e}^{(n)}\leq\sum_{l=1}^{n}2^{n-l+1}\mathbf{e}_{c}^{(l),\top}\mathbf{H}^{(l)}\mathbf{e}_{c}^{(l)}$
(10)

_Proof_ In the proof of [49] uses only one property of ReLU that its second-
order derivatives are all zeros. And it does not require the activation to be
non-negative. So it applies to all piecewise linear activation functions.

###### Lemma 2

The Hessian matrix of activation in layer $l$ can be written as

$\mathbf{H}^{(l)}=\mathbf{W}^{(l),\top}\mathbf{B}^{(l),\top}\mathbf{H}^{(l+1)}\mathbf{B}^{(l)}\mathbf{W}^{(l)}+\mathbf{W}^{(l),\top}\mathbf{C}^{(l)}\mathbf{W}^{(l)}$
(11)

where $W^{(l)}$ is the weight of layer $l$,

$\mathbf{B}^{(l)}=\mathrm{diag}\left\\{\frac{\partial\mathbf{x}^{(l+1)}_{1}}{\partial\mathbf{z}^{(l)}_{1}},\ldots,\frac{\partial\mathbf{x}^{(l+1)}_{n_{l}}}{\partial\mathbf{z}^{(l)}_{n_{l}}}\right\\}$
(12)

is a diagonal matrix, each element on diagonal is the derivative of activation
function, and

$\mathbf{C}^{(l)}=\mathrm{diag}\left\\{\frac{\partial
L}{\partial\mathbf{x}^{(l+1)}_{1}}f^{\prime\prime}(z^{(l)}_{1}),\ldots,\frac{\partial
L}{\partial\mathbf{x}^{(l+1)}_{n_{l}}}f^{\prime\prime}(z^{(l)}_{n_{l}})\right\\}$
(13)

here $x^{(l+1)}$ is the input of layer $(l+1)$ and $z^{(l)}$ is the
unactivated output of layer $(l)$.

_Proof_ Similar to the proof in the [49], we calculate
$\mathbf{H}^{(l)}_{a,b}$ by

$\displaystyle\mathbf{H}^{(l)}_{a,b}$
$\displaystyle=\frac{\partial^{2}L}{\partial\mathbf{x}^{(l)}_{a}\partial\mathbf{x}^{(l)}_{b}}=\frac{\partial}{\partial\mathbf{x}^{(l)}_{b}}\left(\sum_{i}\frac{\partial
L}{\partial\mathbf{x}^{(l+1)}_{i}}\frac{\partial\mathbf{x}^{(l+1)}_{i}}{\partial\mathbf{x}^{(l)}_{a}}\right)$
(14)
$\displaystyle=\sum_{i}\frac{\partial}{\partial\mathbf{x}^{(l)}_{b}}\left(\frac{\partial
L}{\partial\mathbf{x}^{(l+1)}_{i}}\frac{\partial\mathbf{x}^{(l+1)}_{i}}{\partial\mathbf{z}^{(l)}_{i}}\frac{\partial\mathbf{z}^{(l)}_{i}}{\partial\mathbf{x}^{(l)}_{a}}\right)$
$\displaystyle=\sum_{i}\mathbf{W}^{(l)}_{i,a}\frac{\partial}{\partial\mathbf{x}^{(l)}_{b}}\left(\frac{\partial
L}{\partial\mathbf{x}^{(l+1)}_{i}}\frac{\partial\mathbf{x}^{(l+1)}_{i}}{\partial\mathbf{z}^{(l)}_{i}}\right)$
$\displaystyle=\sum_{i}\mathbf{W}^{(l)}_{i,a}\frac{\partial\mathbf{x}^{(l+1)}_{i}}{\partial\mathbf{z}^{(l)}_{i}}\frac{\partial}{\partial\mathbf{x}^{(l)}_{b}}\frac{\partial
L}{\partial\mathbf{x}^{(l+1)}_{i}}+\sum_{i}\mathbf{W}^{(l)}_{i,a}\frac{\partial
L}{\partial\mathbf{x}^{(l+1)}_{i}}\frac{\partial}{\partial\mathbf{x}^{(l)}_{b}}\frac{\partial\mathbf{x}^{(l+1)}_{i}}{\partial\mathbf{z}^{(l)}_{i}}$
$\displaystyle=\sum_{i,j}\mathbf{W}^{(l)}_{i,a}\frac{\partial\mathbf{x}^{(l+1)}_{i}}{\partial\mathbf{z}^{(l)}_{i}}\frac{\partial^{2}L}{\partial\mathbf{x}^{(l+1)}_{j}\partial\mathbf{x}^{(l+1)}_{i}}\frac{\partial\mathbf{x}^{(l+1)}_{j}}{\partial\mathbf{z}^{(l)}_{j}}\mathbf{W}^{(l)}_{j,b}+\sum_{i}\mathbf{W}^{(l)}_{i,a}\frac{\partial
L}{\partial\mathbf{x}^{(l+1)}_{i}}f^{\prime\prime}(\mathbf{z}_{i})\mathbf{W}^{(l)}_{i,b}$

###### Lemma 3

Let $M\in\mathbb{R}^{p\times p}$ and $x_{n}\in\mathbb{R}^{p}$. If
$\norm{x_{n}}<K_{x}$ $\forall n$, $\norm{x_{n}-x_{m}}<\epsilon$, then
$|x_{n}^{\top}Mx_{n}-x_{m}^{\top}Mx_{m}|<2K_{M}K_{x}\epsilon$, where
$K_{M}=\norm{M}\stackrel{{\scriptstyle
d}}{{=}}\sup_{\norm{x}=\norm{y}=1}x^{\top}Ky$.

_Proof_ By the triangular inequality and definition of $K_{M}$,

$\displaystyle|x_{n}^{\top}Mx_{n}-x_{m}^{\top}Mx_{m}|$
$\displaystyle\leq|x_{n}^{\top}Mx_{n}-x_{m}^{\top}Mx_{n}|+|\leq|x_{m}^{\top}Mx_{n}-x_{m}^{\top}Mx_{m}|$
(15) $\displaystyle\leq
K_{M}\norm{x_{n}}\norm{x_{n}-x_{m}}+K_{M}\norm{x_{m}}\norm{x_{n}-x_{m}}$
$\displaystyle<2K_{M}K_{x}\epsilon$

###### Lemma 4

Let $x\in\mathbb{R}^{p}$. If $\norm{x}<\epsilon$, then
$|x^{\top}\mathbf{C}^{(l)}x|<K_{l}\sqrt{L}\epsilon$, where $K_{l}$ is a
constant and $L$ is the training loss.

_Proof_ Let $y=\mathbf{W}^{(l)}x$, then $\exists K_{W}>0$, s.t. $\norm{y}\leq
K_{W}\epsilon$. And since we are using MSE for the loss,

$\frac{\partial L}{\partial\mathbf{x}^{(l+1)}_{i}}=\sum_{j}\frac{\partial
L}{\partial\mathbf{x}^{(n)}_{j}}\frac{\partial\mathbf{x}^{(n)}_{j}}{\partial\mathbf{x}^{(l+1)}_{i}}\leq\sum_{j}K_{J}\sqrt{L}\frac{\partial\mathbf{x}^{(n)}_{j}}{\partial\mathbf{x}^{(l+1)}_{i}}\leq
K\sqrt{L}$ (16)

for some $K_{J},K>0$, because of the boundedness of $f^{\prime}$ and
$\mathbf{W}$. And $\norm{\mathbf{C}^{(n)}}<\infty$ by the boundedness of
$f^{\prime\prime}$, then by definition of $\norm{\mathbf{C}}$ we have the
result.

With the preceding lemmas, now we can prove our main result.

_Proof_ Suppose the activation function is $f$. By the approximation property
of $f$, $\forall\epsilon>0$, $\exists\tilde{f}$, such that

$\norm{f-\tilde{f}}_{\infty}<\epsilon,\quad\norm{f^{\prime}-\tilde{f}^{\prime}}_{\infty}<\epsilon$
(17)

Let $\tilde{\cdot}$ be the values computed by the approximating activation
function, like $\tilde{\mathbf{x}}^{(l)}$ be the output of layer $l$ computed
by the approximating activation function $\tilde{f}$. Define
$\tilde{\mathbf{e}}$ as the error computed with $\tilde{\mathbf{x}}$ instead
of $\mathbf{x}$. Then we can assume

$\norm{\tilde{\mathbf{e}}-\mathbf{e}}<\epsilon,\quad\norm{\tilde{\mathbf{e}_{c}}-\mathbf{e}_{c}}<\epsilon$
(18)

by (17).

First we can approximate
$\mathbf{e}^{(n),\top}\tilde{\mathbf{H}}^{(n)}\mathbf{e}^{(n)}$ by
$\tilde{\mathbf{e}}^{(n),\top}\mathbf{H}^{(n)}\tilde{\mathbf{e}}^{(n)}$. This
follows from Lemma 3 and the fact that
$\tilde{\mathbf{H}}^{(n)}=\mathbf{H}^{(n)}$ as the last layer does not have an
activation function to convert.

Then we need to approximate
$\mathbf{e}_{c}^{(l),\top}\mathbf{H}^{(l)}\mathbf{e}_{c}^{(l)}$.

$|\mathbf{e}_{c}^{(l),\top}\mathbf{H}^{(l)}\mathbf{e}_{c}^{(l)}-\tilde{\mathbf{e}_{c}}^{(l),\top}\tilde{\mathbf{H}}^{(l)}\tilde{\mathbf{e}_{c}}^{(l)}|\leq|\mathbf{e}_{c}^{(l),\top}\mathbf{H}^{(l)}\mathbf{e}_{c}^{(l)}-\tilde{\mathbf{e}_{c}}^{(l),\top}\mathbf{H}^{(l)}\tilde{\mathbf{e}_{c}}^{(l)}|+|\tilde{\mathbf{e}_{c}}^{(l),\top}\mathbf{H}^{(l)}\tilde{\mathbf{e}_{c}}^{(l)}-\tilde{\mathbf{e}_{c}}^{(l),\top}\tilde{\mathbf{H}}^{(l)}\tilde{\mathbf{e}_{c}}^{(l)}|$
(19)

The first term on the right hand side is easy to bound by Lemma 3. To
investigate the second term, we compute $\mathbf{H}^{(l)}$ recursively from
$l=n-1$. By Lemma 2,

$\displaystyle\mathbf{H}^{(n-1)}$
$\displaystyle=\mathbf{W}^{(n-1),\top}\mathbf{B}^{(n-1),\top}\mathbf{H}^{(n)}\mathbf{B}^{(n-1)}\mathbf{W}^{(n-1)}+\mathbf{W}^{(n-1),\top}\mathbf{C}^{(n-1)}\mathbf{W}^{(n-1)}$
(20)
$\displaystyle=\mathbf{W}^{(n-1),\top}(\mathbf{B}^{(n-1),\top}\mathbf{H}^{(n)}\mathbf{B}^{(n-1)}-\tilde{\mathbf{B}}^{(n-1),\top}\mathbf{H}^{(n)}\tilde{\mathbf{B}}^{(n-1)})\mathbf{W}^{(n-1)}$
$\displaystyle\quad+\mathbf{W}^{(n-1),\top}\tilde{\mathbf{B}}^{(n-1),\top}\mathbf{H}^{(n)}\tilde{\mathbf{B}}^{(n-1)}\mathbf{W}^{(n-1)}+\mathbf{W}^{(n-1),\top}\mathbf{C}^{(n-1)}\mathbf{W}^{(n-1)}$
$\displaystyle=\mathbf{W}^{(n-1),\top}\Delta\tilde{\mathbf{B}}^{(n-1)}\mathbf{W}^{(n-1)}+\tilde{\mathbf{H}}^{(n-1)}+\mathbf{W}^{(n-1),\top}\mathbf{C}^{(n-1)}\mathbf{W}^{(n-1)}$

where
$\Delta\tilde{\mathbf{B}}^{(n-1)}=\mathbf{B}^{(n-1),\top}\mathbf{H}^{(n)}\mathbf{B}^{(n-1)}-\tilde{\mathbf{B}}^{(n-1),\top}\mathbf{H}^{(n)}\tilde{\mathbf{B}}^{(n-1)}$.
Then by Lemma 3 and 4, we have
$\mathbf{H}^{(n-1)}-\tilde{\mathbf{H}}^{(n-1)}\approx
K_{n-1}\sqrt{L}\mathbf{I}$ in the sense of bilinear form. For $n=l-2$,

$\displaystyle\mathbf{H}^{(n-2)}$
$\displaystyle=\mathbf{W}^{(n-2),\top}\mathbf{B}^{(n-2),\top}\mathbf{H}^{(n-1)}\mathbf{B}^{(n-2)}\mathbf{W}^{(n-2)}+\mathbf{W}^{(n-2),\top}\mathbf{C}^{(n-2)}\mathbf{W}^{(n-2)}$
(21)
$\displaystyle=\tilde{\mathbf{H}}^{(n-2)}+\mathbf{W}^{(n-2),\top}\mathbf{B}^{(n-2),\top}\mathbf{H}^{(n-1)}\mathbf{B}^{(n-2)}\mathbf{W}^{(n-2)}+\mathbf{W}^{(n-2),\top}\mathbf{C}^{(n-2)}\mathbf{W}^{(n-2)}$
$\displaystyle\quad-\mathbf{W}^{(n-2),\top}\tilde{\mathbf{B}}^{(n-2),\top}\tilde{\mathbf{H}}^{(n-1)}\tilde{\mathbf{B}}^{(n-2)}\mathbf{W}^{(n-2)}$
$\displaystyle=\tilde{\mathbf{H}}^{(n-2)}+\mathbf{W}^{(n-2),\top}\Delta\tilde{\mathbf{B}}^{(n-2)}\mathbf{W}^{(n)}+\mathbf{W}^{(n-2),\top}\tilde{\mathbf{B}}^{(n-2),\top}\mathbf{W}^{(n-1),\top}\Delta\tilde{\mathbf{B}}^{(n-1)}\mathbf{W}^{(n-1)}\tilde{\mathbf{B}}^{(n-2)}\mathbf{W}^{(n-2)}$
$\displaystyle\quad+\mathbf{W}^{(n-2),\top}\tilde{\mathbf{B}}^{(n-2),\top}\mathbf{W}^{(n-1),\top}\mathbf{C}^{(n-1)}\mathbf{W}^{(n-1)}\tilde{\mathbf{B}}^{(n-2)}\mathbf{W}^{(n-2)}+\mathbf{W}^{(n-2),\top}\mathbf{C}^{(n-2)}\mathbf{W}^{(n-2)}$

Similarly, we have $\mathbf{H}^{(n-2)}-\tilde{\mathbf{H}}^{(n-2)}\approx
K_{n-2}\sqrt{L}\mathbf{I}$ in the sense of bilinear form. Thus
$\mathbf{H}^{(l)}-\tilde{\mathbf{H}}^{(l)}\approx K_{l}\sqrt{L}\mathbf{I}$
$\forall l=1,2,\ldots,n-1$ by induction. Hence we can approximate
$\tilde{\mathbf{e}_{c}}^{(l),\top}\tilde{\mathbf{H}}^{(l)}\tilde{\mathbf{e}_{c}}^{(l)}$
by
$\mathbf{e}_{c}^{(l),\top}(\mathbf{H}^{(l)}+K_{l}\sqrt{L}\mathbf{I})\mathbf{e}_{c}^{(l)}$.

Finally, because both sides of the inequality in Lemma 10 are well
approximated, the proof is completed.