# Topologically protected dynamics in three-dimensional nonlinear
antisymmetric Lotka-Volterra systems

Muhammad Umer  0000-0002-1941-1833<EMAIL_ADDRESS>Department of Physics,
National University of Singapore, Singapore 117551, Republic of Singapore.
Jiangbin Gong  0000-0003-1280-6493<EMAIL_ADDRESS>Department of Physics,
National University of Singapore, Singapore 117551, Republic of Singapore.

###### Abstract

Studies of topological bands and their associated low-dimensional boundary
modes have largely focused on linear systems. This work reports robust
dynamical features of three-dimensional ($3$D) nonlinear systems in connection
with intriguing topological bands in $3$D. Specifically, for a $3$D setting of
coupled rock-paper-scissors cycles governed by the antisymmetric Lotka-
Volterra equation (ALVE) that is inherently nonlinear, we unveil distinct
characteristics and robustness of surface polarized masses and analyze them in
connection with the dynamics and topological bands of the linearized Lotka-
Volterra (LV) equation. Our analysis indicated that insights learned from Weyl
semimetal phases with type-I and type-II Weyl singularities based on a
linearized version of the ALVE are still remarkably useful, even though the
system dynamics is far beyond the linear regime. This work indicates the
relevance and importance of the concept of topological boundary modes in
analyzing high-dimensional nonlinear systems and hopes to stimulate another
wave of topological studies in nonlinear systems.

Introduction.— The concept of topological robustness was introduced to
condensed-matter physics in 1980’s [1, 2, 3]. Continued studies on topological
matter have led to many important discoveries of new phases of matter [4,
*Bernevig2006, *Fu2007, 7, *Qi2011, 9, *Wan2011, *Zyuzin2012, *Hosur2013,
*Huang2015, 14, *Bomantara2016a, 16, *Schindler2018, *Li2018, 19], including
topological insulators [4, *Bernevig2006, *Fu2007] and superconductors [7,
*Qi2011, 20], Weyl semimetals [9, *Wan2011, *Zyuzin2012, *Hosur2013,
*Huang2015, 14, *Bomantara2016a], as well as semimetals with linked and/or
knotted nodal lines [21]. Beyond the usual condensed-matter context,
topological phases of matter, especially their robust boundary states
protected by topology, have advanced studies of photonic crystals [22,
*Raghu2008, *Wang2009], acoustic systems [25], mechanical meta-materials [26,
*Suesstrunk2015, *Nash2015], soft matter [29, *Zhou2018, *Pedro2019], and
biological [32, *Yamauchi2020] systems. Remarkably, though topological band
theory is based on the linear Schrödinger (or Schrödinger-like) equation with
translational invariance, it is also relevant to the understanding of robust
boundary behavior in inherently nonlinear systems [34, 35, 36]. This
recognition is expanding the territory of topological physics with a
potentially long-term impact.

In this letter, we report how a well-established concept of topological
gapless phases in $3$D systems, known as Weyl semimetal phases and their
variants, can emerge in, and guide our understanding of a class of nonlinear
systems, the so-called antisymmetric Lotka-Volterra equation (ALVE) system of
coupled rock-paper-scissors (RPS) cycles in a $3$D configuration.

Specifically, polarized masses in $3$D ALVE systems are found to display
distinct behaviors in different parameter regimes, with the time-averaged
masses localized at the surface. Further, the propagation of masses is chiral
at certain surfaces and surface polarized time-averaged masses are robust
against the perturbation of system parameters. Such robust features hint a
possible important role of topology concerning $3$D topological phases. To
reveal the underlying topological physics, we linearize the LV equation, which
then resembles to the Schrödinger equation of spinless particle on a $3$D
lattice. Indeed, the resultant linearized equation of motion is analogous to
that in symmetry class-$A$ of the ten-fold way topological classification in
condensed-matter physics. Surprisingly, the qualitative dynamical behavior of
the nonlinear system is in parallel with that of the linearized LV equation,
such as boundary propagation of probability density and localization of time-
averaged probability density at certain surfaces of the $3$D network. A
topological band analysis hence becomes necessary and useful: there are rich
gapless Weyl semimetal phases, including type-I Weyl semimetal, type-II Weyl
semimetal, and hybrid Weyl semimetal with coexisting type-I and type-II Weyl
nodes. As seen below, these various Weyl semimetal phases are critical to
digest different dynamical behaviors in the nonlinear system. These findings
hence establish a strong link between topological band theory and $3$D
nonlinear systems. This work also indicates that a number of distinct
topological semimetal phases may be effectively realized in nonlinear systems,
even though the dynamics is far beyond the linear regime.

Figure 1: Intra-cycle $t_{a}$ (red) and inter-cycle $t_{b},~{}t_{c}$ (green,
blue) couplings in $x-y$ plane and $z$ direction are shown in panel (a, b).
Magenta coloured path $\tilde{s}$ and arrows in panel (a) indicate the
direction of chiral propagating masses an each layer while (c) shows that each
point along the path $\tilde{s}$ has $N_{z}$ number of nodes stacked in the
$z$ direction. Moreover, magenta path $\tilde{s}$ shows hard boundary for ALVE
system while partially transparent nodes in (a, b) show system under periodic
boundary conditions (PBC) in $x$, $y$, and $z$ direction.

ALVE model system.— ALVE is a nonlinear model that describes the coexistence
of species in game theory [37, *Knebel2015, *Toupo2015, *Geiger2018] and
population dynamics [41]. The ALVE is given as

$\displaystyle\begin{aligned}
\frac{\partial{u_{l}}}{\partial{t}}=u_{l}\sum^{M}_{m=1}\alpha_{lm}u_{m}\;,\end{aligned}$
(1)

where $\alpha_{lm}$'s are elements of the antisymmetric payoff matrix ${\bf
A}$ with $\alpha_{ml}=-\alpha_{lm}$, $u_{l}~{}(u_{m})$ is the mass of the
specie at site $l~{}(m)$, and $M$ is the total number of sites in the system.
Mass on each site $u_{l}$ evolve under nonlinear interaction between masses
$\sim u_{l}u_{m}$, where total mass remains conserved
$\sum_{l=1}^{M}\partial_{t}u_{l}=0$.

To explore the relevance of $3$D topological band theory in higher-dimensional
nonlinear systems, we construct now a $3$D network system composed of stacked
Kagome layers of coupled RPS game cycles. Each layer has $N$ number of sites
arranged on the Kagome geometry in $x-y$ plane, as shown in Fig. 1a. We then
stack $N_{z}$ such layers in the $z$ direction for a $3$D network, as shown in
Fig. 1b. Of more importance is the rules of the $3$D game. Within one or
across different nearest neighbouring RPS cycles, the proposed rules are given
as; R dominates S, S dominates P, and P dominates R, such that each strategy
dominates (dominated by) only one other strategy [42, *Claussen2008,
*Szolnoki2014]. These rules are implemented through intra-cycle $t_{a}$,
inter-cycle payoff elements $t_{b}$ and $t_{c}$ in the $x-y$ plane and $z$
direction, respectively. The resulting ${\bf A}=\\{\alpha_{lm}\\}$ matrix is
then found to be

$\displaystyle\begin{aligned} &{\bf
A}=\sum_{x=1}^{N_{x}}\sum_{y=1}^{N_{y}}\sum_{z=1}^{N_{z}}\Big{(}t_{a}\big{[}\text{S}^{\dagger}_{x,y,z}\text{P}_{x,y,z}+\text{R}^{\dagger}_{x,y,z}\text{S}_{x,y,z}+\text{P}^{\dagger}_{x,y,z}\text{R}_{x,y,z}\big{]}\\\
&+t_{b}\big{[}\text{P}^{\dagger}_{x,y,z}\text{R}_{x+1,y,z}-\text{S}^{\dagger}_{x,y,z}\text{R}_{x,y+1,z}+\text{S}^{\dagger}_{x,y,z}\text{P}_{x-1,y+1,z}\big{]}\\\
&+t_{c}\big{[}\text{P}^{\dagger}_{x,y,z}\text{R}_{x,y,z+1}+\text{P}^{\dagger}_{x,y,z+1}\text{R}_{x,y,z}+\text{R}^{\dagger}_{x,y,z}\text{S}_{x,y,z+1}\\\
&+\text{R}^{\dagger}_{x,y,z+1}\text{S}_{x,y,z}+\text{S}^{\dagger}_{x,y,z}\text{P}_{x,y,z+1}+\text{S}^{\dagger}_{x,y,z+1}\text{P}_{x,y,z}\big{]}-h.c\Big{)},\end{aligned}$
(2)

where $\text{C}^{\dagger}_{x,y,z}$ ($\text{C}_{x,y,z}$) are creation
(annihilation) bosonic operator for specie
$\text{C}\in\\{\text{R},\text{P},\text{S}\\}$ at site $(x,~{}y,~{}z)$ of the
$3$D network (See Supplemental Materials 111see Supplemental Material at
insert link for details for matrix form) and $h.c$ is the hermitian conjugate.

Here we assume that nodes are numbered first within a single layer, and then
the same counting order continues to the second, third layer, etc.
Furthermore, we consider strictly positive initial mass
$u_{l}(t=0)=\sum_{j=1}^{N_{z}}\delta_{l,jN}+1/M$ on each site $l$ such that
mass is localized at $N^{\rm th}$-site (also at the hinge) of each layer in
the system and rest of the nodes have background mass $\approx~{}1/M$. Then,
${\bf u}(t=0)$ can be regarded as a perturbed state away from the strictly
positive steady-state (Nash equilibrium state) ${\bf c}$ that satisfies ${\bf
Ac}=0$, with $c_{l}=1/M$. The dynamical features we explore below are in the
vicinity of the system’s steady-state.

Figure 2: Nonlinear evolution of masses on the surface $\tilde{s}$ of $3$D
network are presented for parameter values (a) $t_{a}=1.0$, $t_{c}=3/8$, (b)
$t_{a}=1.4$, $t_{c}=3/8$, (c) $t_{a}=1.0$, $t_{c}=5/8$, and (d) $t_{a}=3.0$,
$t_{c}=4/10$, for fixed $t_{b}=1.0$.  Figure 3: Time-averaged profile of
masses is shown for parameter values (a) $t_{a}=1.0$, $t_{c}=5/8$,
$\Delta=0.2$ and (b) $t_{a}=3.0$, $t_{c}=4/10$, for fixed $t_{b}=1.0$, $N=204$
and $N_{z}=5$. Here $z$ indicate the layer in $z$-direction and we observe
that for a system with $N_{z}=5$, $z=1~{}(2)$ and $z=5~{}(4)$ have similarly
localization of masses.

Dynamical features in $3$D ALVE system.— The ALVE [Eq. (1)] is computationally
simulated and analysed, with $N=52$ and $N_{z}=5$. First, we consider
parameter values $t_{a}=1.0$ and $t_{c}=3/8$ (justification for this choice to
be elaborated later). To better visualize the chiral time propagation of
masses, we first define a path $\tilde{s}$ (see magenta coloured path in Fig.
1a) that only contains the boundary nodes and passes through $1^{\rm st}$ and
$N^{\rm th}$ nodes of each layer. At each position along the path $\tilde{s}$,
there are $N_{z}$ different nodes due to their different coordinates along the
$z$ direction. Fig. 2a depicts the time-evolving masses along the identified
path $\tilde{s}$ (but also recording different behaviors at different layers),
as a function of evolution time. Clearly, one sees that the movement of the
masses is confined to the boundary of the system and can make a round trip
(with an average velocity of $\nu=0.0430$) if projected onto the $x-y$ plane.
Another parameter selection, namely, $t_{a}=1.4,~{}t_{c}=3/8$, correspond to
analogous behavior of instantaneous masses $u_{l}(t)$ along the coordinate
$\tilde{s}$ as shown in Fig. 2b, where average velocity ($\nu=0.0470$) is
relatively higher such that one complete round trip across the path
$\tilde{s}$ takes less the time needed in Fig. 2a.

