# The kinetic theory of ultra-subsonic fermion systems and applications to
flat band magic angle twisted bilayer graphene

Seth M. Davis<EMAIL_ADDRESS>Condensed Matter Theory Center and Joint
Quantum Institute, Department of Physics, University of Maryland, College
Park, MD 20742, USA  Sankar Das Sarma  Condensed Matter Theory Center and
Joint Quantum Institute, Department of Physics, University of Maryland,
College Park, MD 20742, USA

###### Abstract

The only kinematically-allowed phonon-scattering events for bands of subsonic
fermions ($v_{F}<v_{p}$) are interband transitions, leading to different
low-$T$ transport physics than in the typical supersonic case. We apply a
kinetic theory of phonon-limited transport to a generic two-band system of
subsonic fermions, deriving formulae for relaxation times and resistivity that
are accurate in the limit of close bands and small $v_{F}/v_{p}$. We predict
regimes of $\rho\propto T$, $\rho\propto T^{2}$, and perfect conductivity. Our
theory predicts linear-in-$T$ resistivity down to a crossover temperature that
is suppressed from its supersonic analogue by a factor of $v_{F}/v_{p}$,
offering a new explanation for low-$T$ “strange metal” behavior observed in
flat band systems.

Rapid progress in the fabrication and manipulation of layered two-dimensional
van der Waals heterostructures has lead to an unprecedented ability to
engineer nearly flat band (NFB) electronic systems, which have already
displayed a wealth of exotic phenomena Geim and Grigorieva (2013); Novoselov
et al. (2006); Bistritzer and MacDonald (2011); Morell et al. (2010); Li et
al. (2019); Kim et al. (2017); Cao et al. (2018a, b, 2020a, 2021); Yankowitz
et al. (2019); Kerelsky et al. (2019); Lu et al. (2019); Stepanov et al.
(2020); Sharpe et al. (2019); Chen et al. (2020); Rozen et al. (2021); Zhou et
al. (2022, 2021a, 2021b); Serlin et al. (2020); Wu et al. (2018, 2019a);
Tschirhart et al. (2022); Polshyn et al. (2020); Jaoui et al. (2021); Polshyn
et al. (2019); Cao et al. (2020b); Sarma and Wu (2022); Zhang et al. (2022);
Polski et al. (2022); Xie and MacDonald (2020); Andrei and MacDonald (2020);
Li et al. (2021); Ghiotto et al. (2021); Pan et al. (2020); Pan and Sarma
(2021); Morales-Durán et al. (2021); Ahn and Sarma (2022); Kerelsky et al.
(2021); Khalaf et al. (2019); Arora et al. (2020). However, all solid state
systems contain phonons. To understand observations of novel physics in solid
state NFB systems, it is important to understand how phonons interact with NFB
fermions. This work focuses on one aspect: when the fermions in question are
subsonic ($v_{F}<v_{p}$), kinematics requires that all single-phonon
scattering processes are interband transitions, with consequences on the
low-$T$ transport physics.

The prime example of the NFB systems is magic angle twisted bilayer graphene
(MATBLG). MATBLG has been found to exhibit SC proximate to strongly correlated
insulating states Cao et al. (2018a, b, 2020a, 2021); Yankowitz et al. (2019);
Kerelsky et al. (2019); Lu et al. (2019); Wu et al. (2018) and has been
reported to exhibit a linear-in-$T$ “strange metal”-like resistivity over a
large range of dopings and temperatures Jaoui et al. (2021); Polshyn et al.
(2019); Cao et al. (2020b); Sarma and Wu (2022), sometimes down to
temperatures as low as $50mK$ Jaoui et al. (2022). These phenomena have
inspired analogy between MATBLG and the cuprate high-$T_{c}$ SCs, as well as
speculation that SC in MATBLG might be driven by strong correlation physics.
However, phonon-based theories of SC Wu et al. (2019b); Li et al. (2020) and
high-$T$ transport Wu et al. (2019b); Davis et al. (2023a); Sarma and Wu
(2022) in MATBLG have been put forth that give generally good quantitative
agreement with experiment. These results - along with many others Stepanov et
al. (2020); Saito et al. (2020) \- suggest an alternative physical picture in
which MATBLG hosts standard, phonon-driven BCS-style SC that competes with
interaction-driven insulating orders at commensurate filling fractions of the
NFBs Lu et al. (2019). It is thus imperative to understand whether low-$T$,
linear-in-$T$ transport in MATBLG is indeed arising from a strange metal
state.

The standard kinetic theory of acoustic phonon scattering is exceptionally
accurate in describing transport in layered graphene systems (as well as in
normal metals and semiconductors) at temperatures above a few Kelvin Hwang and
Das Sarma (2008); Wu et al. (2019b); Davis et al. (2023b, a), and accurately
describes TBLG transport away from the magic angle Wu et al. (2019b); Davis et
al. (2023a). In this work, we apply the same framework to NFB systems at
asymptotically low $T$, extending the kinetic theory well beyond its regime of
proven validity. Remarkably, the familiar theory predicts qualitatively
different low-$T$ transport physics for subsonic fermions than it does for the
supersonic alternative, due entirely to the kinematic differences between the
two limits.

We develop a transport theory for an “ultrasubsonic” (USS) two-band fermion
system, which is accurate in the asymptotic double-limit of small
$v_{F}/v_{p}$ and small band separation. Our main results are analytical
formulae for relaxation times and for the resistivity of a general USS system,
with transparent dependencies on temperature and doping level. We show that
the interband nature of subsonic fermion scattering manifests an exponentially
divergent relaxation time at low temperatures, in contrast with the familiar
$\tau\propto T^{-4}$ Bloch-Grüneisen (BG) power law applicable to supersonic
fermions. The divergence in relaxation time is capable of perfectly balancing
the thermodynamic suppression of states away from the Fermi level, leading to
a non-intuitive physical picture in which states far from the Fermi level
contribute meaningfully to transport. This scenario predicts a linear-in-$T$
resistivity over a wide range of $T$, and down to a temperature much lower
than the BG paradigm of supersonic fermion bands would suggest possible. The
divergence in relaxation time can also drive the system to perfect
conductivity at asymptotically low $T$, mimicking a SC transition. We also
note that the geometry of isolated NFBs can manifest a mid-$T$ $\rho\propto
T^{2}$ power law. All these features are consistent with the hitherto
unexplained phenomenology of MATBLG transport.

Figure 1: Key aspects of phonon scattering for subsonic fermions can be
understood from simple schematic drawings. We consider a cross section of a
fermion band structure (blue), with momentum on the horizontal axis and energy
on the vertical axis. We pick an “initial” Bloch state and superimpose on it
the Debye dispersion of the acoustic phonons (orange). Due to energy and
momentum conservation, a fermion in the initial state can only scatter to a
new state at an intersection point of the fermion and phonon bands. In (A), we
depict the standard case of phonon scattering for supersonic (Dirac) fermions.
The maximum allowed energy transfer defines $T_{BG}$, while arbitrarily small
energy transfers are allowed. In (B), we show that a single-band system of
subsonic fermions does not allow single-phonon scattering events. This follows
from the simple fact that the fermion and phonon dispersions only intersect at
the initial state. In (C) we consider a multiband subsonic fermion system,
where phonon scattering is necessarily interband. In addition to a maximum
allowed energy transfer defining a Bloch-Grüneisen scale, there is also a
minimum allowed energy transfer.

There are reasons to suspect that the standard kinetic theory of acoustic-
phonon-limited-transport may not apply to NFB systems. Kinetic theory replaces
the interacting Hamiltonian with a Markovian collision integral, and is thus
incapable of treating emergent physics due to strong correlations, which many
expect to be important in NFB systems. Further, vertex corrections to the
electron-phonon coupling could become important in the absence of Migdal’s
theorem. Nevertheless, it is important to understand the predictions of naive
kinetic theory for NBF systems. The fact that our theory gives a simple and
concrete mechanism for robust linear-in-$T$ resistivity is noteworthy, given
that this is the primary signature of the “strange metal” state. Further, it
is interesting that phonon scattering in subsonic fermion bands can also
generate regimes of $\rho\propto T^{2}$ resistivity - often assumed to arise
from electron-electron scattering in Fermi liquids - and regimes of perfect
conductivity that mimic the behavior of superconductivity. While a realistic
treatment of NFB electron-phonon physics is beyond reach, it is possible that
the atypical kinematics of subsonic fermions is a key ingredient to
understanding transport in NFB systems, which our approach is able to fully
capture.

Although our theory has MATBLG in mind, we work with a general USS model. Our
transport theory could apply to a system with renormalized fermion and phonon
band structures and a renormalized deformation potential coupling in a theory
that has already taken interaction-induced corrections into account. This, in
fact, is the standard transport theory in normal metals where electron-
electron interactions generally do not affect transport in the phonon-
dominated regime. Further, in addition to providing an alternative explanation
for experimental transport observations in MATBLG, our results could find
application generally in all subsonic band systems, including other flat band
systems generated by layered 2D heterostructures and heavy fermion systems. It
is also noteworthy that strange metallicity is ubiquitous in heavy fermion
compounds, which also manifest flat bands (i.e. heavy mass) and multiple bands
near the Fermi surface, indicating a possible role for our proposed subsonic
transport mechanism.

Phonon scattering of supersonic and subsonic fermions.—In an electron-phonon
scattering process, conservation of energy and momentum define a “scattering
manifold” of kinematically allowed Bloch states that an electron can scatter
to. The maximum energy difference between a Bloch state on the scattering
manifold and the initial state defines the Bloch-Grüneisen temperature,
$T_{BG}$ Hwang and Sarma (2019). For example, in an isotropic system, we
simply have $k_{B}T_{BG}\equiv 2v_{p}k_{F}$ [Fig.1 (A)]. When $T\gg T_{BG}$,
enough phonon modes will be populated that all kinematically allowed
scattering events are possible. This is the so-called “equipartition regime”,
which is characterized by a linear-in-$T$ scattering rate for each Bloch
state. This is universal in kinetic theory and usually (but not always Davis
et al. (2023b, a)) gives a linear-in-$T$ resistivity above $T_{BG}$. On the
other hand, when $T\ll T_{BG}$, only low-energy phonon modes are available and
scattering is restricted to a small neighborhood of the initial state. This is
what gives rise to the famous low-$T$ $\tau\propto T^{-(d+2)}$ power law in
the Bloch-Grüneisen regime (where $d$ is the spatial dimension). Crucially,
when the electrons are supersonic, the scattering manifold is smoothly
connected to the original state, allowing arbitrarily small-momentum
scattering events. As a result, the BG regime holds all the way to zero
temperature. In single-layer graphene and in normal metals,
$v_{F}/v_{p}\approx\mathcal{O}(10^{2})$ or larger. (We also note that in many
normal metals, the Debye temperature is lower than $T_{BG}$, in which case the
BG crossover takes place at the Debye temperature Hwang and Sarma (2019).)

However, when the Bloch state in question is subsonic, the scattering manifold
determined by energy-momentum conservation is necessarily disconnected from
the original state. Similarly, in a band made entirely of subsonic Bloch
states, all phonon scattering processes in subsonic fermion systems are
interband processes. This implies that the scattering manifold also has a
minimum allowed energy transfer [see Fig. 1 (C)], defining another temperature
scale applicable for subsonic scattering, $T_{SS}.$ Below $T_{SS}$ we expect
an crossover to a regime in which the entire scattering manifold is
inaccessible. That is, all thermally populated phonons are too low-energy to
satisfy the conservation laws of an electron-phonon scattering event.
Intuitively, we expect that below $T_{SS}$, the relaxation time of the Bloch
state diverges exponentially as $T\rightarrow 0$ (instead of $\tau\propto
T^{-4}$ as in the supersonic BG regime).