We now consider parameter values $t_{a}=1.0,~{}t_{c}=5/8$. In this case, as
shown in Fig. 2c, chiral propagation of masses splits into multiple waves
moving at different velocities, a phenomenon distinctively different from the
previously two situations. Specifically, the characteristic velocity of mass
propagation of the middle wave in Fig. 2c appears to assume two different
values during two time windows. The upper wave 1 and middle wave in Fig. 2c as
bounded by the two green lines, propagate at a characteristic velocity of
$\nu_{1}=0.020$ node per unit time, whereas the lower wave 2 and middle wave
as bounded by red lines, propagate at a characteristic velocity
$\nu_{2}=0.010$ node per unit time such that $\nu_{1}/\nu_{2}\approx 2$.
Overall, it is observed that mass chiral propagation now has two
characteristic velocities. To appreciate boundary localization, we consider
the time-averaged behavior as some quasi-steady feature. In particular, we
present the profile of time-averaged masses $\langle{u_{l}}\rangle_{T}$ in
Fig. 3a, where boundary localization of time-averaged masses in $x-y$ plane
can be observed, a feature shared by different layers.

It should also be highlighted that the above-observed localization behavior of
the mass propagation is robust against the disorder among the system
parameters. For example, disordered payoff matrix elements can be modelled by
$\alpha_{lm}\pm\delta_{l}$ and $\alpha_{ml}\mp\delta_{l}$, where
$\delta_{l}\in(-\Delta,\Delta)$ is derived from a uniform distribution of real
numbers. In Fig. 3a where time-averaged surface polarized masses are
presented, we have already averaged over one hundred realizations of disorder
among the payoff matrix elements.

To conclude our investigations of different parameter regimes, we now look
into the dynamics for parameter values $t_{a}=3.0,~{}t_{c}=4/10$. Fig. 2d
indicates a totally different feature: the mass propagation is not chiral
(without a clear velocity to identify) and the mass distribution is
delocalized. Moreover, the time-averaged masses are shown in Fig. 3b for
different layers. Delocalization of time-averaged masses across the entire
$3$D system in Fig. 3b suggest the lack of a mechanism to confine the mass
propagation at the surface of the system.

Above analysis implies that the $3$D ALVE system in different parameter
regimes may be connected with distinct topological phases. Indeed, the
observed robust boundary localization behavior in Fig. 2a-2c, and Fig. 3a
resembles to the surface states protected by certain bulk topology in linear
lattice models. If that is the case, then the absence of chiral mass
propagation and boundary localization in Fig. 2d and Fig. 3b suggests a
topologically trivial regime. To establish such connections, we linearize the
ALVE system and carry out topological band theory analysis of the linearized
system.

Linearized LV equation and topological band analysis.— We reexpress the mass
at node $l$ as $u_{l}=c_{l}+\delta{u_{l}}$, where $c_{l}$ is the steady-state
mass and treat $\delta{u_{l}}$ as a small quantity. A linear approximation
emerges once we take $\delta{u_{l}}\delta{u_{m}}\approx 0$, such that we have
the following linearized LV equation [35],

$\displaystyle\begin{aligned} i\partial_{t}{\bf\delta{u}}=\frac{i}{M}{\bf
A}{\bf\delta{u}}\;.\end{aligned}$ (3)

It is noted that linearized LV equation becomes the Schrödinger equation for
${\bf H}=i{\bf A}$ (up to a multiplicative factor of $1/M$, which only
rescales the energy eigenvalues). Here the Hermitian matrix ${\bf H}$ can be
taken as the Hamiltonian of a tight-binding system of spinless particle on a
$3$D lattice.

Before we dive into the topological band theory analysis, let us examine the
dynamical features of the linearized system described by Eq. (3).
Interestingly, we observe similar qualitative features (see Supplemental
Materials [45] for numerical results) as in the nonlinear system described by
Eq. (1), namely, surface propagation of probability density and surface
localization of time-averaged probability density. The probability density
wave propagate with some definite average velocity for parameter values
similar to that of Fig. 2a and Fig. 2b. Moreover, probability density wave
split into multiple waves propagating with two distinct velocities for the
parameter values as given for Fig. 2c and these velocities have a ratio of
$\approx 2$ [45]. It is worth mentioning that the absolute velocities of
probability density waves in the linearized systems here are different
(quantitatively) as compared to that of surface propagating masses in
nonlinear system. Such difference can be understood via a self-trapping
mechanism typically induced by nonlinear interactions [46] (note that the
self-trapping mechanism may boost the chiral wave propagation velocity).
Interestingly, if we initiate the nonlinear dynamics by considering the
starting configuration much closer to the steady-state, then the average
velocity of mass propagation does approach that of the linear system (See
Supplemental Materials [45]). Finally, further suggesting the necessity of
connecting the nonlinear dynamics with a linearized version, note that for
parameter values in Fig. 2d, the probability density of the linear system does
not propagate at the surface of the system and the time-averaged probability
densities are also delocalized across the $3$D network.

We now introduce a Fourier transformation, under PBC in all spatial
directions, from position to momentum space such that the bulk Hamiltonian
$H({\bf k})$ in the momentum representation becomes

$\displaystyle\begin{aligned} &H({\bf
k})=t_{b}\sin(k_{x})\gamma_{1}-t_{b}\sin(k_{y})\gamma_{4}+t_{b}\sin(-k_{x}+k_{y})\gamma_{6}\\\
&+[t_{a}+t_{b}\cos(k_{x})+2t_{c}\cos(k_{z})]\gamma_{2}-[t_{a}+t_{b}\cos(k_{y})\\\
&+2t_{c}\cos(k_{z})]\gamma_{5}+[t_{a}+t_{b}\cos(-k_{x}+k_{y})+2t_{c}\cos(k_{z})]\gamma_{7}\;,\end{aligned}$
(4)

where $\gamma_{j}$'s are Gell-Mann matrices [45]. The Hamiltonian Eq. (4) of
linearized LV equation belongs to class-$A$ of ten-fold way symmetry
classification of topological gapped systems [47, *Ryu2010] and does not
possess any internal symmetry such as charge conjugation, inversion, or time-
reversal symmetry. Because symmetry class-$A$ has trivial classification for
gapped system in $3$D, only non-trivial gapless phases can emerge here.

Figure 4: Parameter values are taken as (a, b) $t_{a}=1.0,t_{c}=3/8$, (c, d)
$t_{a}=1.4,t_{c}=3/8$, and (e, f) $t_{a}=1.0,t_{c}=5/8$, where (a, c, e) shows
the slice Chern number $C(k_{z})$, (b, d) shows linear (tilted) dispersion of
type-I (type-II) Weyl node for fixed $k_{y}=2\pi/3$, and (f) shows coexisting
type-I and type-II Weyl nodes for $(k_{x},k_{y})=(-2\pi/3,2\pi/3)$. These
type-I and type-II Weyl nodes have linear dispersion in $x-y$ direction such
that $|\nu_{x}^{I}|=|\nu_{y}^{I}|=0.0022$ and
$|\nu_{x}^{II}|=|\nu_{y}^{II}|=0.0044$ respectively, where we define
$\nu^{I/II}_{i}=\partial{E}/\partial_{k_{i}}$. Moreover, it is noted that
$\nu_{i}^{II}/\nu_{i}^{I}\approx 2$ for $i=x~{}(y)$.

For fixed $t_{b}=1.0$ without loss of generality, one can diagonalize the
system for energy eigenvalues [45] and it can be noted that two out of three
bands touch at distinct momenta points given as
$(k_{x_{0}},k_{y_{0}},k_{z_{0}})=(2\pi/3,-2\pi/3,\pm\cos^{-1}[\frac{\beta-
t_{a}}{2t_{c}}])$, $(-2\pi/3,2\pi/3,\pm\cos^{-1}[\frac{\beta-t_{a}}{2t_{c}}])$
where inequality $|\frac{\beta-t_{a}}{2t_{c}}|<1.0$ for
$\beta\in\\{t_{b}/2,2t_{b}\\}$ captures the band touching along $k_{z}$
momenta component. As we show in the following, these band touching points
carry non-trivial topological charge. We thus refer to these points as Weyl
nodes and corresponding phases as Weyl semimetals.

First, we consider parameter values $t_{a}=1.0$, $t_{c}=3/8$ similar to that
of Fig. 2a and inequality $|\frac{\beta-t_{a}}{2t_{c}}|<1.0$ holds for
$\beta=t_{b}/2$. The linearized LV system exhibit four type-I Weyl nodes,
where Fig. 4b shows one of them for fixed $k_{y}=2\pi/3$ with linear
dispersion in $k_{x}-k_{z}$ momenta plane. Two Weyl nodes appear at two
different energy values, for fixed $k_{z}=k_{z_{0}}$, where upper (lower) band
touches the middle band [45]. Charge and chirality of Weyl node located at
(upper) lower energy can then be inferred from the (negative) change of Chern
number [2] of (upper most) lowest band (-$C_{3}$) $C_{1}$ in Fig. 4a, while
the Chern number of middle band $C_{2}(k_{z})$ is zero and does not change due
to simultaneous touching with lower and upper bands at fixed $k_{z_{0}}$ (see
[45] for surface-state Fermi arc spectrum).

Second, we consider parameter values $t_{a}=1.4$, and $t_{c}=3/8$ similar to
that of Fig. 2b and inequality $|\frac{\beta-t_{a}}{2t_{c}}|<1.0$ holds for
$\beta=2t_{b}$. The system exhibits a total of four type-II Weyl nodes where
one of them is shown in Fig. 4d whose tilted nonlinear (flat) dispersion along
$k_{z}$ can be observed for fixed $k_{y}=2\pi/3$. As we discussed above, the
Chern number of lower $C_{1}$ and negative of upper band $-C_{3}$ captures the
charge and chirality of these Weyl nodes as we vary the momenta $k_{z}$ in
Fig. 4c.

Third, we consider parameter values $t_{a}=1.0,t_{c}=5/8$ similar to that of
Fig. 2c and inequality $|\frac{\beta-t_{a}}{2t_{c}}|<1.0$ now holds for
$\beta=t_{b}/2$ and $\beta=2t_{b}$ simultaneously. The linearized LV system
has a total of eight Weyl nodes [45] as compared to four in the previous two
cases, out of which half (other half) of them are type-I (type-II). Linear and
tilted nonlinear dispersion of these coexisting type-I and type-II Weyl nodes
can be observed in Fig. 4f for fixed $(k_{x},k_{y})=(2\pi/3,-2\pi/3)$, where
$\pm$ sign shows the chirality of corresponding Weyl node which can be duly
inferred from Chern number $C(k_{z})$ in Fig. 4e. Moreover, one can obtain
topological trivial phase for parameter values similar to that of Fig. 2d and
Fig. 3b, where linearized LV system neither has any Weyl node nor surface
states and each band has trivial topology [45].

Discussion.— We have shown that dynamical features of $3$D nonlinear ALVE
system exhibit topological features of polarization and robustness against the
disorder of system parameters. The underlying linearized LV equation also
exhibits similar qualitative features of surface propagation and localization.
Chiral propagation of mass (probability density) wave can be associated with
topological Weyl semimetal phases (of type-I or type-II) for the linearized LV
equation. On the other hand, chiral propagation of multiple mass waves in Fig.
2c can be associated with a hybrid Weyl semimetal phase, where type-I and
type-II Weyl nodes coexist in the linearized system. Non-chiral and
delocalized propagation in Fig. 2d, 3b are explained via a topological trivial
phase. These results clearly indicate the important relevance of topological
band analysis for $3$D nonlinear systems.

With RPS game cycles experimentally studied in social circles [49, 50], a $3$D
construction of such platforms is necessary to uncover the implications of
high-dimensional topological phases in nonlinear systems. Above discussed
linearized LV system may be also realized in electrical circuits [51] due to
the full control over circuit design. As encouraged by the richness of
semimetal phases found from $3$D nonlinear models here, it is also possible to
envision the relevance of higher-order topological phases in nonlinear
systems.

It is a pleasure to acknowledge helpful discussions with Raditya Weda
Bomantara. J. Gong acknowledges support from Singapore NRF Grant No. NRF-
NRFI2017-04 (WBS No. R-144-000-378-281).

Supplementary Materials

This supplementary material consists of three section. In Sec. I, we present
the antisymmetric ${\bf A}$ matrix and discuss the localization of time-
averaged masses in nonlinear system for various parameter values. In Sec. II,
we present the numerical results for the dynamics of linearized system.
Furthermore, we discuss the dynamics of nonlinear system under different
initial state with varying perturbations and show that the dynamics of
nonlinear system approaches that of linear system for smaller perturbations.
In Sec. III, we discuss the topological band analysis, band spectrum and
surface state Fermi arcs.

##  I Antisymmetric Matrix

First, we present a small block of antisymmetric payoff matrix A from Eq. (2)
of the main text, to better illustrate the payoff elements $t_{j}$ 's for
$j\in\\{a,b,c\\}$. This small block in the matrix form is given as,

$\displaystyle{\bf
A}=\setcounter{MaxMatrixCols}{15}\begin{pmatrix}0&-t_{a}&t_{a}&0&-t_{b}&0&0&0&t_{b}&0&\dots&0&-t_{c}&t_{c}&\dots\\\
t_{a}&0&-t_{a}&0&0&0&0&0&0&-t_{b}&\dots&t_{c}&0&-t_{c}&\dots\\\
-t_{a}&t_{a}&0&0&0&0&0&0&0&0&\dots&-t_{c}&t_{c}&0&\dots\\\
0&0&0&0&-t_{a}&t_{a}&0&0&0&0&\dots&0&0&0&\dots\\\
t_{b}&0&0&t_{a}&0&-t_{a}&0&0&0&0&\dots&0&0&0&\dots\\\
0&0&0&-t_{a}&t_{a}&0&0&0&0&0&\dots&0&0&0&\dots\\\
0&0&0&0&0&0&0&-t_{a}&t_{a}&0&\dots&0&0&0&\dots\\\
0&0&0&0&0&0&t_{a}&0&-t_{a}&0&\dots&0&0&0&\dots\\\
-t_{b}&0&0&0&0&0&-t_{a}&t_{a}&0&0&\dots&0&0&0&\dots\\\
0&t_{b}&0&0&0&0&0&0&0&0&\dots&0&0&0&\dots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots\\\
0&-t_{c}&t_{c}&0&0&0&0&0&0&0&\dots&0&-t_{a}&t_{a}&\dots\\\
t_{c}&0&-t_{c}&0&0&0&0&0&0&0&\dots&t_{a}&0&-t_{a}&\dots\\\
-t_{c}&t_{c}&0&0&0&0&0&0&0&0&\dots&-t_{a}&t_{a}&0&\dots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\\
\end{pmatrix}\;,~{}~{}{\bf u}=\begin{pmatrix}R_{x,y,z}\\\ P_{x,y,z}\\\
S_{x,y,z}\\\ R_{x+1,y,z}\\\ P_{x+1,y,z}\\\ S_{x+1,y,z}\\\ R_{x,y+1,z}\\\
P_{x,y+1,z}\\\ S_{x,y+1,z}\\\ S_{x-1,y+1,z}\\\ \vdots\\\ R_{x,y,z+1}\\\
P_{x,y,z+1}\\\ S_{x,y,z+1}\\\ \vdots\end{pmatrix}\;,~{}~{}~{}$ (S1)

where ${\bf u}$ represents an arbitrary state vector.

Second, we present the time-averaged behavior of masses for the parameter
values $t_{a}=1.0,t_{c}=3/8$ and $t_{a}=1.4,t_{c}=3/8$ in Fig. S1a and Fig.
S1b respectively. These two cases correspond to the surface propagation of
masses in a single wave with some definite average velocity as shown in Fig.
2(a-b) of the main text. We present the profile of time averaged masses
$\langle{u_{l}}\rangle_{T}$ in Fig. S1a-b, where boundary localization of
time-averaged masses can be observed for various layers in the $z$ direction.
Here the average is performed over time $t=T$, in which the surface
propagating mass wave completes one round trip around the surface of the
system. Moreover, the above mentioned results are obtained after considering
one hundred realizations of disorder among the payoff matrix elements. Thus,
this indicates that the surface propagation of masses is robustness against
the disorder among the system parameter.

Figure S1: Parameter values are considered as (a) $t_{a}=1.0$, $t_{c}=3/8$,
$\Delta=0.2$, (b) $t_{a}=1.4$, $t_{c}=3/8$, $\Delta=0.2$ for fixed
$t_{b}=1.0$, $N=204$ and $N_{z}=5$. The averaging time $t=T$ is considered in
which the surface propagating mass wave completes one round trip around the
surface of the system. Layer $z=4~{}(5)$ have similar localization as that of
$z=2~{}(1)$ respectively.

##  II Linearized LV equation

In order to construct the linearized LV equation, we follow the scheme
discussed in Ref. [35]. We consider a small perturbation in steady-state
solution ($\partial_{t}{\bf c}=0$ or $\partial_{t}c_{l}=0$) of nonlinear
system such that $u_{l}=c_{l}+\delta{u_{l}}$. Antisymmetric Lotka-Volterra
equation then results in,

$\displaystyle\partial_{t}c_{l}+\partial_{t}\delta{u_{l}}$ $\displaystyle=$
$\displaystyle(c_{l}+\delta{u_{l}})\sum_{m}\alpha_{lm}(c_{m}+\delta{u_{m}}),$
$\displaystyle\partial_{t}\delta{u_{l}}$ $\displaystyle=$ $\displaystyle
c_{l}\sum_{m}\alpha_{lm}c_{m}+c_{l}\sum_{m}\alpha_{lm}\delta{u_{m}}+\delta{u_{l}}\sum_{m}\alpha_{lm}c_{m}+\delta{u_{l}}\sum_{m}\alpha_{lm}\delta{u_{m}}$
$\displaystyle\partial_{t}\delta{u_{l}}$ $\displaystyle=$ $\displaystyle
c_{l}\sum_{m}\alpha_{lm}\delta{u_{m}},$ (S2)

where $\sum_{m}\alpha_{lm}c_{m}=0$ and $\delta{u_{l}}\delta{u_{m}}\approx 0$
under linear approximation. We multiply Eq. (S2) by $i$ on both sides and
obtain,

$\displaystyle i\partial_{t}\delta{u_{l}}$ $\displaystyle=$
$\displaystyle\frac{i}{M}\sum_{m}\alpha_{lm}\delta{u_{m}},$ (S3)

for $c_{l}=1/M$. We recognize the Hamiltonian ${\bf H}=i{\bf A}$ upto a
multiplicative factor of $1/M$ and $\delta{u_{l}}$ is the solution of this
linear Schrödinger (Schrödinger-like) equation such that $|\delta{u_{l}}|^{2}$
describes the probability density.

Figure S2: Dynamics of initial state in linear system is presented for
parameter values (a) $t_{a}=1.0$, $t_{c}=3/8$, (b) $t_{a}=1.4$, $t_{c}=3/8$,
(c) $t_{a}=1.0$, $t_{c}=5/8$ and (d) $t_{a}=3.0$, $t_{c}=4/10$ for
$N_{z}=5,N=52$ and $t_{b}=1.0$. Panel (e,f) shows the dynamics and velocity of
masses in nonlinear system as we vary the perturbation $\lambda$. Parameter
values are taken as $t_{a}=1.0,t_{c}=3/8$ in (e,f).

We computationally solve this linear equation Eq. (S3) for our system
described by ${\bf H}=i{\bf A}$ where ${\bf A}$ is given by Eq. (2) of the
main text. Before explaining the computational results, we address the
question regarding the initial state of the linear system at $t=0$. We
consider the initial state of the system given by $\delta{\bf
u}(t=0)=\sum_{p}\mid\Psi_{\rm surface}^{p}\rangle\langle\Psi_{\rm
surface}^{p}\mid{\bf u}\rangle$, where
$u_{l}=\sum_{j=1}^{N_{z}}\delta_{l,jN}+1/M$ and $\mid\Psi_{\rm
surface}^{p}\rangle$ is the $p^{th}$ surface state of the linear system ${\bf
H}$. By construction, the initial state $\delta{\bf u}(t=0)$ is localized at
the surface of the $3$D system. It is important to highlight that the
$\mid\Psi_{\rm surface}\rangle$ are either right or left moving surface states
such that time evolution results in unidirectional propagation of state [46].
Furthermore. it is worth noting that $\mid{\bf u}\rangle$ is similar to that
considered for the nonlinear system in the main text.

We now turn our attention towards the dynamics of the linear system whose
results are presented in Fig. S2. To visualize the dynamics of the state at
the surface of the system bounded by path $\tilde{s}$ (as shown in Fig. 1 of
the main text), we present the probability density of the time evolving state
in Fig. S2. It is observed that state evolves unidirectionally along the path
$\tilde{s}$ at the surface of the system. For parameter values of
$t_{a}=1.0,t_{c}=3/8$ and $t_{a}=1.4,t_{c}=3/8$, it is observed that the state
evolves as a single wave of probability density as shown in Fig. S2a and Fig.
S2b with an average velocity of $\nu\approx 0.0080$ and $\nu\approx 0.0090$
respectively. However, for parameter values $t_{a}=1.0,t_{c}=5/8$, the
probability density wave split into two, which have two distinct propagation
velocities $\nu_{1}\approx 0.0094$ and $\nu_{2}\approx 0.0052$
($\nu_{1}/\nu_{2}\approx 2$) as shown in Fig. S2c. Finally, Fig. S2d shows the
trivial scenario where probability density is delocalized and does not have
chiral propagation. Qualitatively these dynamical features are the same as
observed for the parent nonlinear system. Quantitatively, however, the
propagation velocities are different by roughly one order of magnitude. This
may be related to a self-trapping mechanism for chiral solitons [46].

To understand the above-mentioned magnitude difference in the wave propagation
velocities between nonlinear and linear systems, we consider the nonlinear
system with a varying initial state. We consider the initial state where mass
on each site $l$ is given as
$u_{l}(t=0)=\sum^{N_{z}}_{j=1}{\lambda}\delta_{l,jN}+1/M$, and $\lambda\leq 1$
is a positive real number that tunes the magnitude of perturbation. By
decreasing $\lambda$ from unity, one can move closer to the steady-state of
the system (in the limit of $\lambda=0$ one approaches the steady-state (Nash
equilibrium state) ) of the nonlinear system. In our system, $\lambda$
controls the interaction between $m^{th}$, and $(m-1)^{th}$ nodes which in
turns affect the velocity of mass propagation. We consider various values of
$\lambda$ and present our results in Fig. S2e-f. In Fig. S2e, we show that the
masses propagate at the surface of $3$D nonlinear system with decreasing
velocities, if we tune the initial state closer to the steady-state by
decreasing the parameter $\lambda$. Compared to $\lambda=1$, the velocity of
mass propagation is much slower for $\lambda=10^{-1}$ and $\lambda=10^{-2}$.
In fact, from Fig. S2f, one can observe that for $\lambda=10^{-3}~{}({\rm
log}(\lambda)\approx-6.9)$ the propagation velocity is $\nu\approx 0.009$, in
agreement with the velocity of probability density $\nu=0.0080$ in linear
system for parameter values $t_{a}=1.0,t_{c}=3/8$. These results strengthen
that main topological features of the linear system are highly useful to
understand the dynamics at the system boundary, even when the dynamics is far
from the linear regime.