Basic kinetic theory.—We first review the fundamental equations for phonon-
limited resistivity in the Boltzmann framework Davis et al. (2023b, a). The
resistivity is given in terms of the relaxation times $\tau_{S}$:

$\displaystyle\frac{\delta^{ij}}{\rho}$
$\displaystyle=\frac{1}{4}\frac{e^{2}}{k_{B}T}\frac{1}{\mathcal{A}}\sum_{S}\frac{\tau_{S}v^{i}_{S}v^{j}_{S}}{\cosh\left(\frac{\varepsilon-\mu}{2k_{B}T}\right)^{2}},$
(1)

where $\mathcal{A}$ is the area of the system and $v_{S}$ is the Fermi
velocity of the state $S$. The relaxation times are determined from a self-
consistent integral equation, derived from the Boltzmann equation. In the case
of longitudinal acoustic phonons, in the Debye approximation, which couple to
the fermions via the deformation potential, this takes the form

$\displaystyle\frac{\pi
D^{2}}{\hbar\rho_{M}v_{p}^{2}}\frac{1}{\mathcal{A}}\sum_{S^{\prime}}\tilde{\Delta}_{S,S^{\prime}}\mathcal{C}_{S,S^{\prime}}\mathcal{F}_{S,S^{\prime}}^{\mu,T}\left[\tau_{S}-\frac{v_{S^{\prime}}}{v_{S}}\tau_{S^{\prime}}\cos\theta_{\bm{\mathrm{v}}}\right]$
$\displaystyle=1.$ (2)

Above, $\rho_{M}$ is the mass density of the system, $D$ is the deformation
potential, and $v_{p}$ is the phonon velocity. The function
$\tilde{\Delta}_{S,S^{\prime}}$ enforces the conservation of energy and
momentum and defines the scattering manifold, while
$\mathcal{C}_{S,S^{\prime}}$ are “matrix elements” encoding geometric
wavefunction overlap data Hwang and Das Sarma (2008); Wu et al. (2019b); Davis
et al. (2023b, a). The dependencies on temperature and chemical potential are
contained in $\mathcal{F}_{S,S^{\prime}}^{\mu,T}$, which encodes the thermal
occupation data for the fermion and phonon states. These are given explicitly
in the SM SM (Please see Supplemental Material for the microscopic model,
additional details on Boltzmann kinetic theory and derivations of
ultrasubsonic transport theory results, an analysis of the applicability of
our theory to MATBLG, and a treatment of the ultrasubsonic Dirac cone). On a
finite momentum grid, Eq. (2) is a matrix inversion problem that can be solved
for the relaxation times, $\tau_{S}$. Equation (1) then gives the resistivity.

Figure 2: In (A) we depict the scenario of an NFB system isolated from the
rest of the band structure by a sizable energy gap. Phonon scattering can
produce a $\rho\propto T^{2}$ resistivity power law for $\Delta
E_{\text{band}}\ll k_{B}T\ll\Delta E_{\text{gap}}$ in this band geometry. In
(B), we depict the band geometry of the ultrasubsonic limit, defined by
$v_{F}\ll v_{p}$ and small band separation.

Equipartition regime: $\rho\propto T^{2}$.—When $T\gg T_{BG}$ (equipartition
regime), the inverse relaxation times scale linearly with temperature:

$\displaystyle\tau_{S}=\frac{c_{S}}{k_{B}T}+\mathcal{O}\left(\frac{\Delta\varepsilon}{(k_{B}T)^{3}}\right),$
(3)

for some $(T,\mu)$-independent proportionality constants $c_{S}$ that depend
on the details of the system. In the case of NFBs separated from all other
bands by a large energy gap, there is a mid-$T$ regime in which $T$ is larger
than the NFB bandwidth, yet $T$ is still small compared to the gap between the
NFBs and the other bands [Fig. 2 (A)]. In this case, we may neglect bands
other than the NFBs, apply the equipartition scaling for the relaxation times
[Eq. (3)], and expand the thermal weighting functions ($\cosh$) in the formula
for the resistivity [Eq. (1)], giving

$\displaystyle\frac{\delta^{ij}}{\rho}$
$\displaystyle=\frac{1}{4}\frac{e^{2}}{(k_{B}T)^{2}}\frac{1}{\mathcal{A}}\sum_{S}c_{S}v^{i}_{S}v^{j}_{S}+\mathcal{O}\left(\frac{1}{(k_{B}T)^{4}}\right).$
(4)

We thus find a mid-$T$, $\rho\propto T^{2}$ scaling regime due entirely to
phonon scattering. This is noteworthy since $\rho\propto T^{2}$ scaling is
generally seen as the hallmark of transport dominated by electron-electron
scattering in a Fermi liquid Ashcroft and Mermin (1976); Ziman (1960); Coleman
(2015), but the “mid-$T$” regime generated by NFB phonon scattering has
exactly the same $T^{2}$ dependence.

Ultrasubsonic kinetic theory.—We now consider electron-phonon scattering in a
limit that allows transparent analytical results. We consider a two-band model
in the double-limit of small $v_{F}/v_{p}$ and small band separation, which we
call the “ultrasubsonic (USS) limit” [see Fig. 2 (B)]. For a given Bloch state
($S$, with energy $\varepsilon$), its scattering manifold is a small loop on
the opposite band, surrounding its “compliment”, the point on the opposite
band with the same momentum as the original state ($\bar{S}$, with energy
$\bar{\varepsilon}$). In the USS limit, the scattering manifold is well-
approximated by a circle of radius $|\varepsilon-\bar{\varepsilon}|/(\hbar
v_{p})$, centered at $\bar{S}$. The energy variation along the scattering
manifold is then roughly $(v_{F}/v_{p})|\varepsilon-\bar{\varepsilon}|,$ which
is double-suppressed in the USS limit. We thus make the approximations
$\mathcal{F}^{\mu,T}(S,S^{\prime})\approx\mathcal{F}^{\mu,T}(S,\bar{S})$ and
$\tau(S^{\prime})\approx\tau(\bar{S})$ for $S^{\prime}$ along the scattering
manifold, approximating the relaxation times and the thermal occupancy
function along the scattering manifold by their values at the compliment
point.

Using the approximations discussed above and taking the thermodynamic limit
(continuum limit in momentum space) of Eq. (2), we have

$\displaystyle 1$
$\displaystyle=\frac{D^{2}}{2\rho_{M}\hbar^{3}v_{p}^{4}}|\varepsilon-\bar{\varepsilon}|\mathcal{F}^{\mu,T}_{S,\bar{S}}\left[X_{S}\tau_{S}-\tilde{X}_{S}\tau_{\bar{S}}\right],$
(5)

where the $(T,\mu)$-independent factors $X_{S}$ and $\tilde{X}_{S}$ encode
wavefunction overlap and band geometry data. They are given in the SM SM
(Please see Supplemental Material for the microscopic model, additional
details on Boltzmann kinetic theory and derivations of ultrasubsonic transport
theory results, an analysis of the applicability of our theory to MATBLG, and
a treatment of the ultrasubsonic Dirac cone). We may combine Eq. (5) for the
state $S$ with its analogue for the compliment state $\bar{S}$, into a 2x2
matrix equation determining both $\tau(S)$ and $\tau(\bar{S})$, which we solve
analytically for the USS relaxation times:

$\displaystyle\tau_{S}$
$\displaystyle=\frac{2\rho_{M}\hbar^{3}v_{p}^{4}}{D^{2}}\frac{2}{|\varepsilon-\bar{\varepsilon}|^{2}}\sinh\left(\frac{|\varepsilon-\bar{\varepsilon}|}{2k_{B}T}\right)\frac{1}{X_{S}X_{\bar{S}}-\tilde{X}_{S}\tilde{X}_{\bar{S}}}\left[X_{\bar{S}}\frac{\cosh\left(\frac{\bar{\varepsilon}-\mu}{2k_{B}T}\right)}{\cosh\left(\frac{\varepsilon-\mu}{2k_{B}T}\right)}+\tilde{X}_{S}\frac{\cosh\left(\frac{\varepsilon-\mu}{2k_{B}T}\right)}{\cosh\left(\frac{\bar{\varepsilon}-\mu}{2k_{B}T}\right)}\right].$
(6)

Combining the relaxation time formula [Eq. (6)] with the resistivity formula
[Eq. (1)] gives the master formula for the resistivity of the ultra-subsonic
fermion system. We present this in the SM SM (Please see Supplemental Material
for the microscopic model, additional details on Boltzmann kinetic theory and
derivations of ultrasubsonic transport theory results, an analysis of the
applicability of our theory to MATBLG, and a treatment of the ultrasubsonic
Dirac cone). In the case of a particle-hole (PH) symmetric system, the
formulae for the USS resistivity simplifies:

$\displaystyle\frac{\delta^{ij}}{\rho}$
$\displaystyle=\frac{1}{4}\frac{e^{2}}{k_{B}T}\frac{\rho_{M}\hbar^{3}v_{p}^{4}}{D^{2}}\left\lgroup\begin{aligned}
&\frac{1}{\mathcal{A}}\sum_{S}\frac{v^{i}_{S}v^{j}_{S}}{X_{S}-\tilde{X}_{S}}\frac{1}{\varepsilon^{2}}\frac{\sinh\left(\frac{|\varepsilon|}{k_{B}T}\right)}{\cosh\left(\frac{\varepsilon+\mu}{2k_{B}T}\right)\cosh\left(\frac{\varepsilon-\mu}{2k_{B}T}\right)}\\\
+&\frac{1}{2}\sinh\left(\frac{|\mu|}{k_{B}T}\right)^{2}\frac{1}{\mathcal{A}}\sum_{S}\frac{v^{i}_{S}v^{j}_{S}X_{S}}{X^{2}_{S}-\tilde{X}^{2}_{S}}\frac{1}{\varepsilon^{2}}\left[\frac{\sinh\left(\frac{|\varepsilon|}{k_{B}T}\right)}{\cosh\left(\frac{\varepsilon+\mu}{2k_{B}T}\right)\cosh\left(\frac{\varepsilon-\mu}{2k_{B}T}\right)}\right]^{3}\end{aligned}\right\rgroup.$
(7)

Figure 3: We depict the difference in the physics of low-$T$ transport in
subsonic and supersonic fermion systems. In each figure, the band structure of
the fermions is depicted in blue, and the orange filling indicates the
energies of the electrons that contribute meaningfully to transport. The
dotted line denotes the Fermi level. In (A), we show the scenario for
supersonic fermions, in which only the states in the immediate vicinity of the
Fermi surface contribute to transport. In (B), we illustrate the different
physics of subsonic fermion scattering, where the states contributing
meaningfully to low-$T$ transport are those satisfying $|\varepsilon|>|\mu|$.
In this picture, states very far from the Fermi level can contribute much more
than states even just barely below it.

In the limit $k_{B}T\gg|\varepsilon-\bar{\varepsilon}|$, expanding the $\sinh$
factor in Eq. (6) gives the expected $\tau\propto T^{-1}$ scattering rate of
the equipartition regime. However, when
$k_{B}T\ll|\varepsilon-\bar{\varepsilon}|$, instead of the usual crossover to
a $\tau\propto T^{-4}$ power law, we find an exponential blow up of the
relaxation times, in line with the physical picture in which the entire
scattering manifold becomes thermally inaccessible.

Extreme low-$T$ limit.—In the PH symmetric case [Eq. 7], the only
$T$-dependent factor in the Brillouin zone summations is is the kernel