##  III Energy spectrum and surface state Fermi arcs

From equation ${\bf H}=i{\bf A}$, we consider periodic boundary conditions in
all spatial directions followed by Fourier transformation from position to
momentum space given as,

$\displaystyle\begin{aligned}
C_{x,y,z}=\frac{1}{\sqrt{N_{x}N_{y}N_{z}}}\sum_{k_{x},k_{y},k_{z}}C_{k_{x},k_{y},k_{z}}e^{-i(k_{x}ax+k_{y}ay+k_{z}az)}\;,\\\
C^{\dagger}_{x,y,z}=\frac{1}{\sqrt{N_{x}N_{y}N_{z}}}\sum_{k_{x},k_{y},k_{z}}C^{\dagger}_{k_{x},k_{y},k_{z}}e^{i(k_{x}ax+k_{y}ay+k_{z}az)}\;,\end{aligned}$
(S4)

where $C\in\\{\text{R},\text{P},\text{S}\\}$, and lattice constant $a=1$. Bulk
Hamiltonian $H({\bf k})$ of the system is given by Eq. (4) of the main text
such that ${\bf H}=\sum_{\bf k}\Psi^{\dagger}_{\bf k}H({\bf k})\Psi_{\bf k}$
and $\Psi^{\dagger}_{\bf k}=(\text{R}^{\dagger}_{\bf
k},\text{P}^{\dagger}_{\bf k},\text{S}^{\dagger}_{\bf k})$. The traceless
Gall-mann matrices in Eq. (4) of the main text are given as,

$\displaystyle\gamma_{1}=\begin{pmatrix}0&1&0\\\ 1&0&0\\\
0&0&0\end{pmatrix},~{}~{}\gamma_{2}=\begin{pmatrix}0&-i&0\\\ i&0&0\\\
0&0&0\end{pmatrix},~{}~{}\gamma_{3}=\begin{pmatrix}1&0&0\\\ 0&-1&0\\\
0&0&0\end{pmatrix},~{}\gamma_{4}=\begin{pmatrix}0&0&1\\\ 0&0&0\\\
1&0&0\end{pmatrix},~{}~{}\gamma_{5}=\begin{pmatrix}0&0&-i\\\ 0&0&0\\\
i&0&0\end{pmatrix},~{}~{}$ $\displaystyle\gamma_{6}=\begin{pmatrix}0&0&0\\\
0&0&1\\\ 0&1&0\end{pmatrix},~{}~{}\gamma_{7}=\begin{pmatrix}0&0&0\\\ 0&0&-i\\\
0&i&0\end{pmatrix},~{}~{}\gamma_{8}=\frac{1}{\sqrt{3}}\begin{pmatrix}1&0&0\\\
0&1&0\\\ 0&0&-2\end{pmatrix}.$

Figure S3: Band structure and surface states Fermi arcs are shown for
parameter values $t_{a}=1.0,t_{c}=3/8$ in (a, d), $t_{a}=1.4,t_{c}=3/8$ in (b,
e), and $t_{a}=1.0,t_{c}=5/8$ in (c, f) for fixed $t_{b}=1.0$. Red (Green)
points shows the type-I (type-II) Weyl nodes while blue curve represent the
surface state Fermi arcs.

We now diagonalize the bulk Hamiltonian for energy eigenvalues by solving the
polynomial cubic equation $-E^{3}+E\theta_{1}+\theta_{2}=0$ for fixed
$t_{b}=1$, whose solutions give the energy eigenvalues,

$\displaystyle\begin{aligned}
E_{1}&=-\frac{2(3^{1/3})\theta_{2}+2^{1/3}(-9\theta_{3}+\sqrt{-12\theta_{2}^{3}+81\theta_{3}^{2}}~{})^{2/3}}{6^{2/3}(-9\theta_{3}+\sqrt{-12\theta_{2}^{3}+81\theta_{3}^{2}}~{})^{1/3}},\\\
E_{2}&=\frac{(-1)^{1/3}\big{(}2(3^{1/3})\theta_{2}+(-2)^{1/3}(-9\theta_{3}+\sqrt{-12\theta_{2}^{3}+81\theta_{3}^{2}}~{})^{2/3}\big{)}}{6^{2/3}(-9\theta_{3}+\sqrt{-12\theta_{2}^{3}+81\theta_{3}^{2}}~{})^{1/3}},\\\
E_{3}&=\frac{-2(-3)^{2/3}\theta_{2}+(-6)^{1/3}(-9\theta_{3}+\sqrt{-12\theta_{2}^{3}+81\theta_{3}^{2}}~{})^{2/3}}{3(2^{2/3})(-9\theta_{3}+\sqrt{-12\theta_{2}^{3}+81\theta_{3}^{2}}~{})^{1/3}},\end{aligned}$

where $\theta_{1}=t_{a}+2t_{c}\cos(k_{z})$,
$\theta_{2}=3+3\theta_{1}^{2}+2\theta_{1}[\cos(k_{x})+\cos(k_{y})+\cos(k_{x}-k_{y})]$,
and
$\theta_{3}=2(\theta_{1}^{2}-\theta_{1})[-\sin(k_{x})+\sin(k_{y})+\sin(k_{x}-k_{y})]$.
One can solve for the relation $E_{i}=E_{j}$ ($i\neq j$) for the band touching
points which results in Weyl nodes
$(k_{x_{0}},k_{y_{0}},k_{z_{0}})=(-2\pi/3,2\pi/3,\pm\cos^{-1}[\frac{\beta-
t_{a}}{t_{c}}])$, and $(2\pi/3,-2\pi/3,\pm\cos^{-1}[\frac{\beta-
t_{a}}{t_{c}}])$ for $\beta=\\{0.5,2.0\\}$. Fig. S3a, S3b, and S3c, show the
band structure where Weyl nodes can be observed at $E_{1}=E_{2}$ and
$E_{2}=E_{3}$ for fixed $k_{z}=k_{z_{0}}$. Type-I and type-II nature of these
Weyl nodes has already been shown in Fig. 4 of the main text.

Surface state Fermi arcs are shown in Fig. S3d, S3e, and S3f for three
representative values of parameters under periodic boundary conditions in
$x-z$ directions and open boundary conditions in $y$ direction where red
(green) points indicate the type-I (type-II) Weyl nodes. Here, blue curve
represent the points in energy-momenta space where $E_{\rm left}=E_{\rm
right}$, with $E_{\rm left}(E_{\rm right})$ indicating the energy of left
(right) moving surface states as a function of $k_{x}-k_{z}$ momenta under
open boundary conditions in $y$ direction. It is worth noting that in hybrid
Weyl semimetal phase, surface state Fermi arc connects a type-I Weyl node to a
type-II Weyl node of opposite chirality.

Finally, we note that all the Chern numbers for parameter values $t_{a}=3.0$,
$t_{b}=1$, and $t_{c}=4/10$ that correspond to Fig. 2d of the main text are
found to be zero as we sweep through $k_{z}$. Thus, there are no Weyl nodes or
surface states Fermi arcs in this case.

## References

  * Klitzing  _et al._ [1980] K. v. Klitzing, G. Dorda, and M. Pepper, Phys. Rev. Lett. 45, 494 (1980).
  * Thouless  _et al._ [1982] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982).
  * Haldane [1988] F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
  * Kane and Mele [2005] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
  * Bernevig  _et al._ [2006] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, science 314, 1757 (2006).
  * Fu  _et al._ [2007] L. Fu, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 98, 106803 (2007).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Burkov and Balents [2011] A. A. Burkov and L. Balents, Phys. Rev. Lett. 107, 127205 (2011).
  * Wan  _et al._ [2011] X. Wan, A. M. Turner, A. Vishwanath, and S. Y. Savrasov, Phys. Rev. B 83, 205101 (2011).
  * Zyuzin  _et al._ [2012] A. A. Zyuzin, S. Wu, and A. A. Burkov, Phys. Rev. B 85, 165110 (2012).
  * Hosur and Qi [2013] P. Hosur and X. Qi, Comptes Rendus Physique 14, 857 (2013).
  * Huang  _et al._ [2015] S.-M. Huang, S.-Y. Xu, I. Belopolski, C.-C. Lee, G. Chang, B. Wang, N. Alidoust, G. Bian, M. Neupane, C. Zhang, _et al._ , Nature communications 6, 1 (2015).
  * Bomantara  _et al._ [2016] R. W. Bomantara, G. N. Raghava, L. Zhou, and J. Gong, Phys. Rev. E 93, 022209 (2016).
  * Bomantara and Gong [2016] R. W. Bomantara and J. Gong, Phys. Rev. B 94, 235447 (2016).
  * Benalcazar  _et al._ [2017] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Science 357, 61 (2017).
  * Schindler  _et al._ [2018] F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. Parkin, B. A. Bernevig, and T. Neupert, Science advances 4, eaat0346 (2018).
  * Li  _et al._ [2018] L. Li, M. Umer, and J. Gong, Phys. Rev. B 98, 205422 (2018).
  * Ghorashi  _et al._ [2020] S. A. A. Ghorashi, T. Li, and T. L. Hughes, Phys. Rev. Lett. 125, 266804 (2020).
  * Tong  _et al._ [2013] Q.-J. Tong, J.-H. An, J. Gong, H.-G. Luo, and C. H. Oh, Phys. Rev. B 87, 201109 (2013).
  * Li  _et al._ [2018] L. Li, C. H. Lee, and J. Gong, Phys. Rev. Lett. 121, 036401 (2018).
  * Haldane and Raghu [2008] F. D. M. Haldane and S. Raghu, Phys. Rev. Lett. 100, 013904 (2008).
  * Raghu and Haldane [2008] S. Raghu and F. D. M. Haldane, Phys. Rev. A 78, 033834 (2008).
  * Wang  _et al._ [2009] Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljačić, Nature 461, 772 (2009).
  * Peng  _et al._ [2016] Y.-G. Peng, C.-Z. Qin, D.-G. Zhao, Y.-X. Shen, X.-Y. Xu, M. Bao, H. Jia, and X.-F. Zhu, Nature communications 7, 1 (2016).
  * Kariyado and Hatsugai [2015] T. Kariyado and Y. Hatsugai, Scientific reports 5, 1 (2015).
  * Süsstrunk and Huber [2015] R. Süsstrunk and S. D. Huber, Science 349, 47 (2015).
  * Nash  _et al._ [2015] L. M. Nash, D. Kleckner, A. Read, V. Vitelli, A. M. Turner, and W. T. Irvine, Proceedings of the National Academy of Sciences 112, 14495 (2015).
  * Delplace  _et al._ [2017] P. Delplace, J. Marston, and A. Venaille, Science 358, 1075 (2017).
  * Zhou  _et al._ [2018] D. Zhou, L. Zhang, and X. Mao, Phys. Rev. Lett. 120, 068003 (2018).
  * Pedro  _et al._ [2019] R. P. Pedro, J. Paulose, A. Souslov, M. Dresselhaus, and V. Vitelli, Phys. Rev. Lett. 122, 118001 (2019).
  * Prodan and Prodan [2009] E. Prodan and C. Prodan, Phys. Rev. Lett. 103, 248101 (2009).
  * Yamauchi  _et al._ [2020] L. Yamauchi, T. Hayata, M. Uwamichi, T. Ozawa, and K. Kawaguchi, arXiv preprint arXiv:2008.10852 (2020).
  * Knebel  _et al._ [2020] J. Knebel, P. M. Geiger, and E. Frey, Phys. Rev. Lett. 125, 258301 (2020).
  * Yoshida  _et al._ [2021] T. Yoshida, T. Mizoguchi, and Y. Hatsugai, Phys. Rev. E 104, 025003 (2021).
  * Yoshida  _et al._ [2022] T. Yoshida, T. Mizoguchi, and Y. Hatsugai, Scientific reports 12, 1 (2022).
  * Knebel  _et al._ [2013] J. Knebel, T. Krüger, M. F. Weber, and E. Frey, Phys. Rev. Lett. 110, 168106 (2013).
  * Knebel  _et al._ [2015] J. Knebel, M. F. Weber, T. Krüger, and E. Frey, Nature communications 6, 1 (2015).
  * Toupo and Strogatz [2015] D. F. P. Toupo and S. H. Strogatz, Phys. Rev. E 91, 052907 (2015).
  * Geiger  _et al._ [2018] P. M. Geiger, J. Knebel, and E. Frey, Phys. Rev. E 98, 062316 (2018).
  * Goel  _et al._ [1971] N. S. Goel, S. C. Maitra, and E. W. Montroll, Rev. Mod. Phys. 43, 231 (1971).
  * Reichenbach  _et al._ [2006] T. Reichenbach, M. Mobilia, and E. Frey, Phys. Rev. E 74, 051907 (2006).
  * Claussen and Traulsen [2008] J. C. Claussen and A. Traulsen, Phys. Rev. Lett. 100, 058104 (2008).
  * Szolnoki  _et al._ [2014] A. Szolnoki, M. Mobilia, L.-L. Jiang, B. Szczesny, A. M. Rucklidge, and M. Perc, Journal of the Royal Society Interface 11, 20140735 (2014).
  * Note [1] See Supplemental Material at insert link for details.
  * Ezawa [2022] M. Ezawa, arXiv preprint arXiv:2201.01891 (2022).
  * Schnyder  _et al._ [2008] A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Phys. Rev. B 78, 195125 (2008).
  * Ryu  _et al._ [2010] S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. Ludwig, New Journal of Physics 12, 065010 (2010).
  * Semmann  _et al._ [2003] D. Semmann, H.-J. Krambeck, and M. Milinski, Nature 425, 390 (2003).
  * Wang  _et al._ [2014] Z. Wang, B. Xu, and H.-J. Zhou, Scientific reports 4, 1 (2014).
  * Lee  _et al._ [2018] C. H. Lee, S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, and R. Thomale, Communications Physics 1, 1 (2018).


# Overview of Solution Methods for Elliptic Partial Differential Equations on
Cartesian and Hierarchical Grids

Damyn M. Chipman

###### Abstract

Elliptic partial differential equations (PDEs) arise in many areas of
computational sciences such as computational fluid dynamics, biophysics,
engineering, geophysics and more. They are difficult to solve due to their
global nature and sometimes ill-conditioned operators. We review common
discretization methods for elliptic PDEs such as the finite difference, finite
volume, finite element, and spectral methods and the linear systems they form.
We also provide an overview of classic to modern solution methods for the
linear systems formed by these discretization methods. These methods include
splitting and Krylov methods, direct methods, and hierarchical methods.
Finally, we show applications that would benefit from fast and efficient
solvers for elliptic PDEs, including projection methods for the incompressible
Navier-Stokes equations and the shallow water wave equations with dispersive
corrections.

##  1 Introduction

Many physical systems can be described using differential equations. Examples
include Newton’s Laws of Motion, almost all conservation laws (mass, momentum,
energy, etc.), many engineering applications, and many more. To simulate or
model such systems, we often use numerical techniques to discretize and solve
the corresponding problem on computers. As computational resources are finite,
researchers must employ efficient algorithms to solve these problems.

Elliptic partial differential equations are a class of PDEs that explain
global and steady state phenomena. Elliptic PDEs arise in fluid dynamics, heat
transfer, electromagnetism, geophysics, biology, and other application areas.
Many elliptic problems can be solved efficiently, but problems on complex
geometries or complicated meshes are more challenging to solve. Solution
methods for elliptic PDEs is a well studied topic, with many papers, books,
and courses detailing their methods. However, fast and efficient ways to solve
elliptic PDEs is an ongoing field of research. Improvements on classical
methods, as well as new ideas, have been recently introduced.

###  1.1 Outline

We will first address Laplace’s and Poisson’s equation in Section 2 as well as
numerical approaches to solving these problems. In Section 3, we will look at
various solution methods for solving the linear system formed by the
discretization methods discussed in Section 2. Next in Section 4, we look at
how elliptic solvers can be used on unstructured and adaptive meshes. In
Section 5, we give a brief overview of the pressure Poisson equation needed to
maintain divergence free velocity fields and a shallow water equation model
with dispersive corrections. Finally, in Section 6 we look at open research
questions in the area of solvers for elliptic PDEs.

##  2 The Elliptic Partial Differential Equation

Partial differential equations are classified by their highest order
derivative terms. A second order, linear differential equation can be written
as

$\displaystyle A(x,y)\frac{\partial^{2}u(x,y)}{\partial
x^{2}}+B(x,y)\frac{\partial^{2}u(x,y)}{\partial x\partial
y}+C(x,y)\frac{\partial^{2}u(x,y)}{\partial y^{2}}+...$
$\displaystyle...+D(x,y)\frac{\partial u(x,y)}{\partial
x}+E(x,y)\frac{\partial u(x,y)}{\partial y}+F(x,y)u(x,y)+G(x,y)$
$\displaystyle=0.$

Second-order linear PDEs are classified according to the value of the
determinant of this expression:

$\displaystyle B^{2}-4AC$ $\displaystyle<0,\ \ \ \text{Elliptic}$
$\displaystyle B^{2}-4AC$ $\displaystyle=0,\ \ \ \text{Parabolic}$
$\displaystyle B^{2}-4AC$ $\displaystyle>0,\ \ \ \text{Hyperbolic}$

In this overview, we will consider common elliptic PDEs such as Laplace’s
Equation

(1) $\displaystyle\nabla^{2}u(x,y)$ $\displaystyle=0,$

Poisson’s Equation

(2) $\displaystyle\nabla^{2}u(x,y)$ $\displaystyle=f(x,y),$

and the variable coefficient Poisson’s Equation

(3) $\displaystyle\nabla\cdot\Big{(}\beta(x,y)\nabla u(x,y)\Big{)}$
$\displaystyle=f(x,y),$

where $\nabla=(\frac{\partial}{\partial x},\frac{\partial}{\partial y})$,
$\nabla^{2}=\nabla\cdot\nabla=\frac{\partial^{2}}{\partial
x^{2}}+\frac{\partial^{2}}{\partial y^{2}}$, $x,y\in\Omega$, and each with the
appropriate boundary conditions on the domain boundary
$\partial\Omega=\Gamma$. Such boundary conditions (BCs) can either be
Dirichlet (Type-I), Neumann (Type-II), or Robin/Mixed (Type-III) BCs.
Dirichlet problems impose the value of $u$ on the boundaries, Neumann problems
impose the flux or normal gradient $\partial_{n}u$ on the boundaries, while
Robin problems impose a linear combination of Dirichlet or Neumann type BCs.

Although analytical solutions exist for some (simple) variations of the
problems above, we are interested in looking at numerical methods to solving
these equations. To use such numerical methods, we first look at various ways
to discretize the domain of elliptic PDEs. Common ways to numerically solve
elliptic PDEs start with a discrete representation of the domain. We will look
at node based, cell based, element based, and spectral based approaches. These
approaches are called the finite difference method, finite volume method,
finite element method, and the spectral method, respectively. We will refer to
Figure 1 for a diagram of each.

Figure 1: Discretization Methods. Top Left: Finite difference grid based on a
mesh of points or nodes. Top Right: Finite volume mesh made of cells (dashed
blue boxes) with cell averages in the center (red points). Bottom Left: Finite
element mesh where each blue box is an element with a shape function defined
at all points within the element. Bottom Right: Spectral mesh made up of a
tensor of Chebyshev nodes.

###  2.1 Finite Difference

If $\Omega$ is on a logically rectangular domain (i.e. 2D Cartesian plane),
perhaps the simplest discretization of the domain $\Omega$ is by collocating a
mesh of points throughout the domain. Given upper and lower bounds in both
directions, $[x_{l},x_{u}],[y_{l},y_{u}]$, the x- and y-location of each point
can be defined as

(4) $\displaystyle x_{i}$ $\displaystyle=x_{l}+i\Delta x,\ \ \
i=0,1,...,N_{x}-1$ (5) $\displaystyle y_{i}$ $\displaystyle=y_{l}+j\Delta y,\
\ \ j=0,1,...,N_{y}-1$

where $N_{x},N_{y}$ is the number of points in the x- and y-direction and grid
spacing is defined as

(6) $\displaystyle\Delta x=\frac{x_{u}-x_{l}}{N_{x}-1}\ \ \ \Delta
y=\frac{y_{u}-y_{l}}{N_{y}-1}.$

These points are shown in the first plot of Figure 1. This leads to the finite
difference method, where the function to approximate is solved for at each
mesh point. In the finite difference approach, the expressions for the
derivatives in a PDE are replaced with Taylor Series approximations. For
example, a second order accurate, central-difference approximation to the
second derivative can be given as

(7) $\displaystyle\frac{\partial^{2}u}{\partial x^{2}}$
$\displaystyle=\frac{u(x-\Delta x,y)-2u(x,y)+u(x+\Delta x,y)}{\Delta
x^{2}}+\mathcal{O}(\Delta x^{2})$

and if we let $u(x_{i},y_{j})=u_{i,j}$, then we can write this as

(8) $\displaystyle\frac{\partial^{2}u}{\partial x^{2}}$
$\displaystyle=\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{\Delta
x^{2}}+\mathcal{O}(\Delta x^{2}).$

Using these Taylor Series expansions, we can replace the continuous Poisson’s
equation 2 with the discrete system of equations:

(9) $\displaystyle\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{\Delta
x^{2}}+\frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{\Delta y^{2}}$
$\displaystyle=f_{i,j}.$

For each grid point $i,j=0,...,N-1$, this expression forms a linear system
with $(N_{x}-1)\times(N_{y}-1)$ unknowns.

###  2.2 Finite Volume

In a finite volume scheme, the domain is broken up into cells, which can be
structured or unstructured polygons (2D) or polyhedrons (3D). The function to
approximate is solved for in terms of a cell average. Finite volume schemes
are often used to solve equations modeling conservation laws, where the cell
quantity is conserved in time through balancing fluxes (what comes in and out
of the cell boundaries) and the cell source (what the cell is generating or
destroying). In a finite volume scheme, the cell average is updated according
to the integral formulation of the conservation law as

(10) $\displaystyle\frac{\partial}{\partial
t}\int_{\Omega_{i}}\textbf{q}d\Omega_{i}+\int_{\Gamma_{i}}\textbf{F}(\textbf{q})\cdot\hat{n}_{i}d\Gamma_{i}+\int_{\Omega_{i}}\textbf{s}d\Omega_{i}=0,$

where $\Omega_{i}$ and $\Gamma_{i}$ are the cell domains and boundaries,
respectively, q is a vector of quantities in each cell, F is a flux function,
and s is the cell’s source term. In the second plot of Figure 1, the blue
dashed lines are the cell boundaries, and the red points are the cell centers.
For up to second order schemes, the cell average is collocated at the center
of the cell; this changes for higher order schemes.

###  2.3 Finite Element

In a finite element scheme, the domain is broken into elements. These elements
can also be structured or unstructured polygons (2D) or polyhedrons (3D) as in
the finite volume method. The PDE to solve is converted into a weak form by
multiplying the PDE by a test function $v(x,y)$ and integrating over the
domain. For example, converting Poisson’s equation into the weak form looks
like this