$\displaystyle
K(\varepsilon,\mu,T)\equiv\frac{\sinh\left(\frac{|\varepsilon|}{k_{B}T}\right)}{\cosh\left(\frac{\varepsilon+\mu}{2k_{B}T}\right)\cosh\left(\frac{\varepsilon-\mu}{2k_{B}T}\right)}.$
(8)

This factor represents the competition at low-$T$ between the divergence of
the relaxation times and the thermodynamic suppression of states far from the
Fermi level. In the extreme low-$T$ limit, where all the $\sinh$ and $\cosh$
functions blow up exponentially, we simply have

$\displaystyle K(\varepsilon,\mu,T\rightarrow 0)\rightarrow
2\Theta[|\varepsilon|-|\mu|].$ (9)

The thermal kernel $K$ determines which Bloch states meaningfully contribute
to transport. Equation (9) suggests that all states with $|\mu|<|\varepsilon|$
contribute equally, while states with $|\mu|>|\varepsilon|$ do not contribute,
even though they may be close to the Fermi energy. We emphasize that this is
radically different than the usual paradigm, where low-$T$ transport is almost
entirely determined by states in the immediate vicinity of the Fermi level.
This is depicted in Fig. 3.

Applying Eq. (9), the two summations in Eq. (7) are simply $T$-independent
constants ($\mathcal{C}_{1}$, $\mathcal{C}_{2}$), and the extreme low-$T$
expression for the resistivity is

$\displaystyle\rho$ $\displaystyle\approx\frac{\hbar
D^{2}}{4e^{2}\rho_{M}(\hbar
v_{p})^{4}}\frac{k_{B}T}{\mathcal{C}_{1}+\mathcal{C}_{2}\sinh[|\mu|/(k_{B}T)]^{2}}.$
(10)

From Eq. (10), it is apparent that at the charge neutrality point ($\mu=0$),
we have purely linear-in-$T$ phonon-induced resistivity down to $T=0$. On the
other hand, when $\mu\neq 0$, then the low-$T$ resistivity is proportional to
the factor $\rho\propto\exp[-2|\mu|/(k_{B}T)]$, and is exponentially
suppressed when $k_{B}T\ll 2|\mu|$, giving a crossover to perfect
conductivity. Comparing this with the physics of supersonic fermions, we find

$\displaystyle
k_{B}T_{\text{crossover}}=2|\mu|=\frac{v_{F}}{v_{p}}k_{B}T^{\text{trad}}_{BG},$
(11)

where $T^{\text{trad}}_{BG}\equiv 2v_{p}k_{F}$ gives the traditional lower
bound for the regime of linear-in-$T$ resistivity based on the usual Bloch-
Grüneisen paradigm. The crossover temperature is parametrically suppressed by
the small parameter, $v_{F}/v_{p}$. The subsonic case is thus expected to host
linear-in-$T$ resistivity scaling down to a significantly lower temperature
than one would estimate based on intuition from supersonic fermion scattering.

Concluding discussion.—Phonon scattering processes in subsonic fermion systems
are necessarily interband transitions, implying the existence of a nonzero,
minimum energy transfer allowed by kinematics. At temperatures below this
energy scale, phonon scattering is suppressed and Bloch states become long-
lived, with important effects on the low-$T$ transport physics. Applying
standard Boltzmann kinetic theory of acoustic phonon scattering to nearly flat
band (NFB) systems, we find that this underlies a robust, linear-in-$T$
scaling of the resistivity down to temperatures far lower than the Bloch-
Grüneisen paradigm of supersonic fermion scattering would suggest. This result
provides a concrete mechanism for the linear-in-$T$ resistivity over a wide
range of temperatures based only on familiar concepts of solid state physics
and the distinct kinematics of NFB systems. In particular, it provides an
alternative theoretical explanation for low-$T$ “strange metal” resistivity
scaling in NFB systems. Our theory also predicts regimes of perfect
conductivity and $\rho\propto T^{2}$ scaling which compete with the
$\rho\propto T$ regime, which are not commonly associated with phonon physics.
Both $\rho\propto T^{2}$ and low-$T$ $\rho\propto T$ behaviors are reported in
MATBLG Jaoui et al. (2022). Thus, low-$T$ $\rho\propto T$, $\rho\propto
T^{2}$, and perfect conductivity regimes may all arise in NFB systems from the
same universal phonon-scattering physics, providing a possible explanation for
these reported observations in MATBLG.

Our theory of USS fermion transport is general and could apply widely to NFB
systems in other 2D heterostructures or heavy fermion materials. However, our
results are particularly germane to the ongoing debate on the presence of a
strange metal state in MATBLG. Using the Bistritzer-MacDonald (BM) model for
the non-interacting band structure of MATBLG, we see that much of the band is
made up of subsonic states ($\approx 92\%$). Large sections of the NFBs have
small band separation ($\leq 0.0005\ eV$) and small $v_{F}/v_{p}$ ($\leq
0.02$), rendering USS theory highly applicable. We discuss this more
quantitatively in the SM SM (Please see Supplemental Material for the
microscopic model, additional details on Boltzmann kinetic theory and
derivations of ultrasubsonic transport theory results, an analysis of the
applicability of our theory to MATBLG, and a treatment of the ultrasubsonic
Dirac cone). We emphasize that a direct quantitative application of USS theory
to the BM model could be of limited utility since the bands are expected to be
renormalized by interaction effects, strain, and twist disorder.

###### Acknowledgements.

We thank Fengcheng Wu for helpful discussions. This work is supported by the
Laboratory for Physical Sciences.

## References

  * Geim and Grigorieva (2013) A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013), URL https://doi.org/10.1038%2Fnature12385.
  * Novoselov et al. (2006) K. S. Novoselov, E. McCann, S. V. Morozov, V. I. Fal’ko, M. I. Katsnelson, U. Zeitler, D. Jiang, F. Schedin, and A. K. Geim, Nature Physics 2, 177 (2006), URL https://doi.org/10.1038%2Fnphys245.
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011), URL https://doi.org/10.1073%2Fpnas.1108174108.
  * Morell et al. (2010) E. S. Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Physical Review B 82 (2010), URL https://doi.org/10.1103%2Fphysrevb.82.121407.
  * Li et al. (2019) X. Li, F. Wu, and A. H. MacDonald, _Electronic structure of single-twist trilayer graphene_ (2019), URL https://arxiv.org/abs/1907.12338.
  * Kim et al. (2017) K. Kim, A. DaSilva, S. Huang, B. Fallahazad, S. Larentis, T. Taniguchi, K. Watanabe, B. J. LeRoy, A. H. MacDonald, and E. Tutuc, Proceedings of the National Academy of Sciences 114, 3364 (2017), URL https://doi.org/10.1073%2Fpnas.1620140114.
  * Cao et al. (2018a) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, et al., Nature 556, 80 (2018a), URL https://doi.org/10.1038%2Fnature26154.
  * Cao et al. (2018b) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b), URL https://doi.org/10.1038%2Fnature26160.
  * Cao et al. (2020a) Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Phys. Rev. Lett. 124, 076801 (2020a), URL https://link.aps.org/doi/10.1103/PhysRevLett.124.076801.
  * Cao et al. (2021) Y. Cao, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 595, 526 (2021).
  * Yankowitz et al. (2019) M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019), URL https://doi.org/10.1126%2Fscience.aav1910.
  * Kerelsky et al. (2019) A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, et al., Nature 572, 95 (2019), URL https://doi.org/10.1038%2Fs41586-019-1431-9.
  * Lu et al. (2019) X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, et al., Nature 574, 653 (2019), URL https://doi.org/10.1038%2Fs41586-019-1695-0.
  * Stepanov et al. (2020) P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe, T. Taniguchi, F. H. Koppens, J. Lischner, L. Levitov, and D. K. Efetov, Nature 583, 375 (2020).
  * Sharpe et al. (2019) A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019), URL https://doi.org/10.1126%2Fscience.aaw3780.
  * Chen et al. (2020) G. Chen, A. L. Sharpe, E. J. Fox, Y.-H. Zhang, S. Wang, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, et al., Nature 579, 56 (2020), URL https://doi.org/10.1038%2Fs41586-020-2049-7.
  * Rozen et al. (2021) A. Rozen, J. M. Park, U. Zondiner, Y. Cao, D. Rodan-Legrain, T. Taniguchi, K. Watanabe, Y. Oreg, A. Stern, E. Berg, et al., Nature 592, 214 (2021).
  * Zhou et al. (2022) H. Zhou, L. Holleis, Y. Saito, L. Cohen, W. Huynh, C. L. Patterson, F. Yang, T. Taniguchi, K. Watanabe, and A. F. Young, Science 375, 774 (2022), URL https://doi.org/10.1126%2Fscience.abm8386.
  * Zhou et al. (2021a) H. Zhou, T. Xie, T. Taniguchi, K. Watanabe, and A. F. Young, Nature 598, 434 (2021a), URL https://doi.org/10.1038%2Fs41586-021-03926-0.
  * Zhou et al. (2021b) H. Zhou, T. Xie, A. Ghazaryan, T. Holder, J. R. Ehrets, E. M. Spanton, T. Taniguchi, K. Watanabe, E. Berg, M. Serbyn, et al., Nature 598, 429 (2021b), URL https://doi.org/10.1038%2Fs41586-021-03938-w.
  * Serlin et al. (2020) M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. F. Young, Science 367, 900 (2020), URL https://doi.org/10.1126%2Fscience.aay5533.
  * Wu et al. (2018) F. Wu, A. MacDonald, and I. Martin, Physical Review Letters 121 (2018), URL https://doi.org/10.1103%2Fphysrevlett.121.257001.
  * Wu et al. (2019a) F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. H. MacDonald, Phys. Rev. Lett. 122, 086402 (2019a), URL https://link.aps.org/doi/10.1103/PhysRevLett.122.086402.
  * Tschirhart et al. (2022) C. L. Tschirhart, E. Redekop, L. Li, T. Li, S. Jiang, T. Arp, O. Sheekey, T. Taniguchi, K. Watanabe, K. F. Mak, et al., _Intrinsic spin hall torque in a moire chern magnet_ (2022), URL https://arxiv.org/abs/2205.02823.
  * Polshyn et al. (2020) H. Polshyn, J. Zhu, M. A. Kumar, Y. Zhang, F. Yang, C. L. Tschirhart, M. Serlin, K. Watanabe, T. Taniguchi, A. H. MacDonald, et al., Nature 588, 66 (2020), URL https://doi.org/10.1038%2Fs41586-020-2963-8.
  * Jaoui et al. (2021) A. Jaoui, I. Das, G. Di Battista, J. Díez-Mérida, X. Lu, K. Watanabe, T. Taniguchi, H. Ishizuka, L. Levitov, and D. K. Efetov, _Quantum critical behavior in magic-angle twisted bilayer graphene_ (2021), URL https://arxiv.org/abs/2108.07753.
  * Polshyn et al. (2019) H. Polshyn, M. Yankowitz, S. Chen, Y. Zhang, K. Watanabe, T. Taniguchi, C. R. Dean, and A. F. Young, Nature Physics 15, 1011 (2019).
  * Cao et al. (2020b) Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Phys. Rev. Lett. 124, 076801 (2020b), URL https://link.aps.org/doi/10.1103/PhysRevLett.124.076801.
  * Sarma and Wu (2022) S. D. Sarma and F. Wu, Physical Review Research 4 (2022), URL https://doi.org/10.1103%2Fphysrevresearch.4.033061.
  * Zhang et al. (2022) Y. Zhang, R. Polski, A. Thomson, . Lantagne-Hurtubise, C. Lewandowski, H. Zhou, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, _Spin-orbit enhanced superconductivity in bernal bilayer graphene_ (2022), URL https://arxiv.org/abs/2205.05087.
  * Polski et al. (2022) R. Polski, Y. Zhang, Y. Peng, H. S. Arora, Y. Choi, H. Kim, K. Watanabe, T. Taniguchi, G. Refael, F. von Oppen, et al., _Hierarchy of symmetry breaking correlated phases in twisted bilayer graphene_ (2022), URL https://arxiv.org/abs/2205.05225.
  * Xie and MacDonald (2020) M. Xie and A. MacDonald, Physical Review Letters 124 (2020), URL https://doi.org/10.1103%2Fphysrevlett.124.097601.
  * Andrei and MacDonald (2020) E. Y. Andrei and A. H. MacDonald, Nature Materials 19, 1265 (2020), URL https://doi.org/10.1038%2Fs41563-020-00840-0.
  * Li et al. (2021) T. Li, S. Jiang, L. Li, Y. Zhang, K. Kang, J. Zhu, K. Watanabe, T. Taniguchi, D. Chowdhury, L. Fu, et al., Nature 597, 350 (2021), URL https://doi.org/10.1038%2Fs41586-021-03853-0.
  * Ghiotto et al. (2021) A. Ghiotto, E.-M. Shih, G. S. S. G. Pereira, D. A. Rhodes, B. Kim, J. Zang, A. J. Millis, K. Watanabe, T. Taniguchi, J. C. Hone, et al., Nature 597, 345 (2021), URL https://doi.org/10.1038%2Fs41586-021-03815-6.
  * Pan et al. (2020) H. Pan, F. Wu, and S. D. Sarma, Physical Review B 102 (2020), URL https://doi.org/10.1103%2Fphysrevb.102.201104.
  * Pan and Sarma (2021) H. Pan and S. D. Sarma, Physical Review Letters 127 (2021), URL https://doi.org/10.1103%2Fphysrevlett.127.096802.
  * Morales-Durán et al. (2021) N. Morales-Durán, A. H. MacDonald, and P. Potasz, Physical Review B 103 (2021), URL https://doi.org/10.1103%2Fphysrevb.103.l241110.
  * Ahn and Sarma (2022) S. Ahn and S. D. Sarma, Physical Review B 105 (2022), URL https://doi.org/10.1103%2Fphysrevb.105.115114.
  * Kerelsky et al. (2021) A. Kerelsky, C. Rubio-Verdú, L. Xian, D. M. Kennes, D. Halbertal, N. Finney, L. Song, S. Turkel, L. Wang, K. Watanabe, et al., Proceedings of the National Academy of Sciences 118 (2021).
  * Khalaf et al. (2019) E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath, Phys. Rev. B 100, 085109 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.085109.
  * Arora et al. (2020) H. S. Arora, R. Polski, Y. Zhang, A. Thomson, Y. Choi, H. Kim, Z. Lin, I. Z. Wilson, X. Xu, J.-H. Chu, et al., Nature 583, 379 (2020), URL https://doi.org/10.1038%2Fs41586-020-2473-8.
  * Jaoui et al. (2022) A. Jaoui, I. Das, G. D. Battista, J. Díez-Mérida, X. Lu, K. Watanabe, T. Taniguchi, H. Ishizuka, L. Levitov, and D. K. Efetov, Nature Physics 18, 633 (2022), URL https://doi.org/10.1038%2Fs41567-022-01556-5.
  * Wu et al. (2019b) F. Wu, E. Hwang, and S. Das Sarma, Phys. Rev. B 99, 165112 (2019b), URL https://link.aps.org/doi/10.1103/PhysRevB.99.165112.
  * Li et al. (2020) X. Li, F. Wu, and S. Das Sarma, Phys. Rev. B 101, 245436 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.101.245436.
  * Davis et al. (2023a) S. M. Davis, F. Wu, and S. D. Sarma, _Acoustic phonon contribution to the resistivity of twisted bilayer graphene_ (2023a), eprint 2302.07275.
  * Saito et al. (2020) Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, and A. F. Young, Nature Physics 16, 926 (2020).
  * Hwang and Das Sarma (2008) E. H. Hwang and S. Das Sarma, Phys. Rev. B 77, 115449 (2008), URL https://link.aps.org/doi/10.1103/PhysRevB.77.115449.
  * Davis et al. (2023b) S. M. Davis, Y.-Z. Chou, F. Wu, and S. Das Sarma, Phys. Rev. B 107, 045426 (2023b), URL https://link.aps.org/doi/10.1103/PhysRevB.107.045426.
  * Hwang and Sarma (2019) E. H. Hwang and S. D. Sarma, Physical Review B 99 (2019), URL https://doi.org/10.1103%2Fphysrevb.99.085105.
  * SM (Please see Supplemental Material for the microscopic model, additional details on Boltzmann kinetic theory and derivations of ultrasubsonic transport theory results, an analysis of the applicability of our theory to MATBLG, and a treatment of the ultrasubsonic Dirac cone) (Please see Supplemental Material for the microscopic model, additional details on Boltzmann kinetic theory and derivations of ultrasubsonic transport theory results, an analysis of the applicability of our theory to MATBLG, and a treatment of the ultrasubsonic Dirac cone).
  * Ashcroft and Mermin (1976) N. W. Ashcroft and N. D. Mermin, _Solid State Physics_ (Harcourt College Publishers, 1976), ISBN 9780030839931\.
  * Ziman (1960) J. M. Ziman, _Electrons and Phonons_ (Oxford University Press, 1960), ISBN 9780198507796.
  * Coleman (2015) P. Coleman, _Introduction to Many-Body Physics_ (Cambridge University Press, 2015), ISBN 9780521864886\.