(11) $\displaystyle\int_{\Omega}\frac{\partial^{2}u}{\partial
x^{2}}vd\Omega+\int_{\Omega}\frac{\partial^{2}u}{\partial y^{2}}vd\Omega$
$\displaystyle=\int_{\Omega}fvd\Omega$ (12)
$\displaystyle\Rightarrow\int_{\Omega}\frac{\partial u}{\partial
x}\frac{\partial v}{\partial x}d\Omega+\int_{\Omega}\frac{\partial u}{\partial
y}\frac{\partial v}{\partial y}d\Omega-\Big{[}\frac{\partial u}{\partial
x}v\Big{]}\Big{|}_{x_{l}}^{x^{u}}-\Big{[}\frac{\partial u}{\partial
y}v\Big{]}\Big{|}_{y_{l}}^{y^{u}}$ $\displaystyle=\int_{\Omega}fvd\Omega,$

where integration by parts is used to transfer a derivative to the test
function. The idea behind the finite element method is to use shape or basis
functions inside each element (defined over the entire element)

(13) $\displaystyle u(x,y)\approx\bar{u}=\sum_{k=1}^{N_{k}}c_{k}\phi(x,y)$

and to use the same basis $\phi$ as the test function $v$ (this is called the
Galerkin principle). Upon substitution of the above into the weak form, and
integrating the basis function with either quadrature or with analytical
expressions, a linear system is formed for the coefficients $c_{k}$. Often,
the coefficients are the actual function values $c_{k}=u_{i,j}$.

There are several variations of the finite element method. The method
explained above forms a global linear system involving all elements. To avoid
this, mapping a physical element to some reference element, makes the
contributions from a single element only influence itself and its neighbors.
This makes the coefficient matrix in the linear system much more sparse and
structured, making it easier to solve. Additionally, by not imposing
continuity across element boundaries, one arrives at the discontinuous
Galerkin method, where the discontinuities are handled through element fluxes.
By choosing certain shape or basis functions, one can derive different
variations of the finite element method, including the B-spline finite element
method [15] and the spectral element method [26].

###  2.4 Spectral Methods

In spectral methods, the approach is to approximate the solution $u$ as a
linear combination of a finite set of orthogonal basis functions

(14) $\displaystyle
u(x,y)=\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}c_{i,j}\phi_{i}(x)\phi_{j}(y)$

where the coefficients are chosen to minimize a norm, such as the $L^{2}$ norm
of the residual $r(x,y)=\nabla^{2}u(x,y)-f(x,y)$. On a discrete domain, this
is similar to requiring $\nabla^{2}u(x_{i},y_{j})=f(x_{i},y_{j})$ at all
interior grid points. As this acts like an interpolation scheme, increasing
the number of discretization points on with a fixed interval actually leads to
highly oscillatory results. Thus, in spectral methods, it is common to use
grid points that are clustered near the ends of the interval, such as
Chebyshev points, defined on an interval $[a,b]$ as

(15) $\displaystyle x_{i}$
$\displaystyle=a+\frac{1}{2}(b-a)\Big{(}1+\cos\big{(}\pi(1-\frac{i}{N+1})\big{)}\Big{)},i=0,...,N+1.$

These points are shown on the last plot of Figure 1. With a good basis of grid
points, spectral methods can achieve very fast convergence. The linear system
formed by this method will be dense as it functions like a high-order
interpolation scheme. However, as convergence is much faster than high-order
finite difference methods, one can use far fewer points on a grid, so the size
of the system is kept small. If one uses Fourier series as the basis functions
$\phi$, one can accelerate the solution of the linear system using fast
Fourier Transform algorithms. Though, due to the periodic nature of the
Fourier series, it is more difficult to implement non-periodic boundary
conditions ([18], [31]).

###  2.5 Other Methods and Summary

These methods are not all the possible ways to solve an elliptic PDE. Other
methods include using quadrature methods on the integral formulation of the
PDE, as well as collocation methods such as radial basis function
approximations. However, the methods we reviewed show some of the most classic
approaches to solving elliptic PDEs.

##  3 Solution Methods for Elliptic Partial Differential Equations

The discretization methods described above will generally lead to a linear
system of equations. The next step is to take advantage of the structure of
these linear systems to solve the system efficiently. We will focus on finite
difference and finite volume approaches. To generally talk about these
solution methods, we assume that we form the linear system

(16) $\displaystyle\textbf{A}\textbf{u}=\textbf{b}$

where $\textbf{A}\in\mathbb{R}^{N\times N}$ is a coefficient matrix from one
of the discretization methods, $\textbf{u}=u_{i}=u(x_{i})$ for
$i\in\textbf{I}_{x}$ and index set $\textbf{I}_{x}$ for the discrete domain,
and b is a vector with the right hand side information $f$ and encoded
boundary information. The goal here is to solve for u where we take advantage
of the structure of A. We organize the various methods into the following
three categories: 1) iterative methods, 2) direct methods, and 3) hierarchical
methods.

###  3.1 Iterative Methods

Iterative methods start with an initial guess of a solution to 16 and correct
the iterate until convergence to a specified tolerance. The simplest iterative
methods are called splitting methods where the linear system is modified
according to

(17)
$\displaystyle\textbf{A}=\textbf{M}-\textbf{N}\Rightarrow\textbf{M}\textbf{u}=\textbf{N}\textbf{u}+\textbf{b}.$

This suggests the following recursion relationship for the next iteration:

(18) $\displaystyle\textbf{M}\textbf{u}^{k+1}$
$\displaystyle=\textbf{N}\textbf{u}^{k}+\textbf{b}$ (19)
$\displaystyle\textbf{u}^{k+1}$
$\displaystyle=\textbf{M}^{-1}\textbf{N}\textbf{u}^{k}+\textbf{M}^{-1}\textbf{b}.$

The idea is to choose M that captures as much of A as possible, but is still
easy and quick to invert. As A is either banded or sparse, classical splitting
methods for elliptic PDEs split A into
$\textbf{A}=\textbf{D}-\textbf{L}-\textbf{U}$, where D is the diagonal
components of A, and L and U are the lower and upper pieces, respectively.
Classical choices for M and N are summarized in Table 1.

Jacobi | $\textbf{M}=\textbf{D}$ | $\textbf{N}=\textbf{L}+\textbf{U}$
---|---|---
Gauss-Sidel | $\textbf{M}=\textbf{D}-\textbf{L}$ | $\textbf{N}=\textbf{U}$
Successive Over Relaxation | $\textbf{M}=\frac{1}{\omega}(\textbf{D}-\omega\textbf{L})$ | $\textbf{N}=\frac{1}{\omega}\big{(}(1-\omega)\textbf{D}+\omega\textbf{U}\big{)}$
Table 1: Iterative Methods: Splitting Methods

Another class of iterative methods are called Krylov subspace methods. The
Krylov space is defined as

(20)
$\displaystyle\mathcal{K}_{k}=\text{span}\\{\textbf{r}_{0},\textbf{A}\textbf{r}_{0},\textbf{A}^{2}\textbf{r}_{0},...,\textbf{A}^{k-1}\textbf{r}_{0}\\}$

based on the initial residual
$\textbf{r}_{0}=\textbf{b}-\textbf{A}\textbf{u}^{(0)}$. The goal is to take
the next iteration from this particular space. Two common Krylov methods are
the Conjugate Gradient method [13] and the Generalized Minimal Residual
(GMRES) method [29]. In the conjugate gradient method, the approximation is
adjusted by a conjugate direction, or a vector that is conjugate with respect
to A. This vector is called the search direction p and is scaled by $\alpha$,
which is computed by solving a quadratic minimization problem. The GMRES
method builds up an orthogonal matrix Q through a process called the Arnoldi
iteration. The Arnoldi iteration forms
$\textbf{A}=\textbf{Q}\textbf{H}\textbf{Q}^{*}$ for orthogonal matrix Q and
Hessenberg matrix H. The next iteration is found via Qy where y is found from
a least squares problem involving H. These methods are summarized in Table 2.

Conjugate Gradient | $\textbf{u}^{(k+1)}=\textbf{u}^{(k)}+\alpha^{(k)}\textbf{p}^{(k)}$
---|---
GMRES | $\textbf{u}^{(k)}=\textbf{Q}^{(k)}\textbf{y}$
Table 2: Iterative Methods: Krylov Subspace Methods

Most iterative methods are considered “matrix-free” methods. A “matrix-free”
method is a method that does not explicitly form the matrix A, but is rather
applied to a vector. For example, if implementing a Conjugate Gradient method
for a finite difference discretization, one has to compute the product
$\textbf{A}\textbf{u}^{(k)}$. Instead of doing the full matrix-vector
calculation, one can write a function that takes u and returns the 2nd-order,
central difference operator as computed in 9.

As finite difference discretization schemes for elliptic PDEs lead to sparse
matrices, the application of A to a vector can be done in $\mathcal{O}(N)$
operations, where $N$ is the number of unknowns in the vector. Thus, the
performance for most iterative methods is approximately $\mathcal{O}(N\times
N_{iter})$, where $N_{iter}$ is the number of iterations required for a
specified tolerance. However, as $N$ gets larger, often so does $N_{iter}$,
leading to poor scaling, as noted in [23]. In addition, iterative methods may
not always converge for a given initial guess or structure of A, which make
them unfavorable for “black-box” implementations for linear solvers.

###  3.2 Direct Methods

Motivation for direct solvers stems from wanting to improve upon the
disadvantages of iterative methods. Martinsson notes in [25] some advantages
to using direct methods over iterative ones:

  * •

Direct methods can be applied to multiple right-hand side vectors b or
multiple boundary conditions once a factorization or solution operator is
built, whereas iterative methods must be solved anew for each right-hand side
or for different boundary conditions.

  * •

Direct methods can take advantage of “close” matrices (i.e. if we have an
inverse or factorization of A and perturb it by $\epsilon$, we could adjust
the inverse to account for it instead of recompute the inverse).

  * •

Most direct methods can take advantage of fast and efficient algorithms for
matrix factorization such as the singular value decomposition, LU
decomposition, QR decomposition, etc.

We will look at how direct methods are useful as we consider some common
direct methods from [18] and [32].

Many direct methods are based on matrix factorizations. Perhaps the most well-
known is the LU decomposition. LU decomposition factors the coefficient matrix
into a lower and upper triangular matrix: $\textbf{A}=\textbf{L}\textbf{U}$.
The idea is to use Gaussian elimination to eliminate entries below the main
diagonal, and then use back-substitution to solve for each entry in u. In
general, LU decomposition requires $\mathcal{O}(N^{3})$ floating point
operations and thus is impractical for large matrices. There are banded
solvers for Gaussian elimination that can take advantage of the sparsity of a
matrix. The Cholesky decomposition is a variant of Gaussian elimination for
symmetric matrices. Other matrix factorizations include the QR-decomposition,
$\textbf{A}=\textbf{Q}\textbf{R}$, and the singular value decomposition,
$\textbf{A}=\textbf{U}\boldsymbol{\Sigma}\textbf{V}^{*}$.

For a finite difference discretization of elliptic PDEs, A is banded and
sparse, and more efficient algorithms for LU decomposition exist. For 1D
problems, A is diagonally dominant and sparse with the bandwidth (the number
of entries off the main diagonal in a matrix) dependent on the order of the
stencil. For the stencil shown in the second order discretization in 9, A is
tridiagonal. This allows us to use algorithms such as Thomas’s algorithm for
solving a tridiagonal system in $\mathcal{O}(N)$ steps ([14]). In higher
dimensions, block versions of Thomas’s algorithm exist ([28]).

Compared to iterative methods, direct methods will terminate in a finite
number of steps. Direct methods are more fit for “black-box” implementations.
However, because most direct methods need to explicitly form A, they are more
difficult to implement with limited computing resources. Indeed, going to
higher dimensions and higher orders often dramatically increases memory and
compute requirements.

###  3.3 Hierarchical Methods

The methods grouped under hierarchical methods attempt to accelerate some of
the ideas from iterative and direct methods by breaking the problem into a
hierarchy of subproblems. By recursively breaking the original problem into
smaller subproblems, significant improvements can be made in complexity and
performance. We’ll talk about three hierarchical methods here: the multigrid
method, nested dissection, and the Hierarchical Poincaré-Steklov (HPS) method.

####  3.3.1 The Multigrid Method

The multigrid method was introduced by Brandt in [6]. It has been widely used
in various applications and for all types of solution methods. Briggs in [7]
gives an overview and tutorial of the multigrid method.

In the multigrid method, the idea is to use multiple levels of grids and
solution methods on each level to solve a larger problem. To look at the
multigrid method, we define the error in the linear system to be
$\textbf{e}=\textbf{u}^{(k)}-\textbf{u}_{exact}$ (the difference between the
exact solution and the $kth$ iteration). After a few iterations of a smoother
(typically Jacobi’s method), because of the local averaging, any high-
frequency error is quickly dampened away. What takes longer to eliminate is
the low-frequency error associated with the global problem. By coarsening the
grid, the lower frequency error is dampened quicker. Multigrid combines the
ability of iterative methods to locally reduce error and a coarsening grid
technique to accelerate convergence.

In the multigrid method, one starts with the solution on a fine grid, for
example, with grid spacing $h$, and performs a few iterations of a smoother.
After a few iterations, the $h$-level grid is projected onto a coarser
$2h$-level grid. On this level, one performs a few more iterations of a
smoother. Because the grid is coarser, this step is faster. Again, after a few
iterations, the solution is projected onto an even coarser $4h$-level grid and
a few more iterations are performed. This is done a specified number of times,
and then the solution is interpolated back up the levels to the original grid.
This is shown in Figure 2.

Figure 2: The 1D Multigrid Method. Starting with the finest grid, perform a
few iterations of an iterative method. This is called relaxing the solution.
Then project onto a coarser grid, and relax again. Do this down the levels in
the grid to a desired precision. Once the solution is converged to on the
coarsest level, interpolate back up the levels to obtain the solution on the
finest level.

In what is called the “full multigrid” method, the process starts on the
coarsest level instead of the finest. The solution is relaxed on this level,
and then interpolated down onto a finer mesh. The interpolated solution is
used as an initial guess for solving the problem on the finer mesh. The
relaxed solution on a coarser grid is often an ideal initial guess for the
problem on the finer mesh, resulting in quick convergence on that level.

The multigrid method allows one to accelerate an iterative solver. Multigrid
methods are very effective as pre-conditioners for iterative methods. This is
the general pattern of hierarchical methods: the ability to use classical
iterative and direct methods on smaller grids where they perform well, and
then “scale” them up to larger problem sizes.

####  3.3.2 Nested Dissection

The nested dissection method formulated by George in [11] is a direct method
that builds upon Gaussian elimination for problems on a grid. It is also the
basis for forming what is called the multifrontal method. By taking advantage
of the ordering of points on a grid, one can permute A to first eliminate
points that split the mesh into two unconnected meshes. This permutation takes
the form of $\textbf{P}^{*}\textbf{A}\textbf{P}$, where the goal is to form P
to reorganize A in a way that eliminates points in an efficient manner.

To detail nested dissection better, consider an $N\times N$ mesh such as a
finite difference mesh discussed in Section 2. Assume that the initial
ordering of points corresponds to an index set that iterates over the points
in the mesh row-by-row. Now, the idea of nested dissection is to reorganize
the points such that we first eliminate points down the middle of the mesh
(i.e. the points at $x=0$ in the first plot of Figure 1). Once these points
are solved for, it spilts the mesh into two equally sized pieces that are
disconnected. Each of the disconnected meshes now are smaller, and thus easier
to solve. This idea of splitting the mesh into disconnected pieces by first
eliminating points along an interface can be recursively applied to each
split. This means that after dividing the mesh into two, one can divide those
two meshes into four meshes, and so on. Recursive splitting is a common
characteristic in hierarchical methods.

Nested dissection was first introduced by George in [11], and further
generalized by Lipton et al. in [19]. Martinsson has a tutorial on nested
dissection in [23]. In fact, nested dissection served as motivation for
another hierarchical method proposed by Martinsson and Gillman called the
Hierarchical Poincaré-Steklov method.

####  3.3.3 The Hierarchical Poincaré-Steklov Method

The work done by Gillman and Martinsson in [25], [21], and [12] (with a
practical tutorial found in [22]) culminate in what they call the Hierarchical
Poincaré-Steklov (HPS) method. It is a direct solver for elliptic PDEs that is
based on a binary tree of rectangular patches where the solution operator to A
is built by recursively merging child patches. Like direct methods, the HPS
method involves a factorization step and a solve step. The chief advantage of
the HPS method over other direct methods is that it does not require the
explicit formulation and storage of A.

The HPS method starts with an original problem domain, and recursively divides
the domain in half. This creates a binary tree of patches as shown in Figure
4. Once the domain has been decomposed into this tree of patches, two
operators are defined on the lowest level, called the leaf level. These
operators are the solution operator S and Dirichlet-to-Neumann (DtN) operator
T. The solution operator maps boundary data to solution data on the interior
of the patch (i.e. solves the local boundary value problem), and the DtN
operator maps Dirichlet data on the boundary to Neumann data on the boundary.
These operators can be formed using any elliptic PDE solver, including fast
solvers like spectral methods. After forming these operators, the next step is
to recursively merge each sibling patch up the tree. the merge step is
demonstrated in Figure 3. This results in a global solution operator that can
be stored and used multiple times (at different time steps or with varying
boundary conditions, etc.), and is similar to a direct method matrix
factorization. The final step is applying the solution operator to each level
down the tree to obtain the solution everywhere in the domain. This step is
just a matrix-vector multiplication and is very fast.

Similar to other direct methods, the HPS method forms an in-memory solution
operator that can be applied to several right-hand side vectors. This property
makes it ideal for problems where several elliptic solves are necessary. While
most iterative methods have better asymptotic performance than other direct
methods, the HPS method can be accelerated using hierarchically block
seperable (HBS) matrix algebra to achieve near linear asymptotic performance
([12]).

Figure 3: HPS Merge Operation. The merged patch $\Omega_{\tau}$ is the union
of children $\Omega_{\alpha}$ and $\Omega_{\beta}$, i.e.
$\Omega_{\tau}=\Omega_{\alpha}\cup\Omega_{\beta}$. Red, green, and blue nodes
correspond to index sets $\textbf{I}_{1}$, $\textbf{I}_{2}$, and
$\textbf{I}_{3}$, respectively. The merge operation eliminates the nodes on
the interface of the children patches. Figure 4: HPS Solve Stage. Once
$\textbf{S}_{0}$ is formed, apply it to the top level Dirichlet data to get
boundary (and solution) data on the interface of the children. Apply the patch
solution operator down the tree until each leaf has it’s local boundary
information. Then apply the solution operator to get the solution data in the
interior of each leaf.

##  4 Adaptive and Unstructured Mesh Methods for Elliptic Partial Differential
Equations

In many applications of elliptic solvers where the dynamics of the problem are
localized (i.e. shock waves and fronts) or the geometry of the domain is
complicated, more complex meshing is used. These types of meshes include
dynamically adaptive meshes, where the mesh is refined/coarsened through time,
and unstructured meshes, where the cells or elements are typically polygons or
polyhedrons and not logically ordered like Cartesian grids. While solving
elliptic PDEs on these types of meshes is more difficult, the methods we have
talked about can be extended to work on complex meshes.

###  4.1 Unstructured Meshes

When working with unstructured meshes, the domain is broken into polygons (2D)
or polyhedrons (3D). Finite difference methods are typically difficult to
implement on unstructured grids, but other methods like finite volume and
finite element methods are often built around these complex geometries. Just
as discretization methods on unstructured grids have been developed, so have
more robust solvers. For example, in [20], Luo et al. present an unstructured
multigrid method for the discontinuous Galerkin method. They use a sequence of
solution approximations of different polynomial orders. The ultraspherical
spectral element method (ultraSEM) presented by Fortunato et al. is another
highly practical extension of the HPS method. In [9], they apply their
ultraspherical spectral element method from [24] to the HPS method. Fortunato
et al. use spectral discretization on mapped quadrilaterals and triangles and
use the same merge process from the HPS method. They show that the HPS method
also works well on mapped, unstructured meshes.

###  4.2 Adaptive Meshes

Hierarchical meshes, like the meshes used in the HPS method, are also well
suited for adaptivity. On adaptive, hierarchical meshes like quadtrees (2D)
and octrees (3D), a patch can be recursively split to a desired level of
refinement. Areas of the domain that have increased dynamics or need higher
refinement (say at boundaries or corners) can be resolved with finer meshes,
while other areas can be of lower resolution. This adaptivity increases
accuracy without significantly increasing runtime.

In [27], Popinet et al. use a quadtree approach to decompose the domain into a
hierarchical grid. In locations where more refinement is necessary, they
refine the mesh more. For example, in [27], they are solving the Serre-Green-
Naghdi model (we will show this as an application in Section 5) which models
shallow water equations. As they are looking at a tsunami traveling across the
surface of the earth, refinement near the location of the wave fronts is
necessary, but without having to have very dense meshing elsewhere. They show
that they can implement the multigrid method on a quadtree based mesh.

Additionally, the HPS method has been modified to be used on an adaptive mesh.
In [10], Geldermans et al., they show how to use the HPS method on a quadtree
of patches. They demonstrate the merge operation between patches on different
levels, as well as discuss proper interpolation schemes for their choice of
mesh (a Chebyshev tensor product). In addition, the global solution operator
formed by this adaptive HPS method can still be applied to multiple right-hand
side problems. The study of fast, direct solvers for elliptic problems on
adaptively refined quadtrees and octrees is an active area of research.

##  5 Applications of Elliptic Partial Differential Equations

We will look at two examples where the use of a fast elliptic PDE solver would
greatly benefit the solution process. These are cases where an elliptic solve
is required multiple times, such as every time step or for multiple problems.

###  5.1 Navier-Stokes: The Projection Method

The incompressible Navier-Stokes equations is a nonlinear PDE that arises from
conservation of momentum and mass in an incompressible fluid. They are
expressed as

(21) $\displaystyle\frac{\partial\textbf{u}}{\partial
t}+\textbf{u}\cdot\nabla\textbf{u}$ $\displaystyle=-\frac{1}{\rho}\nabla
P+\nu\nabla^{2}\textbf{u}.$

In [8], Chorin first computes an intermediate velocity $\textbf{u}^{*}$ by
time stepping and ignoring the pressure term

(22) $\displaystyle\frac{\textbf{u}^{*}-\textbf{u}^{n}}{\Delta t}$
$\displaystyle=-\textbf{u}^{n}\cdot\nabla\textbf{u}^{n}+\nu\nabla^{2}\textbf{u}^{n}$

and then “projecting” the intermediate velocity to the next time step via

(23) $\displaystyle\textbf{u}^{n+1}=\textbf{u}^{*}-\frac{\Delta t}{\rho}\nabla
P^{n+1}$

In order to solve for the right-hand side of the second step, the pressure
field at the $n+1$ time step must be known. This is found by solving the
following Poisson’s equation arising from taking the divergence of the second
step and using the conservation of mass to eliminate the
$\nabla\cdot\textbf{u}^{n+1}$ term:

(24) $\displaystyle\nabla^{2}P^{n+1}$ $\displaystyle=\frac{\rho}{\Delta
t}\nabla\cdot\textbf{u}^{*}$

This projection method requires an elliptic solve at every time step to solve
for the pressure field. In applications like this where an elliptic solve is
required multiple times, direct methods where one can pre-compute the solution
operator and then apply it each time step greatly speeds up the computation.

###  5.2 The Shallow Water Equations and the Serre-Green-Naghdi Model

Figure 5: Shallow Water Equation Reference