# DeepGraphONet: A Deep Graph Operator Network to Learn and Zero-shot Transfer
the Dynamic Response of Networked Systems

Yixuan Sun, Christian Moya, Guang Lin, Meng Yue
Y. Sun and G. Lin are with the School of Mechanical Engineering, Purdue
University, West Lafayette, IN. E-mail: {yixuan-sun<EMAIL_ADDRESS>Moya and G. Lin are with the Department of Mathematics, Purdue University,
West Lafayette, IN. e-mail: {cmoyacal<EMAIL_ADDRESS>Yue is with the
Interdisciplinary Science Department, Brookhaven National Laboratory, Upton,
NY. e-mail<EMAIL_ADDRESS>
###### Abstract

This paper develops a Deep Graph Operator Network (DeepGraphONet) framework
that learns to approximate the dynamics of a complex system (_e.g.,_ the power
grid or traffic) with an underlying sub-graph structure. We build our
DeepGraphONet by fusing the ability of (i) Graph Neural Networks (GNN) to
exploit spatially correlated graph information and (ii) Deep Operator Networks
(DeepONet) to approximate the solution operator of dynamical systems. The
resulting DeepGraphONet can then predict the dynamics within a given
short/medium-term time horizon by observing a finite history of the graph
state information. Furthermore, we design our DeepGraphONet to be resolution-
independent. That is, we do not require the finite history to be collected at
the exact/same resolution. In addition, to disseminate the results from a
trained DeepGraphONet, we design a zero-shot learning strategy that enables
using it on a different sub-graph. Finally, empirical results on the (i)
transient stability prediction problem of power grids and (ii) traffic flow
forecasting problem of a vehicular system illustrate the effectiveness of the
proposed DeepGraphONet.

###### Index Terms:

Deep learning, operator regression, graph neural networks, networked dynamical
systems

##  I Introduction

Networked dynamical systems are ubiquitous in science and engineering, _e.g.,_
the power grid or traffic networks. To simulate and predict the dynamics of
such systems, researchers have developed sophisticated, high-fidelity
numerical schemes that can accurately solve the corresponding governing
equations. However, for tasks requiring multiple forward simulations, _e.g.,_
optimization, uncertainty quantification, and control, these high-fidelity
schemes may become prohibitively expensive [1].

Deep learning techniques can potentially address this computational cost
problem by acting as a fast surrogate model trained using data. As a result,
many deep learning techniques have been proposed in the literature to simulate
and predict complex dynamical systems. For example, one can use traditional
neural networks, such as fully connected neural networks or recurrent neural
networks (RNN), for time-dependent prediction tasks [2, 3, 4]. Such networks
have the goal of learning (i) the underlying governing equations [5, 6, 7] or
(ii) the future (_i.e.,_ the next state) dynamic response [8, 9]. However,
such traditional neural networks often require a massive amount of data (with
a fixed resolution) to learn the system’s dynamic response to a particular
operating condition. For other operating conditions (_e.g.,_ different initial
conditions), one often needs to retrain the model with a different dataset.
Such a limitation has prevented the wider use of traditional neural networks
for solving complex dynamical systems in science and engineering, where the
data is scarce and expensive to collect.

To address such a limitation, several works have proposed to learn the
solution operator (_i.e.,_ a mapping between infinite-dimensional spaces) of
complex dynamical systems using, for example, Deep Operator Networks
(DeepONet) [10], Graph Neural Operators (GNO) [11], or Fourier Neural
Operators (FNO) [12]. In particular, the DeepONet, introduced in the seminal
paper [10] and developed based on the universal approximation theorem for
nonlinear operators [13], has demonstrated remarkable accuracy and
generalization capability for learning the solution operator of non-autonomous
systems.

We used the DeepONet framework in a recent study [14] to learn the dynamic
response of the power grid after a disturbance without using recurrence and
fixed resolution, requirements typically adopted by RNNs. However, our work
only considered the dynamic response at the bus level. Thus, it neglects the
possibly rich information from the interaction between a bus and its
neighbors. Such a spatial correlation is critical for learning the solution
operator of networked dynamical systems, _e.g.,_ the power grid or traffic
networks.

To enable learning the solution operator of complex networked dynamical
systems, we propose, in this paper, the Deep Graph Operator Network
(DeepGraphONet) framework. DeepGraphONet is a muti-input multi-output Deep
Operator Network that learns the dynamic response of complex networked systems
by exploiting spatially correlated information from the given underlying graph
or sub-graph structure.

###  I-A Related Work

Learning dynamical systems. Many works have proposed using machine and deep
learning to learn unknown dynamical systems from time-series data. In
particular, we classify such works into learning the system’s (i) governing
equations [5, 7, 15] and (ii) future response [8, 9]. For example, in [5], the
authors use a dictionary of functions to learn a sparse representation of the
system’s governing equations. Zhang _et al._ [15] adopted Bayesian sparse
regression to identify differential equation terms from a large pool of
candidates with error bars.

On the other hand, Qin  _et al._ [8] trains a neural network to learn the next
state response of the system given the current state. The framework then can
predict the system’s future response by recursively using the trained network.
In [9], the authors proposed a multistep method with a feed-forward neural
network to approximate the dynamical system response. Most of the above works
can effectively learn the dynamical system for a single operating condition.
However, the above works will require a prohibitive amount of data and
training resources to learn the system’s response to many operating
conditions. To alleviate such limitations, we will use the novel framework of
deep operator learning in this work.

Deep operator learning. Traditional deep learning techniques [16] focus on
approximating the mapping between Euclidean spaces. However, these traditional
techniques may not be adequate for learning the solution operator of a complex
dynamical system. To learn such an operator, several works have proposed
designing operator learning frameworks based on deep neural networks [10, 11,
12, 17].

We build our work on the Deep Operator Network (DeepONet) framework introduced
in [10]. DeepONet is a neural network architecture that can learn nonlinear
operators by using (i) a Branch network to process input information and (ii)
a Trunk network to process query locations within the output function domain.
Then, one computes the output function at a query location by merging the
features from both the Branch and Trunk nets using a dot product. DeepONet has
demonstrated its exceptional approximation and generalization capabilities
using a low amount of data for problems in power engineering [14],
electroconvection, and multi-physics tasks [18], or material science tasks
[19]. DeepONet, however, is only a single-input single-output operator
framework. Thus, one cannot directly use DeepONet for networked dynamical
systems.

In [20], the authors proposed MIONet, a multi-input function single-output
function DeepONet to alleviate the multiple-input problem. Nevertheless, the
MIONet framework fails to consider the spatial correlation among input
functions of a networked system and only produces a single output function.
Our work effectively alleviates both limitations by incorporating a Graph
Neural Network within the DeepONet framework.

Graph Neural Networks. In traditional deep learning, the training data,
_e.g.,_ time series, tabular-like data, and images, is well-structured, and
one assumes it belongs to the Euclidean space. However, we cannot naturally
assume that other forms of data, such as power grid and traffic data, belong
to the Euclidean space. As a result, researchers developed Graph Neural
Networks (GNN) [21] to capture the spatial correlation of data with an
underlying graph representation.

At the early stage, GNNs were used within recursive schemes [22, 23] to reach
steady node states for subsequent tasks. These schemes, however, suffered from
high-computational costs and were limited to some specific graph structures.
To alleviate such issues, other works [24, 25, 26] build on graph signal
processing to develop graph spectral convolution networks and their localized
variants. In graph spectral convolution, one performs the graph convolution
operation in the Fourier domain. Such spectral convolution requires the
eigendecomposition of the graph Laplacian. Spectral convolution, however, is
expensive and does not allow transferring to another graph with a different
topology [27, 28]. This limitation prevents the trained GNNs from achieving
zero-shot learning.

On the other hand, similar to standard convolutional neural networks, several
works [29, 30] proposed non-spectral and spatial convolution directly defined
on the graph. In particular, the spatial convolution locally operates on
nodes, and, thus, one can use it on different graphs. Furthermore, by stacking
spatial convolutional layers, we can effectively process information from
nodes located more than one hop away.

To learn with graph data efficiently and effectively, the authors in [31, 32,
33] employed the spatial graph convolution within a message passing framework.
The main idea of such a framework is to update the node state information by
aggregating information from neighboring nodes, followed by a neural network-
based transformation. In this work, we adopt such a message passing framework
and, in particular, the framework from GraphSage [33], to process the input
information to the Deep Operator Network.

This work aims at learning the dynamic evolution of graph signals using a
fixed graph structure. Many studies [34, 35, 36, 37, 38] have combined deep
learning techniques for time series (_e.g.,_ RNNs [36, 35] and CNNs [34, 37,
38]) with GNNs to produce time-dependent predictions that can take into
account spatial correlations. However, such methods may (i) lack the
flexibility and (ii) require high computational-cost to learn the solution
operator of a networked dynamical system. Thus, in our work, we will employ
the deep operator learning framework to learn the parametric mapping between
graph functions directly.

###  I-B Our Work

The objectives of this work are two-fold.

  1. 1.

Approximating the solution operator. We aim to derive a deep learning-based
method to approximate the solution operator of a system in which the dynamics
evolve within an underlying sub-graph structure. For instance, we aim to use
our method to learn the dynamics of a control area within a large-scale power
grid or the traffic dynamics of a city’s neighborhood.

  2. 2.

Zero-shot learning. We aim to apply the trained proposed model directly to
approximate the solution operator on unseen graphs or sub-graphs with
different structures. In other words, we aim to achieve zero-shot learning.
For instance, we aim to train the model on a small sub-graph and then use it
to make predictions on a larger graph with high accuracy.

We detail the contributions of this work next.

  * •

We first build (in Section III) a Deep Graph Operator Network (DeepGraphONet)
framework that learns to approximate the solution operator of a dynamical
system with an underlying connected sub-graph structure. DeepGraphONet fuses
the ability of (i) Graph Neural Networks (GNN) to exploit the graph
information and (ii) DeepONets [10] to approximate the solution operator of
nonlinear systems. The proposed DeepGraphONet takes as inputs (i) a finite
history of the graph state information and (ii) the desired query location
with an arbitrary resolution for short/medium term prediction. Compared to the
vanilla DeepONet [10], the proposed DeepGraphONet is (i) a multi-input multi-
output framework and (ii) resolution-independent in the input function; that
is, we do not require the Branch sensors to be fixed.

  * •

We then propose (in Section IV) a zero-shot learning strategy that exploits
the property of message passing GNNs to enable directly using a trained
DeepGraphONet on a different graph or sub-graph for the same task.

  * •

Finally, we verify the efficacy of the DeepGraphONet on (i) the transient
stability prediction problem of the IEEE 16-machine 68-bus system and (ii) the
traffic dynamics forecasting problem using the METR-LA dataset, which contains
traffic information measured using loop detectors in the highway of Los
Angeles County [39].

We organize the rest of this paper as follows. Section II introduces the
problem of approximating the solution operator of a dynamical system with an
underlying sub-graph structure. The DeepGraphONet framework to approximate
such a solution operator is detailed in Section III-A. We present the
resolution-independent DeepGraphONet and our zero-shot learning strategy in
Section III-B and Section III-C, respectively. Numerical experiments are used
to illustrate the efficacy of DeepGraphONet in Section IV. Finally, Section V
discusses our future work and Section VI concludes the paper.

##  II Problem Settings

We consider a complex networked dynamical system. We model its (i) networked
structure using the undirected graph $G=(V,E)$, where $V$ is the set of $|V|$
nodes and $E$ is the set of $|E|$ edges, and (ii) dynamics using the initial
value problem (IVP):

$\displaystyle\begin{aligned} \frac{d}{dt}\mathbf{x}(t)&=f(\mathbf{x}(t);G)\\\
\mathbf{x}(t_{0})&=\mathbf{x}_{0}\end{aligned}$ (1)

Here $\mathbf{x}(t)\in\mathcal{X}\subset\mathbb{R}^{d\times|V|}$ is the graph-
valued state function, $\mathbf{x}_{0}\in\mathcal{X}$ the initial condition,
and $f:\mathcal{X}\to\mathcal{X}$ the unknown vector field. For simplicity, we
assume each node’s state has the same dimension $d=1$ and, with some abuse of
notation, we make explicit the dependency of the dynamics $f$ on the
underlying graph structure $G$. We also let for future work the case when we
know an approximate model of the dynamics $f_{\text{approx}}\approx f$.

Solution operator. Our goal is to approximate the dynamic response of (1)
described via the solution operator (also known as flow map):

$\displaystyle\mathcal{F}(\mathbf{x}(t_{0});G)(t)$
$\displaystyle\equiv\mathbf{x}(t)=\mathbf{x}(t_{0})+\int^{t}_{t_{0}}f(\mathbf{x}(\beta);G)d\beta,$

for all $t\in[t_{0},t_{f}]$, where $[t_{0},t_{f}]\subset[t_{0},\infty)$
denotes a finite time interval where the solution operator exists and is
unique.

For large-scale networked dynamical systems, _i.e.,_ when $|V|\gg 1$, learning
the above solution operator may be prohibitively expensive. Moreover, one may
not have the access to all graph states and may be interested in learning the
dynamics of a sub-region within the network, _i.e.,_ a sub-graph. For example,
for power grids, we may want to learn the dynamic response of a Control Area,
or, for traffic networks, we may want to approximate the dynamic response of a
city or a district.

To this end, let us split the graph-valued state function as
$\mathbf{x}=\left(\mathbf{x}_{S},\mathbf{x}_{S^{c}}\right)$. Here,
$S=(V_{S},E_{S})$ denotes the sub-graph region of interest within the network,
where $V_{S}\subset V$ is a connected set of nodes within $V$ and
$E_{s}\subset E$, and $S^{c}$ denotes the complement of $S$. To describe the
dynamics of $\mathbf{x}_{S}(t)$, we employ the Mori-Zwanzig formalism [40,
41], which yields

$\displaystyle\frac{d}{dt}\mathbf{x}_{S}(t)=R(\mathbf{x}_{S}(t);S)+\int_{t_{0}}^{t}\mathcal{K}(\mathbf{x}_{S}(t-\beta),\beta)d\beta+\mathcal{N}(\mathbf{x}_{0})$

In the above, the first term $R$ is the Markovian term, which depends on the
current value of $\mathbf{x}_{S}$; the second term is the memory integral,
which depends on the values of the state and input variables from the initial
time $\beta=t_{0}$ to the current time $\beta=t$. This memory integral
involves $\mathcal{K}$, which is commonly known as the memory kernel. Finally,
the third term is the “noise” term, which depends on the initial condition
$\mathbf{x}_{0}$.

Decaying memory assumption. Learning the memory integral from the initial time
$\beta=t_{0}$ is a challenging task. To alleviate such a challenge, we assume
the memory decays over time, that is, there exists $t_{M}\in(t_{0},t)$ (we
treat $t_{M}$ as a hyperparameter in this paper) such that we can neglect the
effect of the memory for all $\beta\in[t_{0},t-t_{M})$. Formally, this
corresponds to the following approximation of the memory integral. For any
given $\epsilon>0$, there exists $t_{0}<t_{M}<t$ such that:

$\displaystyle\left\lvert\int_{t_{0}}^{t}\mathcal{K}(\mathbf{x}_{S}(t-\beta),\beta)d\beta-\int_{t_{0}}^{t_{M}}\mathcal{K}(\mathbf{x}_{S}(t-\beta),\beta)d\beta\right\rvert<\epsilon$
(2)

Discrete representation of the integral. In practice, one computes the
truncated integral using some discretized scheme (_e.g.,_ quadrature) such
that for any given $\kappa>0$, we have

$\displaystyle\begin{split}&\left|\int_{t_{0}}^{t_{M}}\mathcal{K}(\mathbf{x}_{S}(t-\beta),\beta)d\beta\right.\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \qquad\leavevmode\nobreak\
\qquad\left.-Q(\mathbf{x}_{S}(t-\tau_{1}),\dots,\mathbf{x}_{S}(t-\tau_{m}))\vphantom{\int_{t_{0}}^{t_{M}}}\right|<\kappa,\end{split}$
(3)

where the memory partition of size $m$ (we treat this memory resolution $m$ as
a hyperparameter in this paper) is arbitrary and satisfies
$0\leq\tau_{1},\ldots,\tau_{m}\leq t_{M}$ and
$\\{t-\tau_{i}\\}_{i=1}^{m}\subset[t-t_{M},t]$.

These decaying memory and discrete representation assumptions lead to the
following approximate dynamics of $\mathbf{x}_{S}$:

$\displaystyle\frac{d}{dt}\mathbf{x}_{S}(t)=R(\mathbf{x}_{S}(t))+Q(\mathbf{x}_{S}(t-\tau_{1}),\dots,\mathbf{x}_{S}(t-\tau_{m})).$
(4)

In this paper, we employ the local solution operator within the arbitrary
interval $[t,t+h_{n}]$, $h_{n}\in[0,h]$, where $h$ is the predicting horizon,
and we treat it as a hyperparameter in this paper. For all $t>t_{0}$, the
local memory sensor locations satisfy
$\\{t-\tau_{i}\\}^{m}_{i=1}\subset[t-t_{M},t]$. With this local memory, the
local solution operator is:

$\displaystyle\begin{split}\tilde{\mathcal{F}}(\mathbf{x}_{S}(t),&\\{\mathbf{x}_{S}(t-\tau_{i})\\}^{m}_{i=1};S)(h_{n})=\mathbf{x}_{S}(t)\\\
&+\int_{t}^{t+h_{n}}R(\mathbf{x}_{S}(\beta);S)d\beta\\\ &\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+Q(\mathbf{x}_{S}(t-\tau_{1}),\dots,\mathbf{x}_{S}(t-\tau_{m})).\end{split}$

In the above, we emphasized that the solution operator depends on the finite
history of the sub-graph state information $\mathbf{x}_{S}(t)$ within the time
interval $[t-t_{M},t]$.

Our goal in this paper is to design a Deep Graph Operator Network
(DeepGraphONet) $\mathcal{F}_{\theta}$, with trainable parameters
$\theta\in\mathbb{R}^{p}$, to approximate the solution operator
$\tilde{\mathcal{F}}$ for all $h_{n}\in[0,h]$. To control the approximation
power of the proposed DeepGraphONet, we will use two hyperparameters: the
memory’s (i) size $t_{M}$ and (ii) predicting horizon $h$.

##  III Methods

This section describes the proposed Deep Graph Operator Network
(DeepGraphONet) $\mathcal{F}_{\theta}$ to approximate the solution operator
$\tilde{\mathcal{F}}$.

Fig. 1: The DeepGraphONet architecture. The Branch Net (i) is a Graph Neural
Network that uses message passing to learn the nodes’ representation on the
networked system’s graph $S$ and (ii) takes as input the graph state signals
$\\{\mathbf{x}_{S}(t-\tau_{i})\\}_{i=1}^{m}$, $\tau\in[0,t_{M}]$. The Trunk
Net is a feed-forward neural network that takes as input the query location
$h_{n}$ within the prediction horizon $[0,h]$. We obtain the DeepGraphONet’s
output
$\mathcal{F}_{\theta}(\\{\mathbf{x}_{S}(t-\tau_{i})\\}_{i=1}^{m};S)(h_{n})$ by
applying the GNN-based dot product between the GNN-Branch Net and the Trunk
Net. The GNN-based dot product performs the standard dot product between each
node in the Branch Net-produced graph and the Trunk Net.

###  III-A The Deep Graph Operator Network

Building on the DeepONet, introduced in [10], we propose DeepGraphONet
$\mathcal{F}_{\theta}$, a multi-input multi-output Deep Graph Operator
Network, which consists of two neural networks (see Fig. 1): the Branch
network and the Trunk network.

The Branch network. The Branch net is a Graph Neural Network (GNN) that
processes the finite history of the graph-state information. The Branch maps
the graph-state information $\\{\mathbf{x}_{S}(t-\tau_{i})\\}_{i=1}^{m}$,
where $\mathbf{x}_{S}(t-\tau_{i})\in\mathbb{R}^{|S|}$, to the graph-based
coefficients $\mathbf{b}\in\mathbb{R}^{q\times|S|}$.

To build the proposed GNN-based Branch network, we use a collection of message
passing convolutional graph layers, which we adopted from [33]. For all $i\in
S$, the forward pass of the corresponding message passing layer reads:

$\displaystyle\begin{split}x_{i}^{\prime}(t)=\sigma(W_{1}x_{i}(t)+W_{2}\cdot\text{mean}_{j\in\mathcal{N}_{i}}(x_{j}(t))),\end{split}$

where $x_{i}(t)$ is the $i$th node’s state at time $t$, $\mathcal{N}_{i}$ is
the set of nodes adjacent to the $i$th node, $W_{1}$, $W_{2}$ are trainable
weights, and $\sigma$ is a non-linear activation function.

By stacking these message passing layers, we can transport information from
nodes located two or more hops away to the $i$th node. Such a process allows
the network to learn automatically spatial dependencies among nodes.
Furthermore, since our goal is to obtain the dynamic response of all nodes, we
do not use a readout (pooling) in the proposed GNN-based Branch net.

######  Remark.

Fixed v.s. resolution-independent sensors. If the sensor locations used to
discretize the input representing the finite history of graph-state
information $\\{\mathbf{x}_{S}(t-\tau_{i})\\}_{i=1}^{m}$ are fixed (as in
[10]), then we denote the proposed framework as the standard DeepGraphONet.
If, on the other hand, the sensor locations can change over time, then we
denote the proposed framework as the resolution-independent DeepGraphONet (see
Section III-B for more details).

The Trunk network. The Trunk net is a multi-layer perceptron (MLP) that maps a
given query location $h_{n}$, within the prediction horizon $[0,h]$, to a
collection of trainable Trunk basis functions:

$\varphi:=(\varphi_{1}(h_{n}),\varphi_{2}(h_{n}),\ldots,\varphi_{q}(h_{n}))^{\top}\in\mathbb{R}^{q}.$

Finally, we compute the DeepGraphONet’s output for each node $i\in S$ by
merging the corresponding Branch coefficients $\mathbf{b}^{(i)}$ with the
Trunk basis functions $\varphi$ using the dot product:

$\mathcal{F}_{\theta}^{(i)}(\\{\mathbf{x}_{S}(t-\tau_{i})\\}_{i=1}^{m};S)(h_{n})=\sum_{j=1}^{q}\mathbf{b}_{j}^{(i)}\cdot\varphi_{j}(h_{n}).$

Training the DeepGraphONet $\mathcal{F}_{\theta}$. To train the parameters
$\theta$ of the proposed DeepGraphONet, we minimize, using gradient-based
optimization schemes (_e.g.,_ Adam [42]), the loss function:

$\mathcal{L}(\theta;\mathcal{D})=\frac{1}{N}\sum_{k=1}^{N}\left\lvert\mathbf{x}_{S}^{k}(t+h_{n}^{k})-\mathcal{F}_{\theta}(\\{\mathbf{x}_{S}^{k}(t-\tau_{i})\\}_{i=1}^{m})(h_{n}^{k})\right\rvert_{1}.$

using the dataset $\mathcal{D}$ of $N:=|\mathcal{D}|$ triplets:

$\mathcal{D}:=\left\\{\\{\mathbf{x}_{S}^{k}(t-\tau_{i})\\}_{i=1}^{m},h_{n}^{k},\mathbf{x}_{S}^{k}(t+h_{n}^{k})\right\\}_{k=1}^{N}.$

generated using the true solution operator $\mathcal{F}$.

###  III-B Resolution-independent DeepGraphONet

Fig. 2: The input for the Branch Net in standard and resolution-independent
DeepGraphONet. a. shows the input for the Branch Net in the standard
DeepGraphONet, where the network takes the input function values at a fixed
set of sensors. b. presents, in a resolution-independent DeepGraphONet, the
Branch Net takes both input function values and their locations. The fixed set
of sensors requirement is eliminated.

The vanilla DeepONet, introduced in [10], effectively and accurately
approximates single-input single-output solution operators $\mathcal{F}$.
Thus, DeepONet can predict at arbitrary locations within the output function
domain $[0,h]$. However, the main limitation of the vanilla DeepONet is that
it requires a fixed set of sensors within the input function domain
$[t-t_{M},t]$. This requirement forces the network to take new data points
with the same resolution, which, in turn, prevents the flexible application of
the trained network, _e.g.,_ in scenarios where (i) sampling may be inexact or
(ii) recursive prediction is needed.

To address such a limitation, we propose a simple, intuitive solution that
allows using resolution-independent sensors, as shown in Fig. 2. More
specifically, we enable the GNN-based Branch network to take as inputs not
only the graph-state memory input
$\\{\mathbf{x}_{S}(t-\tau_{i})\\}_{i=1}^{m}$, but also the corresponding,
possibly time-variant, sensor locations. To this end, for each node $i\in S$,
we construct the input function by concatenating the node-state function
values $\\{\mathbf{x}_{S}(t-\tau_{i})\\}_{i=1}^{m}$ with its corresponding
sensor location $\\{t-\tau_{i}\\}_{i=1}^{m}$. These function-location value
pairs provide complete information about the input, enabling us to relax the
fixed sensors’ requirement. In section IV, we will demonstrate that the
resolution-independent DeepGraphONet does not cause performance degradation
when compared to the standard DeepGraphONet.

###  III-C Zero-Shot Learning Scheme

In this paper, we build the proposed DeepGraphONet framework to be versatile.
That is, once we have trained the DeepGraphONet $\mathcal{F}_{\theta^{*}}$
using the graph $S$, we expect to use it with high accuracy on a different
graph $S^{\prime}\neq S$ (with possibly different number of nodes) for the
same task.

Such a zero-shot learning capability of the proposed DeepGraphONet is achieved
immediately due to the message-passing nature of the GNN-based Branch network.
In particular, we can directly apply DeepGraphONet to $S^{\prime}$ because the
trainable weights of the GNN-based Branch depend on the node feature dimension
and not on the adjacency matrix. Thus, the proposed DeepGraphONet is
inductive, and we may use it beyond the graph $S$ used during training. In
Section IV, we will show that DeepGraphONet can effectively approximate the
solution operator of networked dynamical systems within sub-graph structures
$S^{\prime}$ not used during training; thus, achieving zero-shot learning.

##  IV Experiments

To evaluate the accuracy and effectiveness of the proposed DeepGraphONet, we
test it on two different tasks: (i) predicting the transient response of power
grids (Section IV-A) and (ii) forecasting traffic on highway networks (Section
IV-B). In the experiments, we use $t_{M}$, $m$, $h$, and $h_{n}$ to denote the
memory length, number of sensors within the memory, prediction horizon, and
query location within the prediction horizon, respectively. We treat $t_{M}$
and $h$ as hyperparameters in this paper and explore the model performance.

Neural networks and implementation protocols. We use the same neural network
architecture to build two DeepGraphONets: the (i) standard DeepGraphONet,
which has fixed sensors as in [10], and (ii) resolution-independent
DeepGraphONet, which we proposed in this paper (see Fig. 2) to alleviate these
fixed sensors constraint. To build the DeepGraphONets’ architecture, we
proceed as follows. The Branch Net uses twenty message passing graph
convolution layers, and the Trunk Net uses five feed-forward layers. The
output dimension is 100 neurons for both the Branch and Trunk Nets.

We prepared the data to align with the DeepGraphONet’s formulation, where we
set a memory partition and an output function domain bound $h$. With regard to
the resolution-independent DeepGraphONet, with the memory partition size $m$,
we randomly sampled $\lfloor\frac{m}{2}\rfloor$ sensors and stacked them with
their locations as additional feature channels for the Branch Net.

Algorithm 1 shows the procedure of using trained DeepGraphONet to predict the
complete future trajectories. With the trained network, we use the initial
memory and partition to produce predictions covering the entire prediction
horizon. We then repeat this process with the newly observed ground truth as
the memory for the next prediction horizon until reaching the end of the
trajectory.

For example, in the power grid problem, we have a trained DeepGraphONet with
memory length $t_{M}=200\leavevmode\nobreak\ ms$ and prediction horizon
$h=20\leavevmode\nobreak\ ms$. In testing time, we want to produce predictions
of the entire trajectory of size $700\leavevmode\nobreak\ ms$ after the
initial memory. We start with using the initial memory and the trained network
to make predictions for the next 20 $ms$ and save the predictions to a list.
We then use the new observations to produce predictions for another 20 $ms$
until reaching the end of the trajectory. Finally, we concatenate the saved
predictions as the predicted complete trajectory.

We implemented the DeepGraphONet in PyTorch and trained and tested it with an
Nvidia GPU. The source code of the model can be accessed on GitHub 111The
repository will be public upon the acceptance of the manuscript. We measured
and compared the performance of all trained models in terms of $L_{1}$
relative error (%), which can be interpreted as the percentage absolute
deviation from the true values.

Require: system graph $S$; trained DeepGraphONet $\mathcal{F}_{\theta^{*}}$;
hyperparameters ($t_{M}$, $m$, $h$, $h_{n}$); initial local memory
$\\{\mathbf{x}_{S}(t-\tau_{i})\\}_{i=1}^{m}$; complete trajectory length $L$;

initialize an empty list for saving the predictions $C$;

let $N=\lceil\frac{(L-t_{M})}{h}\rceil$;

for  _$n=0,\ldots,N-1$_ do

for  _each $h_{n}\in[0,h]$_ do

use the DeepGraphONet’s forward pass to compute the prediction
$\mathbf{\bar{x}}(t+h_{n})=\mathcal{F}_{\theta^{*}}(\\{\mathbf{x}_{S}(t-\tau_{i})\\}_{i=1}^{m},\\{\tau\\}_{i=1}^{m};S)(h_{n}).$

append $\mathbf{\bar{x}}(t+h_{n})$ to $C$;

end for

update the current time  $t\leftarrow t+h;$ new observations come in
$\\{\mathbf{x}_{S}(t-\tau_{i})\\}_{i=1}^{m};$

end for

concatenate $C$;

Return the concatenated predicted trajectories;

Algorithm 1 Complete Trajectory Prediction Scheme

###  IV-A Experiment 1: Power transient stability

Fig. 3: One-line diagram of the New York-New England 16-generator 68-bus power
grid. We train our DeepGraphONet using the arbitrary subgraph of buses
$S:=\\{15,16,17,19,21,24\\}\subset G$, highlighted with red dots.

Dataset. The power grid dataset consisted of 1830 independent trajectories on
a sub-graph with 6 nodes (red buses in Fig. 4). For each trajectory, there
were $701$ time-steps with a resolution of 1 $ms$. For more details about the
data generation process, one can refer to [14]. Among the available
trajectories, $60\%$, $20\%$, and $20\%$ of them were used for training,
validation, and testing, respectively. We then applied the trained model to
the entire region (marked in blue in Figure 4), achieving zero-shot learning.

Fig. 4: Zero-shot learning on a larger subgraph. We train the DeepGraphONet on
an arbitrary small subgraph (marked in red) and directly apply the trained
network to a much larger subgraph (marked in blue) to achieve high predicting
performance (0.686% $L_{1}$ relative error).

Testing Accuracy. We show the model testing results in Fig. 5, Table I and II.
Overall, our trained both standard and resolution-independent DeepGraphONets
accurately approximated the solution operator for the power grid problem. The
$L_{1}$ relative errors of all short- and medium-term forecasting settings
were within 1%. Moreover, the resolution-independent DeepGraphONet, taking the
input function representation with arbitrary resolution, did not show
performance degradation compared to the standard (Figure 6b).

Varying Memory Lengths. We explored the impact and sufficiency of local memory
lengths. We trained the DeepGraphONet with different memory lengths,
$t_{M}=50\leavevmode\nobreak\ ms,t_{M}=100\leavevmode\nobreak\
ms,t_{M}=200\leavevmode\nobreak\ ms$, and evaluated the networks on the
testing dataset. The prediction results suggest that with larger $t_{M}$, the
average model performance improves (Fig. 5(d)(e)(f)). We observed the same
trend in both original and resolution-independent DeepGraphONets, as shown in
Table II.

Varying Horizons To investigate the performance change of DeepGraphONet with
respect to the change of the output function domain, we varied the prediction
horizons $h$. In particular, we trained the networks with
$h=2\leavevmode\nobreak\ ms,h=10\leavevmode\nobreak\
ms,h=20\leavevmode\nobreak\ ms,h=200\leavevmode\nobreak\ ms$ with the same
memory length, $t_{M}=50\leavevmode\nobreak\ ms$. The results (Fig. 5abc, and
Table I) show for all prediction horizons, both the standard and resolution-
independent DeepGraphONet made accurate approximations to the solution
operators. Furthermore, the network performance increases as $h$ decreases.
For short- and medium-term forecasting, the networks led to $L_{1}$ relative
errors of less than 1%, whereas when $h=200\leavevmode\nobreak\ ms$, the
relative error was over 2%. However, with longer local memory
($t_{M}=200\leavevmode\nobreak\ ms$), the trained DeepGraphONet accurately
predicted the longer horizon ($h=500\leavevmode\nobreak\ ms$), shown in Fig.
6a.

Fig. 5: Comparison of the DeepGraphONet predicted transient trajectories (dashed red lines) with the actual trajectories (solid black lines) on the subgraph of buses $S:=\\{15,16,17,19,21,24\\}\subset G$. DeepGraphONet generated all predictions at the same time given the memory trajectory (dashed gray lines). We test multiple prediction horizons $h$ and memory resolutions $t_{M}$. In particular, a. $h=2\leavevmode\nobreak\ ms$ with $t_{M}=50\leavevmode\nobreak\ ms$; b. $h=10\leavevmode\nobreak\ ms$ with $t_{M}=50\leavevmode\nobreak\ ms$; c. $h=20\leavevmode\nobreak\ ms$ with $t_{M}=50\leavevmode\nobreak\ ms$; d. $h=200\leavevmode\nobreak\ ms$ with $t_{M}=50\leavevmode\nobreak\ ms$; e. $h=20\leavevmode\nobreak\ ms$ with $t_{M}=100\leavevmode\nobreak\ ms$; and f. $h=20\leavevmode\nobreak\ ms$ with $t_{M}=200\leavevmode\nobreak\ ms$. Table I: The mean and standard deviation (st.dev.) of the $L_{1}$-relative error (%) between the predicted and actual transient response trajectories for multiple time horizons $h$ and using the (i) standard DeepGraphONet and (ii) resolution-independent DeepGraphONet. $h\leavevmode\nobreak\ (ms)$ |  | $2$ | $10$ | $20$ | $200$
---|---|---|---|---|---
standard | mean | 0.21 | 0.64 | 0.71 | 2.57
st.dev. | 0.74 | 2.71 | 2.62 | 4.89
res.ind. | mean | 0.37 | 0.62 | 0.72 | 1.71
st.dev. | 1.62 | 2.57 | 2.61 | 3.58
Table II: The mean and standard deviation (st.dev.) of the $L_{1}$-relative error (%) between the predicted and actual transient response trajectories for multiple memory resolutions $m$ and using the (i) standard DeepGraphONet and (ii) resolution-independent DeepGraphONet. $t_{M}\leavevmode\nobreak\ (ms)$ |  | $50$ | $100$ | $200$
---|---|---|---|---
standard | mean | 0.71 | 0.70 | 0.66
st.dev. | 2.62 | 2.82 | 3.14
res.ind. | mean | 0.72 | 0.71 | 0.68
st.dev. | 2.60 | 3.46 | 2.98
Fig. 6: Comparison of the DeepGraphONet predicted transient trajectories
(dashed red lines) with the actual trajectories (solid black lines) on the
subgraph of buses $S:=\\{15,16,17,19,21,24\\}\subset G$ using the proposed
resolution-independent DeepGraphONet. We test two time horizons $h$ and memory
resolutions $t_{M}$. a. $h=500\leavevmode\nobreak\ ms$ with
$t_{M}=200\leavevmode\nobreak\ ms$ and b. $h=2\leavevmode\nobreak\ ms$ with
$t_{M}=50\leavevmode\nobreak\ ms$. Fig. 7: Zero-shot learning results.
Comparison of the DeepGraphONet predicted transient trajectories (dashed red
lines) with the actual trajectories (solid black lines) on a sub-graph of
buses $S^{\prime}:=\\{2,3,4,5,6,7,18,25,27\\}\subset G$ different from the
training subgraph $S$.

Zero-shot Learning Finally, we applied the trained model to a larger sub-graph
from the same power grid model. The new sub-graph has 34 nodes and contains
the graph used during training. By directly applying the trained model to the
larger sub-graph without further training, we achieved zero-shot learning with
high accuracy. Fig. 7 depicts the prediction results on 9 randomly selected
nodes of the new sub-graph. For these nodes, the predictions aligned with the
ground truth well with a $L_{1}$ relative error of 0.686%.

With the above, we conclude the proposed DeepGraphONet can successfully learn
the solution operator with local memories for predicting the transient
stability of power grids. Particularly, the resolution-independent
DeepGraphONet, although having the flexibility of taking the discrete input
function representation with arbitrary resolution, does not show performance
degradation. In addition, one can directly apply the trained DeepGraphONet to
solve the power grid transient stability problem with a different underlying
graph structure, achieving zero-shot learning.

###  IV-B Experiment 2: Traffic forecasting

Fig. 8: METR-LA traffic network (adopted from [36]). The black dots represent
the loop detectors, which we regard as the nodes of the graph $G$.

Dataset. We used the METR-LA dataset [39], containing 179 loop detectors,
shown in Fig. 8. The dataset contained speed measurements of the year 2018,
with a resolution of 5 min. Due to the lack of independent trajectories, we
split the data based on the time for training, validation, and testing
purposes. We used the data from January to June 2018 for training, July to
September 2018 for validation and October to December 2018 for testing. We
trained and evaluated the network using the first 20 nodes in the graph and
directly applied the trained network to the entire graph for zero-shot
learning.

Varying Horizons. We investigated the DeepGraphONet’s performance in traffic
forecasting under the commonly used horizons [43], $h=15\leavevmode\nobreak\
min,h=30\leavevmode\nobreak\ min,h=60\leavevmode\nobreak\ min$ with the local
memory length fixed as $t_{M}=60\leavevmode\nobreak\ min$. Fig. 9 shows the
prediction results from the randomly picked 4 nodes in the graph where the
predicting horizon was $h=60\leavevmode\nobreak\ min$. DeepGraphONet captured
the future speed of all nodes at the same time accurately using the history
observation, while at some sharp turning points, the model did not match the
oscillation exactly but captured the average behavior. The results suggest
that as increasing horizon, the model performance decreased overall and the
error range became larger. Table III lists the mean and standard deviation of
the $L_{1}$ relative error of the models, both for original and resolution-
independent DeepGraphONet, with different predicting horizons on the testing
set. The results imply the great ability of the proposed DeepGraphONet in
predicting the traffic dynamics for all nodes present in the network
simultaneously, with $\sim 6$% error for short-term (5 min) prediction and
$\sim 11$% for long-term (1 hour).

Table III: The mean and standard deviation (st.dev.) of the $L_{1}$-relative error (%) between the predicted and actual traffic trajectories for multiple time horizons $T$ and using the (i) standard DeepGraphONet and (ii) resolution-independent DeepGraphONet. $h\leavevmode\nobreak\ (min)$ |  | $15$ | $30$ | $60$
---|---|---|---|---
standard | mean | 5.52 | 7.60 | 10.79
st.dev. | 9.49 | 14.74 | 22.05
res.ind. | mean | 6.84 | 8.62 | 11.87
st.dev. | 12.56 | 16.79 | 23.57
Fig. 9: Comparison of the DeepGraphONet predicted trajectories (dashed red
lines) with the actual trajectories (solid black lines) on the subgraph
$S:=\\{1,6,8,13\\}$ of the highway network $G$ using $t_{M}=60$ (minutes)
memory length. Fig. 10: Zero-shot learning results. Comparison of the
DeepGraphONet predicted trajectories (dashed red lines) with the actual
trajectories (solid black lines) on all the 179 nodes of the METR-LA highway
network $G$. We illustrate the predicted trajectories for the following four
nodes $\\{21,121,165,172\\}$, which we selected uniformly at random.

#### Zero-shot Learning

To test the DeepGraphONet’s potential of zero-shot learning for traffic
forecasting, we directly applied the trained network to the entire METR-LA
network with 179 nodes. Fig. 10 shows the visualized model prediction results
on 4 randomly picked nodes among the 179. With 1 hour predicting horizon, the
trained model closely predicted the future on the unseen graph structure for
all nodes at the same time, resulting in the $L_{1}$ relative error of 10.79%.

The experiment results for the traffic forecasting problem strongly suggest
our proposed DeepGraphONet was powerful in learning the solution operator for
the traffic dynamics. Also, the DeepGraphONet can be better implemented for
practical applications for it is easy to train, can produce predictions at
arbitrary time-step within the horizon, and does not require forward
dependency to obtain the value at this step.

##  V Discussion and Future Work

We used the proposed DeepGraphONet to learn the solution operator of networked
dynamical systems and demonstrated its effectiveness in power grid transient
stability and traffic forecasting tasks. In particular, the resolution-
independent DeepGraphONet demonstrated excellent accuracy in approximating the
solution operator without performance degradation while having the ability of
flexible discrete input function representation.

We demonstrated that the trained DeepGraphONet can produce predictions at
arbitrary query locations within the output function domain and does not
require forward dependency in the output sequence. In addition, we showed that
one could directly apply the DeepGraphONet, trained on a small sub-graph
region, to a larger graph region for the same task, achieving zero-shot
learning with high accuracy. With the promising wide-ranged applications of
the proposed DeepGraphONet, we provide next a preamble for our future work:

  * •

On Physics-Informed DeepGraphONet. Our current work focuses on developing a
physics-informed version of DeepGraphONet. We plan to use first-principles
models to (i) train the DeepGraphONet without data or (ii) use them as prior
information during data-driven training.

  * •

On Anomaly Detection with DeepGraphONet. We will focus part of our future work
on using DeepGraphoNet for detecting anomalies. To this end, we must endow the
proposed DeepGraphONet with the ability to estimate the network’s edges status
using graph trajectory data.

  * •

On Learning Networked Control Systems. We plan to extend the proposed
DeepGraphONet for learning networked control dynamical systems. This problem
is more challenging because we need to learn the system’s response using graph
trajectory and its response to external inputs. If successful, our
DeepGraphONet for control may become an essential tool for implementing model-
based multi-agent continuous reinforcement learning.

  * •

On Predicting Long-term Trajectories. We plan to develop effective training
strategies to alleviate error accumulation such that we can recursively use
the trained DeepGraphONet to forecast trajectories for long-term horizons with
high accuracy.

##  VI Conclusion

In this paper, we introduced the Deep Graph Operator Network (DeepGraphONet)
framework to predict the dynamics of complex systems with underlying sub-graph
structures. By fusing the ability from (1) Graph Neural Networks to correlate
graph trajectory information and (2) DeepONet to approximate nonlinear
operators, we achieved significant results on complex problems such as
predicting the (i) transient stability of a control area of a power grid and
(ii) the traffic flow within a city or district. Moreover, we built our
DeepGraphONet to be resolution-independent, _i.e.,_ we do not require a set of
fixed sensors to encode the graph trajectory history. Finally, we designed a
zero-shot learning strategy that enables using the proposed DeepGraphONet on a
different sub-graph with high predicting performance.

## Acknowledgment

This work was supported by the Advanced Grid Modeling Program, Office of
Electricity Delivery and Energy Reliability of the U.S. Department of Energy.
GL, CM, and YS gratefully acknowledge the support of the National Science
Foundation (DMS-1555072, DMS-2053746, and DMS-2134209), Brookhaven National
Laboratory Subcontract 382247, and U.S. Department of Energy (DOE) Office of
Science Advanced Scientific Computing Research program DE-SC0021142).

## References

  * [1] A. Iserles, A first course in the numerical analysis of differential equations.  No. 44, Cambridge university press, 2009.
  * [2] J. Li, M. Yue, Y. Zhao, and G. Lin, “Machine-learning-based online transient analysis via iterative computation of generator dynamics,” in 2020 IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids (SmartGridComm), pp. 1–6, IEEE, 2020.
  * [3] H. Wilms, M. Cupelli, A. Monti, and T. Gross, “Exploiting spatio-temporal dependencies for rnn-based wind power forecasts,” in 2019 IEEE PES GTD Grand International Conference and Exposition Asia (GTD Asia), pp. 921–926, IEEE, 2019.
  * [4] Y. Sun, T. Mallick, P. Balaprakash, and J. Macfarlane, “A data-centric weak supervised learning for highway traffic incident detection,” arXiv preprint arXiv:2112.09792, 2021.
  * [5] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” Proceedings of the national academy of sciences, vol. 113, no. 15, pp. 3932–3937, 2016.
  * [6] E. Kaiser, J. N. Kutz, and S. L. Brunton, “Sparse identification of nonlinear dynamics for model predictive control in the low-data limit,” Proceedings of the Royal Society A, vol. 474, no. 2219, p. 20180335, 2018.
  * [7] H. Schaeffer, “Learning partial differential equations via data discovery and sparse optimization,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 473, no. 2197, p. 20160446, 2017.
  * [8] T. Qin, K. Wu, and D. Xiu, “Data driven governing equations approximation using deep neural networks,” Journal of Computational Physics, vol. 395, pp. 620–635, 2019.
  * [9] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Multistep neural networks for data-driven discovery of nonlinear dynamical systems,” arXiv preprint arXiv:1801.01236, 2018.
  * [10] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis, “Learning nonlinear operators via deeponet based on the universal approximation theorem of operators,” Nature Machine Intelligence, vol. 3, no. 3, pp. 218–229, 2021\.
  * [11] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, “Neural operator: Graph kernel network for partial differential equations,” arXiv preprint arXiv:2003.03485, 2020.
  * [12] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, “Fourier neural operator for parametric partial differential equations,” arXiv preprint arXiv:2010.08895, 2020.
  * [13] Tianping Chen and Hong Chen, “Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems,” IEEE Transactions on Neural Networks, vol. 6, pp. 911–917, July 1995.
  * [14] C. Moya, S. Zhang, M. Yue, and G. Lin, “Deeponet-grid-uq: A trustworthy deep operator framework for predicting the power grid’s post-fault trajectories,” arXiv preprint arXiv:2202.07176, 2022.
  * [15] S. Zhang and G. Lin, “Robust data-driven discovery of governing physical laws with error bars,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 474, no. 2217, p. 20180305, 2018.
  * [16] I. Goodfellow, Y. Bengio, and A. Courville, Deep learning.  MIT press, 2016.
  * [17] N. Winovich, K. Ramani, and G. Lin, “Convpde-uq: Convolutional neural networks with quantified uncertainty for heterogeneous elliptic partial differential equations on varied domains,” Journal of Computational Physics, vol. 394, pp. 263–279, 2019.
  * [18] S. Cai, Z. Wang, L. Lu, T. A. Zaki, and G. E. Karniadakis, “Deepm&mnet: Inferring the electroconvection multiphysics fields based on operator approximation by neural networks,” Journal of Computational Physics, vol. 436, p. 110296, 2021.
  * [19] S. Goswami, M. Yin, Y. Yu, and G. Karniadakis, “A physics-informed variational deeponet for predicting the crack path in brittle materials,” arXiv preprint arXiv:2108.06905, 2021.
  * [20] P. Jin, S. Meng, and L. Lu, “Mionet: Learning multiple-input operators via tensor product,” arXiv preprint arXiv:2202.06137, 2022.
  * [21] Z. Wu, S. Pan, F. Chen, G. Long, C. Zhang, and S. Y. Philip, “A comprehensive survey on graph neural networks,” IEEE transactions on neural networks and learning systems, vol. 32, no. 1, pp. 4–24, 2020.
  * [22] A. Sperduti and A. Starita, “Supervised neural networks for the classification of structures,” IEEE Transactions on Neural Networks, vol. 8, no. 3, pp. 714–735, 1997.
  * [23] M. Gori, G. Monfardini, and F. Scarselli, “A new model for earning in raph domains,” vol. 2, pp. 729 – 734 vol. 2, 01 2005.
  * [24] J. B. Estrach, W. Zaremba, A. Szlam, and Y. LeCun, “Spectral networks and deep locally connected networks on graphs,” in 2nd international conference on learning representations, ICLR, vol. 2014, 2014.
  * [25] D. K. Hammond, P. Vandergheynst, and R. Gribonval, “Wavelets on graphs via spectral graph theory,” Applied and Computational Harmonic Analysis, vol. 30, no. 2, pp. 129–150, 2011.
  * [26] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” arXiv preprint arXiv:1609.02907, 2016.
  * [27] Z. Wu, S. Pan, F. Chen, G. Long, C. Zhang, and P. S. Yu, “A Comprehensive Survey on Graph Neural Networks,” IEEE Transactions on Neural Networks and Learning Systems, vol. 32, no. 1, pp. 4–24, 2021.
  * [28] J. Zhou, G. Cui, S. Hu, Z. Zhang, C. Yang, Z. Liu, L. Wang, C. Li, and M. Sun, “Graph neural networks: A review of methods and applications,” AI Open, vol. 1, no. April, pp. 57–81, 2020.
  * [29] J. Atwood and D. Towsley, “Diffusion-convolutional neural networks,” Advances in neural information processing systems, vol. 29, 2016.
  * [30] M. Niepert, M. Ahmed, and K. Kutzkov, “Learning convolutional neural networks for graphs,” in International conference on machine learning, pp. 2014–2023, PMLR, 2016.
  * [31] J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl, “Neural message passing for quantum chemistry,” in International conference on machine learning, pp. 1263–1272, PMLR, 2017.
  * [32] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Lio, and Y. Bengio, “Graph attention networks,” arXiv preprint arXiv:1710.10903, 2017.
  * [33] W. Hamilton, Z. Ying, and J. Leskovec, “Inductive representation learning on large graphs,” Advances in neural information processing systems, vol. 30, 2017.
  * [34] S. Guo, Y. Lin, N. Feng, C. Song, and H. Wan, “Attention based spatial-temporal graph convolutional networks for traffic flow forecasting,” in Proceedings of the AAAI conference on artificial intelligence, vol. 33, pp. 922–929, 2019.
  * [35] Y. Li, R. Zemel, M. Brockschmidt, and D. Tarlow, “Gated graph sequence neural networks,” 4th International Conference on Learning Representations, ICLR 2016 - Conference Track Proceedings, no. 1, pp. 1–20, 2016.
  * [36] Y. Li, R. Yu, C. Shahabi, and Y. Liu, “Diffusion convolutional recurrent neural network: Data-driven traffic forecasting,” arXiv preprint arXiv:1707.01926, 2017.
  * [37] L. Zhao, Y. Song, C. Zhang, Y. Liu, P. Wang, T. Lin, M. Deng, and H. Li, “T-gcn: A temporal graph convolutional network for traffic prediction,” IEEE Transactions on Intelligent Transportation Systems, vol. 21, no. 9, pp. 3848–3858, 2019.
  * [38] J. Zhu, Q. Wang, C. Tao, H. Deng, L. Zhao, and H. Li, “Ast-gcn: Attribute-augmented spatiotemporal graph convolutional network for traffic forecasting,” IEEE Access, vol. 9, pp. 35973–35983, 2021.
  * [39] H. V. Jagadish, J. Gehrke, A. Labrinidis, Y. Papakonstantinou, J. M. Patel, R. Ramakrishnan, and C. Shahabi, “Big data and its technical challenges,” Communications of the ACM, vol. 57, no. 7, pp. 86–94, 2014.
  * [40] A. J. Chorin and O. H. Hald, Stochastic tools in mathematics and science, vol. 1.  Springer, 2009.
  * [41] A. J. Chorin, O. H. Hald, and R. Kupferman, “Optimal prediction with memory,” Physica D: Nonlinear Phenomena, vol. 166, no. 3-4, pp. 239–257, 2002.
  * [42] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
  * [43] Z. Wu, S. Pan, G. Long, J. Jiang, X. Chang, and C. Zhang, “Connecting the dots: Multivariate time series forecasting with graph neural networks,” in Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 753–763, 2020.