The shallow water equations (SWE) form a model for free surface flow derived
from potential theory. The key assumption in the SWEs is that the horizontal
scale is much larger than the vertical scale. Because a “shallow limit” is
assumed, the velocity field in the vertical direction is a depth-averaged
velocity. These types of equations are used to model tsunami waves and debris
flow, among other things.

Extensions to this model done by Bonneton et al. in [5], [16], and [17]
introduce a dispersive term correction in the form of a source term, leading
to

(25) $\displaystyle\frac{\partial\textbf{q}}{\partial
t}+\nabla\cdot\textbf{F}(\textbf{q})=\textbf{b}+\textbf{s}$

where

(26) q $\displaystyle=\begin{bmatrix}h\\\ hu\\\ hv\\\ \end{bmatrix}$ (27)
$\displaystyle\textbf{F}(\textbf{q})$ $\displaystyle=\begin{bmatrix}hu&hv\\\
hu^{2}+\frac{1}{2}gh^{2}&huv\\\ huv&hv^{2}+\frac{1}{2}gh^{2}\\\ \end{bmatrix}$
(28) b $\displaystyle=\begin{bmatrix}0\\\ -gh\frac{\partial b}{\partial x}\\\
-gh\frac{\partial b}{\partial y}\\\ \end{bmatrix}$ (29) s
$\displaystyle=\begin{bmatrix}0\\\
\frac{\alpha-1}{\alpha}gh\frac{\partial\eta}{\partial x}+d_{x}\\\
\frac{\alpha-1}{\alpha}gh\frac{\partial\eta}{\partial y}+d_{y}\\\
\end{bmatrix}.$

The water column height is denoted by $h(x,y)$, with surface velocity
$(u(x,y),v(x,y))$, $b(x,y)$ is a bathymetry term (elevation), $\eta(x,y)$ is
the free surface height and dispersive term correction $(d_{x},d_{y})$. Figure
5 shows the typical setup for this kind of problem.

In order to solve for the dispersive source term, one must solve the following
elliptic equation

(30)
$\displaystyle\big{(}\textbf{I}+\alpha\textbf{T}_{diag}^{b}\big{)}\textbf{d}=\textbf{f}(h,\textbf{u},b,\eta)$

where the details of the right-hand side function can be found in [17]. The
operator $\textbf{T}_{diag}^{b}$ is an elliptic-like operator that does not
change with time. So, for a static mesh, the dispersive operator can be
factorized via a direct method and applied at each time step. This reduces the
need to solve the elliptic PDE every time step to just an application of the
solution operator every time step. If the mesh is adapted, the solution
operator would need to be reformed. However, direct methods like the HPS
method can allow for a local adjustment to the solution operator, further
accelerating the solution.

##  6 Open Research Questions

There remain topics for potential research in the area of fast solvers for
elliptic PDEs. In this paper, we looked at classic approaches to approximating
solutions of PDEs including finite difference, finite volume, finite element,
and spectral methods. Each of these approaches have been analyzed and
expanded, with high order methods and fast and parallel algorithms being the
primary focus of recent advancements. Although other discretization methods
such as collocation methods and meshless methods remain an ongoing research
topic, we will look at some open research topics in the area of solution
methods, as it is an area with more potential questions and applications. We
have already addressed classic solution methods, which we categorized as
iterative, direct, or hierarchical. Iterative methods such as splitting or
Krylov methods and direct methods such as matrix factorizations are well
established in the literature, with textbooks and courses addressing these
topics. Indeed, these topics form the foundation for numerical analysis for
PDEs. The methods we classified as hierarchical methods (i.e., the multigrid
method, nested dissection, and the HPS method) are more open topics. Multigrid
methods remain a popular solution technique for a variety of discretization
methods, as well as on adaptive meshes ([2, 30]). The HPS method has also been
introduced and is being improved upon, as in ([9, 10, 12, 22]). However, open
questions exist with regard to parallel implementations and further
applications of these methods.

Applying the multigrid or HPS method to a particular problem opens up a wide
range of new ideas. For example, in the area of computational fluid dynamics
(CFD), the multigrid method has been used to help solve the Navier-Stokes
equations ([2]). Wherever the multigrid method has been applied, the HPS
method can also be applied for some competitive results. It would be of
interest to see how the HPS method and multigrid method compare when applied
to the same problem with the same discretization scheme. For example, the
multigrid method has been applied to the Serre-Green-Naghdi model with a
quadtree discretization ([27]). The HPS method works well on quadtrees and can
also be applied here, with the advantage of a direct method for faster
solutions.

Both the multigrid and HPS method have potential for high scalability. The
multigrid method has been implemented in parallel in a variety of papers
([4]). The main approach is to use a domain decomposition approach to divide
the domain among processors. The HPS method, however, is fairly new and has
limited study in parallelization. Beams et al. in [3] describe a shared-memory
implementation of the HPS method. However, they do not address distributed
memory paradigms or implementations on a GPU. Additionally, while their
results show promising scaling results, they use a relatively small number of
cores (up to 28 physical cores). It would be of interest to see how the HPS
method scales with communication cost and on much larger machines using a
distributed memory paradigm or heterogeneous architecture. The motivation for
highly scalable algorithms comes from the advances in high performance
computing. Modern supercomputers have millions of cores, and their speed comes
from the distribution of work. State of the art supercomputers need to use
advanced algorithms that scale well enough to take advantage of the compute
power. This is the mission of the United States’ Department of Energy’s (DoE)
Exascale Computing Project (ECP); to combine exascale-capable machines with
highly scalable algorithms for “codesign, modeling and simulation, data
analytics, machine learning, and artificial intelligence” ([1]).

##  7 Conclusion

The discretization of elliptic partial differential equations leads to linear
systems that can be solved with efficient solvers, taking advantage of the
structure of the problem. Discretization methods include the finite difference
method, finite volume method, finite element method, spectral methods, and
more. These discretizations lead to linear systems of the standard form
$\textbf{A}\textbf{u}=\textbf{f}$, where A is sparse and structured. Efficient
methods such as splitting and Krylov methods, matrix factorizations, and
hierarchical domain decomposition methods can take advantage of the structure
of the linear system. These solution methods can also be expanded upon in
order to achieve near linear performance and are well suited for high-
performance applications.

In addition, we presented various areas of study that could greatly benefit
from an efficient, fast, and scalable solver for elliptic PDEs. This includes
many applications in computational fluid dynamics like the incompressible
Navier-Stokes and the shallow water wave equations. Using projection methods
for the Navier-Stokes equations and the Serre-Green-Naghdi dispersive model
leads to elliptic terms that must be solved at each time step. Fast, direct
methods like the ones currently being studied and implemented can accelerate
these areas as well as any area that solves an elliptic partial differential
equation.

## References

  * [1] April 2021, https://www.exascaleproject.org.
  * [2] R. Babich, J. Brannick, R. C. Brower, M. Clark, T. A. Manteuffel, S. McCormick, J. Osborn, and C. Rebbi, Adaptive multigrid algorithm for the lattice Wilson-Dirac operator, Phys. Rev. Lett., 105 (2010), p. 201602.
  * [3] N. N. Beams, A. Gillman, and R. J. Hewett, A parallel shared-memory implementation of a high-order accurate solution technique for variable coefficient Helmholtz problems, Comput. Math. Appl., 79 (2020), pp. 996–1011.
  * [4] B. Bergen, T. Gradl, F. Hulsemann, and U. Rude, A massively parallel multigrid method for finite elements, Comput. Sci. Eng., 8 (2006), pp. 56–62.
  * [5] P. Bonneton, F. Chazel, D. Lannes, F. Marche, and M. Tissier, A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model, J. Comput. Phys., 230 (2011), pp. 1479–1498.
  * [6] A. Brandt, Multi-level adaptive solutions to boundary-value problems, Math Comput., 31 (1977), pp. 333–390.
  * [7] W. L. Briggs, V. E. Henson, and S. F. McCormick, A multigrid tutorial, SIAM, 2000.
  * [8] A. J. Chorin, The numerical solution of the Navier-Stokes equations for an incompressible fluid, B. Am. Math Soc., 73 (1967), pp. 928–931.
  * [9] D. Fortunato, N. Hale, and A. Townsend, The ultraspherical spectral element method, J. Comput. Phys., (2020), p. 110087.
  * [10] P. Geldermans and A. Gillman, An adaptive high order direct solution technique for elliptic boundary value problems, SIAM J. Sci. Comput., 41 (2019), pp. A292–A315.
  * [11] A. George, Nested dissection of a regular finite element mesh, SIAM J Numer. Anal., 10 (1973), pp. 345–363.
  * [12] A. Gillman and P.-G. Martinsson, A direct solver with O(N) complexity for variable coefficient elliptic PDEs discretized via a high-order composite spectral collocation method, SIAM J. Sci. Comput., 36 (2014), pp. A2023–A2046.
  * [13] M. R. Hestenes, E. Stiefel, et al., Methods of conjugate gradients for solving linear systems, vol. 49, NBS Washington, DC, 1952.
  * [14] N. J. Higham, Accuracy and stability of numerical algorithms, SIAM, 2002.
  * [15] P. Kagan, A. Fischer, and P. Z. Bar-Yoseph, New B-Spline Finite Element approach for geometrical design and mechanical analysis, Int. J. Numer. Meth. Eng., 41 (1998), pp. 435–458.
  * [16] D. Lannes and P. Bonneton, Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation, Phys. Fluids, 21 (2009), p. 016601.
  * [17] D. Lannes and F. Marche, A new class of fully nonlinear and weakly dispersive Green–Naghdi models for efficient 2D simulations, J. Comput. Phys., 282 (2015), pp. 238–268.
  * [18] R. J. LeVeque, Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems, SIAM, 2007\.
  * [19] R. J. Lipton, D. J. Rose, and R. E. Tarjan, Generalized nested dissection, SIAM J Numer. Anal., 16 (1979), pp. 346–358.
  * [20] H. Luo, J. D. Baum, and R. Löhner, A p-multigrid discontinuous Galerkin method for the Euler equations on unstructured grids, J. Comput. Phys., 211 (2006), pp. 767–783.
  * [21] P. Martinsson, A direct solver for variable coefficient elliptic PDEs discretized via a composite spectral collocation method, J. Comput. Phys., 242 (2013), pp. 460–479.
  * [22] P. Martinsson, The hierarchical Poincaré-Steklov (HPS) solver for elliptic PDEs: A tutorial, arXiv preprint arXiv:1506.01308, (2015).
  * [23] P.-G. Martinsson, Fast direct solvers for elliptic PDEs, SIAM, 2019\.
  * [24] S. Olver and A. Townsend, A fast and well-conditioned spectral method, SIAM Rev., 55 (2013), pp. 462–489.
  * [25] V. R. P. Martinsson, A fast direct solver for boundary integral equations in two dimensions, J. Comput. Phys., (2004).
  * [26] A. T. Patera, A spectral element method for fluid dynamics: laminar flow in a channel expansion, J. Comput. Phys., 54 (1984), pp. 468–488.
  * [27] S. Popinet, A quadtree-adaptive multigrid solver for the Serre–Green–Naghdi equations, J. Comput. Phys., 302 (2015), pp. 336–358.
  * [28] A. Quarteroni, R. Sacco, and F. Saleri, Numerical mathematics, vol. 37, Springer Science & Business Media, 2010.
  * [29] Y. Saad and M. H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comp., 7 (1986), pp. 856–869.
  * [30] M. Thompson and J. H. Ferziger, An adaptive multigrid technique for the incompressible Navier-Stokes equations, J. Comput. Phys., 82 (1989), pp. 94–121.
  * [31] A. Townsend and S. Olver, The automatic solution of partial differential equations using a global spectral method, J. Comput. Phys., 299 (2015), pp. 106–123.
  * [32] L. N. Trefethen and D. Bau III, Numerical linear algebra, vol. 50, SIAM, 1997.