pnasresearcharticle Datseris Dynamical systems are ubiquitous to model natural
phenomena ranging from power grids, climate, ecosystems, and more. Global
stability analysis is the study of their attracting states over the full state
space. It goes beyond (linear) bifurcation analysis and provides insight on
the resilience of the states, i.e., how close they are to “tip” to a
different, perhaps undesirable stable state such as population death in
ecosystems or circulation shutdown in climate. In this work we develop an
accurate and flexible framework for global stability analysis of arbitrary
dynamical systems, that we believe will accelerate multistability and tipping
point research in many fields. We also provide a software implementation that
is fast and easy to use in the open source DynamicalSystems.jl library. Author
contributions: G.D. coordinated the project. All authors contributed in
compiling the results, creating the code, writing and revising the manuscript.
Authors declare no conflicts of interest. 1To whom correspondence should be
addressed. E-mail<EMAIL_ADDRESS>
# Framework for global stability analysis of dynamical systems

George Datseris  Department of Mathematics and Statistics, University of
Exeter, Exeter, United Kingdom  Kalel Luiz Rossi  Theoretical Physics/Complex
Systems, ICBM, Carl von Ossietzky University Oldenburg, Oldenburg, Lower
Saxony, Germany  Alexandre Wagemakers  Nonlinear Dynamics, Chaos and Complex
Systems Group, Departamento de Física, Universidad Rey Juan Carlos, Madrid,
Spain

###### Abstract

Dynamical systems, that are used to model power grids, the brain, and other
physical systems, can exhibit coexisting stable states known as attractors. A
powerful tool to understand such systems, as well as to better predict when
they may “tip” from one stable state to the other, is global stability
analysis. It involves identifying the initial conditions that converge to each
attractor, known as the basins of attraction, measuring the relative volume of
these basins in state space, and quantifying how these fractions change as a
system parameter evolves. By improving existing approaches, we present a
comprehensive framework that allows for global stability analysis on any
dynamical system. Notably, our framework enables the analysis to be made
efficiently and conveniently over a parameter range. As such, it becomes an
essential complement to traditional continuation techniques, that only allow
for linear stability analysis. We demonstrate the effectiveness of our
approach on a variety of models, including climate, power grids, ecosystems,
and more. Our framework is available as simple-to-use open-source code as part
of the DynamicalSystems.jl library.

###### keywords:

critical transitions $|$ multistability $|$ tipping $|$ continuation $|$
global stability

###### doi:

www.pnas.org/cgi/doi/10.1073/pnas.XXXXXXXXXX

This manuscript was compiled on

Multistable dynamical systems exhibit two or more co-existing stable states,
formally called _attractors_. Multistability is ubiquitous in nature and in
mathematical models (1), with examples ranging from power grids (2, 3, 4), the
climate (5, 6), ecosystems like the Amazon rain forest (7, 8), the brain and
neuronal circuits therein (9, 10, 11), or metabolic systems (12, 13, 14). Some
attractors of these systems can be desirable, such as synchronized
oscillations in power grids, crucial for their proper functioning (15). But
they can also be undesirable, as for example the extinction of a certain
species in ecological models or the collapse of circulation in climate models
(16). In a multistable system, the attractor at which the system ends up
depends on the initial conditions, but perturbations of the state may enforce
switching between attractors, a phenomenon called “tipping” (17, 1, 8).
Alterations in the parameters of a dynamical system can trigger tipping.
Hence, it becomes important to evaluate how “resilient” attractors are to
perturbations, either to parameters or to the system’s variables. This is a
crucial problem of practical importance in several areas of research (8, 18).

A traditional solution to this problem is the _continuation-based bifurcation
analysis_ (CBA). It identifies fixed points, limit cycles (under some
requirements), and describes their linear stability dependence on a system
parameter. One major downside is that, by definition, it cannot be applied to
chaotic attractors (19). In the majority of cases, this analysis must be done
numerically via one of several software, e.g., AUTO (20), MATCONT (21), CoCo
(22), or Bifurcationkit.jl (23). The information provided by these frameworks
is useful, but incomplete: rigorously speaking, linear stability only conveys
information about the system’s response to infinitesimally local
perturbations. It cannot yield insight on the response to finite perturbations
in the state space, which are predominant in practice.

For such responses, it is necessary to study the _global stability_ (24) of
the system’s attractors, which involves the nonlinear dynamics over the full
state space of the system111The term “global stability” is also used when a
dynamical system has a single global attractor, which is different from our
use of the term here.. A proxy for global stability of an attractor is the
portion of all possible initial conditions ending up at this attractor, i.e.,
the _fraction_ of the state space that is in the _basin_ of said attractor.
When the state space is infinite, the concept of the state space fraction
becomes a pragmatic one: we need to define a finite-volume box of _physically
plausible_ initial conditions for the system under study, and we are concerned
about the fractions of these plausible initial conditions. In this analysis,
attractors with larger basins fractions are globally more stable because
stronger perturbations are typically needed to switch the system to another
attractor (24). Frequently, this measure is also a much better indicator of
the loss of stability as a system parameter is varied, when compared to the
linear analysis of the system (see e.g., (24) or (19, Chap. 12)).

Analyzing global stability as a function of a parameter demands extensive
effort from researchers, as it requires the creation of algorithms that can
find system attractors and their global stability, “continue” them across a
parameter, and also perform the expensive numerical simulations required for
such algorithms. In the literature, the only framework so far that can aid
this analysis is the _featurizing and grouping_ approach, proposed first by
Gelbrecht et al (25) as MCBB (Monte Carlo Basin Bifurcation Analysis) and then
later very similarly by Stender et al (26) as bSTAB (basin stability). The
method integrates randomly sampled initial conditions (ICs) of a dynamical
system for a preset time span. The trajectories of these ICs are then mapped
onto _features_ , numbers describing the trajectories, such as the mean or
standard deviation of some of the system variables. All the feature vectors
are clustered using the DBSCAN algorithm (27) so that ideally each cluster
corresponds to an attractor of the system. The fractions of ICs in each
cluster approximate the basin fractions and hence the global stability of each
attractor. More details on the method in the Materials and Methods (MatMeth)
§A.

This method works well in a variety of circumstances and can also be applied
across a parameter range. However, it comes with significant downsides. One is
that it is not clear a-priori which features should be chosen to correctly
separate the attractors into clusters, requiring a lot of trial and error.
Another downside is that the method cannot guarantee that the clusters of
features really correspond to unique attractors, and that it is not mixing two
or more attractors together. An alternative method for finding attractors and
their basins of attraction is the _recurrence-based_ approach recently
proposed in Ref. (28). The method locates attractors by finding recurrences in
the system’s state space, assuming the Poincaré recurrence theorem holds for
the system attractors. The input to this method is a state space box, and its
tessellation, defining a grid to search for recurrences. Hence, the method
will only find attractors within the given box, although the box can initially
be arbitrarily large. We describe the method in more detail in MatMeth§B and
provide a comparison between the two techniques in §22.2. The main advantage
of the recurrences method is that it locates the actual system attractors and
only requires as an input a state space box that may contain the attractors.
So far however it has not been clear how to “continue” attractors across a
parameter range with this method.

In this work, in §11.1, we present a novel global stability analysis and
_continuation_ algorithm that utilizes the recurrences-based method for
finding attractors (28). In §11.3 we apply it to exemplary models of climate,
ecosystem dynamics, and more. As detailed in §2, this novel continuation
algorithm is the most accurate in finding the actual attractors of a dynamical
system, the most transparent in matching attractors across parameter values,
and requires the least amount of guesswork from the researcher. We believe
that this novel continuation of global stability, much like global stability
analysis itself (24), is a crucial new tool for the analysis of dynamical
systems. In some cases it supersedes CBA, and in others it can be complemented
by CBA, as we discuss in §22.1.

This continuation method is part of a novel automated framework that performs
global stability analysis and continuation, which we present in §11.2. Our
framework significantly advances existing methodology, including the
featurizing methods, thereby including all upsides of current literature while
addressing most downsides (MatMeth §C and D describe the improvements in
detail). Its design is based on modular components that can be configured or
extended independently. This allows researchers to simply compose the
methodology that is best suited to their problem, and then let an automated
algorithm execute the process. The framework is accompanied by an extensively
tested, well documented, and highly optimized open source software
implementation, part of the DynamicalSystems.jl (29) general purpose library
for nonlinear dynamics (see MatMeth for code example and documentation).

##  1 Results

###  1.1 Novel global stability continuation algorithm

A major contribution of our framework is the novel algorithm for global
stability analysis and continuation that we name _recurrences-based attractor
find-and-match continuation_ , RAFM for short. As illustrated in Fig. 1A, it
works as follows. Step 0: the starting point of the algorithm. Attractors and
their basins, or basins fractions, are already known at a parameter $p=p_{1}$
using the recurrences-based algorithm of Ref. (28). Step 1: new initial
conditions are seeded from the existing attractors. Then, we set the system
parameter to $p=p_{2}$. Step 2: evolve the seeded initial conditions according
to the dynamic rule of the system. The seeds are evolved until they converge
to an attractor using the recurrences-based method (the grid reflects the
tessellation of the state space that is decided by the user; the finer the
grid, the more accurate the results (28)). The main performance bottleneck of
the recurrences-based method is finding the attractors. Once found,
convergence of other initial conditions is generally much faster (28). To
address this in the continuation, we use the observation that, unless a
bifurcation is occurring, attractor size, shape, and position, typically
depend smoothly on $p$. Hence, the seeded initial conditions at each new
parameter will most likely converge the fastest to the new attractors. Step 3:
match attractors in current parameter $p_{2}$ to those in parameter $p_{1}$.
Matching is arguably the most sophisticated part of the algorithm. It is
flexible and attractors can be matched by their “distance” in state space, see
§11.4 for more details. In this illustration the “red” attractor is matched to
the “purple” one of Step 0, while the “yellow” attractor doesn’t map to the
previous “teal” one, because their state space distance is beyond a pre-
defined threshold. Step 4: with the main bottleneck of the algorithm (finding
the attractors) being taken care of, now compute the basins fractions by
formally applying the recurrence-based algorithm of Ref. (28). Importantly,
the algorithm may still find new attractors (here the “dark green” one) that
didn’t exist before. The end result is the system attractors, and their
basins, as functions of the parameter. The attractors and basins are labelled
with the positive integers (enumerating the different attractors), and the
basins always sum to 1. Note that because RAFM works on a parameter-by-
parameter basis, it can be used to perform continuation across any number of
parameters, not just one (we present one here as the simplest conceptual
example).

Figure 1: A: The recurrences-based seed-and-match algorithm for global
stability continuation described in §11.1. B: Schematic illustration of the
modular framework for global stability analysis and continuation described in
§11.2.

###  1.2 Global stability continuation framework

To perform global stability analysis, several tasks need to be taken in
sequence, see Fig. 1B for an overview. We have abstracted and generalized the
tasks to allow researchers different possibilities of how to achieve them.

The first task is the creation of a dynamical system for the global stability
analysis. For our framework, this is achieved for free simply by making the
implementation part of the DynamicalSystems.jl library (29) (see MatMeth E).

The second task is the creation of a mechanism to find attractors and map
initial conditions to them. Possibilities for this mechanism are: (1)
_featurizing and grouping_ initial conditions into attractors (as discussed in
the introduction), (2) finding attractors using the _recurrences algorithm_
(28), or (3) mapping initial conditions to _previously known attractors by
proximity_ : once the evolution of an initial condition comes close enough to
a pre-determined attractor, the initial condition is mapped to that attractor.
Note that in (1) or (2) attractors are found via random sampling in the state
space, and the probability to find an attractor with basin fraction $f$ after
$n$ samples is $1-(1-f)^{n}$.

Mechanism (1) is paired with instructions on how to group features. Currently,
the possibilities are: (a) _clustering_ features into groups using DBSCAN (as
done in MCBB or bSTAB), (b) grouping features by _histograms in the feature
space_ , so that features that end up in the same histogram bin belong to the
same group (novel grouping approach), or (c) mapping features to their
_nearest feature_ in feature space, from a set of pre-defined features (as
done in bSTAB). Having more grouping possibilities than only clustering can be
useful, as discussed in §11.4.

At this point one can analyze global stability at a given parameter, and also
analyze the basin boundaries of attractors for fractal properties. From here,
the third task is to “continue” the attractors and basins across a parameter
range. Our framework has two continuation methods. The first is what has been
employed so far by the MCBB or bSTAB algorithms, but with significantly
increased accuracy (see MatMeth), and extended to allow any kind of
instructions for how to group features. In this approach, trajectories from
the dynamical system are generated by sampling random ICs across all parameter
values of interest. All these trajectories are mapped to features, that are
then grouped, and each group is representing an attractor. The grouped ICs are
then re-distributed into the parameter slices they came from, providing the
fractions of each group at each parameter value. The alternative is the RAFM
algorithm that we described in §11.1.

###  1.3 Application on exemplary systems

In Fig. 2 we apply RAFM on some exemplary systems. We stress that we could
characterize the different attractors in accurate detail because RAFM finds
the actual system attractors, not some incomplete representation of them
(i.e., features used in the featurizing-and-grouping approach).

Figure 2: Basins fraction continuation for exemplary dynamical systems using
the novel recurrences-based continuation algorithm. The fractions of the
basins of attraction are plotted as stacked band plots (hence, summing to 1).
Each color corresponds to a unique attractor that is found and continued (but
not plotted here). Each simulation scanned 101 parameter values, and in each
it sampled randomly 100 initial conditions. The fractions fluctuate strongly
versus parameter not due to lack of convergence, but because the basin
boundaries are fractal in all systems considered. The systems used are: a)
3-dimensional paradigmatic chaotic model by Lorenz (Lorenz84 (30)) with a co-
existence of a fixed point, limit cycle, and chaotic attractor, undergoing a
crisis with the chaotic attractor merging into the limit cycle; b)
33-dimensional climate toy model (31) featuring bistability of chaotic
attractors; c) 3-dimensional multistable cell-division model (32) where each
cell type is considered to be a distinct attractor in the gene activity state
space; d) 9-dimensional model for turbulent shear flow model in which the
fluid between two walls experiences sinusoidal body forces (33); e)
8-dimensional ecosystem competition dynamics model (34) featuring extreme
multistability (due to the number of attractors, we made no effort to label
them further).

###  1.4 Matching and grouping attractors

Traditional CBA has a rigid “matching” procedure: it always matches the next
point found along a “continuation curve” to the previous point. This is often
correct for infinitesimal perturbations of fixed points, but becomes
problematic for global stability analysis, which attempts to find all
attractors in a state space box and then continue them. In this case, matching
attractors from one parameter to the next becomes a crucial part of the
algorithm. For instance, the analysis presented in Fig. 2 is only coherent
because of the powerful matching procedure implemented in our framework.
Without it, the colors would alternate arbitrarily at each parameter value.

In the featurize-and-group algorithm, matching and grouping are the same
process. In RAFM, matching operates on a parameter-by-parameter basis. Each
time the parameter is incremented, and the new attractors are found, a
matching sub-routine is launched. The distance between attractors before and
after the parameter change is estimated, with “distance” being any symmetric
positive-definite function defined on the space of state space sets. By
default the Euclidean distance of the attractor centroids is used, but the
Hausdorff distance (35) is also provided out of the box. After the distance is
computed between all new-old attractor pairs, the new attractor labels are
matched to the previous attractor labels that have the smallest distance to
them, prioritizing pairs with smallest distance. The matching respects
uniqueness, so that an attractor from one parameter cannot be matched with
more than one attractor from another parameter. Additionally, a distance
threshold value can be provided, so that attractors whose distance is larger
than this threshold are guaranteed to get assigned different IDs. Note that in
principle finding the attractors and matching them are two completely
independent processes. If after the continuation process is finished the user
decides that the chosen matching procedure was unsuitable, they can launch a
“re-matching” algorithm with different matching “distance” function, without
having to re-do any computations for finding the attractors or their fractions
(matching only renames attractor labels).

The last thing to highlight in this section is the desirable post-processing
of grouping similar enough attractors. This happens automatically if one uses
the featurize-and-group continuation method. However, taking as an example
Fig. 2(e), the RAFM method finds countless individual attractors. For the
researcher, the individual attractors may be useful for careful analysis, but
it is sometimes desirable to group similar enough attractors. In our framework
it is possible to use exactly the same grouping infrastructure utilized by the
featurizing-and-grouping continuation, but now applied to the outcome of RAFM
as a post-processing step. In Fig. 3 we highlight examples that utilize the
powerful matching and/or grouping components offered by our framework.

Figure 3: Highlights of the matching or grouping components of the framework.
a) Matching attractors of the Hénon map based on their period. In the chosen
parameter range, an attractor is transformed from chaotic, to period 3, 7, and
14. The attractor stays in approximately the same state space location, so
whether we consider the centroid distance or the Hausdorff distance, the
attractor would be matched to itself in all parameter values due to the very
small distance evaluation. Here however we use as distance
$f(A,B)=|\log_{2}(\mathrm{len}(A))-\log_{2}(\mathrm{len}(B))|$, with
$\mathrm{len}$ measuring the amount of cells (of the state space tessellation)
the attractor covers, and threshold $t=0.9\bar{9}$. This effectively means
that matched attractors must have periods with ratio less than 2. b) Grouping
attractors of Fig. 2(e) so that attractors are grouped into those whose 3rd
species has population less than 0.01, or more. (c) A replication of the MCBB
(25) results for a 2nd-order Kuramoto oscillator network representing a power
grid, using the featurize-and-group continuation implementation from our
framework. Features extracted from sampled trajectories are the means of the
frequencies. (d) Same system as (c), but using the recurrence continuation and
matching attractors by their centroid distance (i.e., as in Fig. 2). The only
extra step was to post-process the results so that all attractors with basins
fractions less than 4% are aggregated (as was done in Ref. (25) and in panel
(c)). (e) Attractor basin fractions for a network of 1st order Kuramoto
oscillators; the attractors here are found and matched using the recurrences
continuation, and then grouped via a histogram of their synchronization order
parameter $R$ (19, Chap. 9). Attractors whose order parameter $R$ falls in the
same histogram bin are aggregated.

##  2 Discussion

###  2.1 Comparison with traditional continuation-based bifurcation analysis

In Table 1 and Fig. 4 we provide a careful comparison between CBA and RAFM. A
direct comparison of the two approaches is difficult, since their operational
basis, and main output, are fundamentally different. On one hand, CBA finds
and continues the curves of individual fixed points or limit cycles across the
joint state-parameter space on the basis of Newton’s method. On the other
hand, RAFM first finds attractors at all parameter values using the original
system equations, and then _matches_ appropriately similar attractors across
different parameters, giving the illusion of continuing them individually.
Additionally, the curves of stable fixed points in the joint parameter space
(Fig. 4) are only a small part of the information our framework provides.
Important provided information are the basin fractions and how they change,
which is completely absent in CBA.

Based on this comparison, we argue that the RAFM algorithm and the global
stability framework we provide is an essential tool for stability analysis of
dynamical systems. We further believe that in some application scenarios RAFM
will supersede CBA, especially given the difference in required user
expertise, required interventions, and steepness of the learning curve that
CBA has over RAFM. In other scenarios, we envision that RAFM can be used as
the default analysis method, providing the majority of information, and CBA
then becomes a more in-depth analysis of fixed points, limit cycles, and
bifurcations, if such analysis is required.

# | Traditional continuation-based bifurcation analysis (CBA) | Recurrences-based attractor find-and-match continuation (RAFM)
---|---|---
1 | provides curves of unstable fixed points / limit cycles | only finds attracting, not repelling, sets
2 | several possibilities for how to continue bifurcation curves | does not explicitly detect bifurcation points
3 | does not put limits on state space extent | needs as an input a state space box that may contain attractors
4⋄ | likely to find fixed points / limit cycles with small or even zero basins | probability to find attractor is proportional to its basins fraction
5 | detects and classifies local bifurcation points | does not compute Jacobian eigenvalues at all
6 | finds and continues fixed points and periodic orbits | finds and continues any kind of attractors, including quasiperiodic or chaotic
7 | requires a computable Jacobian of the dynamic rule | works for any dynamical system including Poincaré maps or projections
8¶ | user must manually search for multistability | different attractors are automatically detected (via random sampling)
9 | does not compute the basins of attraction or their fractions | computes the fractions and, if computationally feasible, also the full basins
10 | uses the local, linearized dynamics | uses the full nonlinear dynamics
11‡ | limited use in indicating loss of stability | more likely to indicate loss of stability as the basin fraction approaches 0
12∗ | parameter change may not affect linear stability of all fixed points | parameter change is more likely to affect global stability of all attractors
13 | no sophistication on matching fixed points | sophisticated, user-configurable matching of attractors
14† | requires expertise and constant interventions | conceptually straightforward even in advanced use-cases

Table 1: A comparison between CBA and RAFM as tools for analyzing the
stability of a dynamical system versus a parameter. The table is colored blue
or green for when to prefer CBA or RAFM respectively.

⋄Newton’s method transforms the dynamical system into a discrete time system
with different basins, making attractors with very small or zero basin sizes
have much larger ones instead. ¶In RAFM, attractors that are not being
continued from a previously found one, are found via random sampling of
initial conditions in the given state space box. The probability to find an
attractor is equal to $1-(1-f)^{n}$ with $f$ the basins fraction of the
attractor and $n$ the amount of sampled initial conditions. ‡Changing a
parameter often does not meaningfully increase the unstable eigenvalues of the
Jacobian matrix, which would indicate loss of stability (19, Chap. 12). On the
other hand, basin fractions typically decrease smoothly towards zero as an
attractor loses stability (24), although this is not guaranteed to be the case
(36), in which scenario, neither method indicates loss of stability. ∗Change
of a parameter may affect the linear stability of a single fixed point, not
all, providing flat lines in the bifurcation diagram for the unaffected fixed
points. On the contrary, loss of global stability of any attractor affects
(typically increases) the global stability of all other attractors. †Advanced
applications of traditional bifurcation analysis software require several
manual interventions during the process, and tuning of several configuration
options, many of which do not have an immediately transparent role, requiring
an expert user to make several decisions. The simplicity of our approach comes
in part because of the brute-force nature of mapping individual initial
conditions to attractors to collect the fractions, the intuitive nature of how
attractors are matched (which is also user configurable), and the lack of
necessity of interventions: after the configuration is decided, the framework
runs automatically.

Figure 4: Stability analysis of a 3-dimensional neural dynamics model (37),
plotting as information the maximum of the first variable of the dynamical
system. For the parameters used, the model undergoes saddle-node and Hopf
bifurcations and features bistability of a limit cycle and a fixed point.
Left: analysis using the BifurcationKit.jl (38) Julia package for CBA. The
analysis must happen in two steps; first the branch of a fixed point is found
and continued, and then, using a different algorithm, the branch of the limit
cycle is continued from the Hopf bifurcation. We scatter-plot special points
found and labelled by the process. The computation required $\sim$19 seconds
on an average laptop. Right: analysis of the same model using our framework.
The analysis happens in one fully automated step, after deciding the state
space box and other meta-parameters. For the analysis, we purposefully
demanded unnecessarily high accuracy, using a 9-th order ODE solver (39, 40)
with tolerances of 10-9, and requiring 1000 recurrences before claiming
convergence in the recurrence algorithm (28). The process integrated in total
101,000 initial conditions, yet required $\sim$16 seconds on the same laptop.
Attractor matching utilized a threshold: attractors whose distance in terms of
their maximum value of the first variable (i.e., same information plotted in
figure) exceeding 3.0 are not matched.

###  2.2 Comparison between attractor-finding methods

Our framework provides two radically different methods for finding attractors:
the recurrence-based and the featurize-and-group approach. Generally speaking,
the recurrence-based method should be preferred when possible, because of its
accuracy (finding the actual attractors), and the possibility for follow-up
analysis of found attractors. However, the featurize-and-group method should
be preferred when the recurrence-based method fails, because e.g., the
provided state space tessellation is ill-defined, or because computational
demands exceed what is available. In Table 2 we provide a comprehensive
comparison between the two methods.

Aspect |  Recurrences method |  Featurizing method
---|---|---
Accuracy |  Highly accurate: finds actual attractors using their unique property of their state space location. |  Less accurate, as trajectories are transformed into features, and attractors correspond to groups of features. The correspondence is not guaranteed to be unique or reversible.
Info stored |  Stores samples of points on the found attractors |  Stores a user-provided function of the group of features (by default: the centroid of the group).
Speed |  Very fast for low dimensional systems and for systems whose attractors are fixed points or periodic orbits. Becomes slow for attractors with long recurrence times (high dimensional chaotic attractors or very fine state space tessellations). |  Performance is independent of system attractors. It linearly scales with the amount of initial conditions, the transient integration time, and the total integration time. Parallelizable. In addition is the cost of the grouping process, which is huge for clustering but trivial for histograms or nearest-feature. See also the benchmark comparison in the Supplementary Information §3
Memory |  Memory allocation scales as $(1/\varepsilon)^{\Delta}$, with $\varepsilon$ the state space tessellation size and $\Delta$ the capacity dimension of the attractor, which is often much lower than the state space dimension (19). |  Memory allocation of the trajectories scales linearly with integration time and sampling rate. Additionally, specifically for clustering used as the grouping mechanism, the total memory allocated is proportional to the square of [initial conditions $\times$ parameter values] which, if one attempts to obtain accurate results, is often beyond the available memory on a typical computer.
Necessary input (guesswork) |  A state space box that may contain the attractors, and a state space tesselation that is fine enough to differentiate the location of attractors. |  A state space box that may contain the attractors; a function mapping attractors to features, such that different attractors produce different features; how much transient time to discard from time evolution; how much time to evolve and record the trajectory for, after transient.
Meta-parameters |  The parameters of the finite state machine of the recurrences algorithm discussed in Ref. (28), and the time stepping $\Delta t$. All are crucial, but all are conceptually straightforward. |  The integration time and sampling rate and all parameters of the grouping procedure. Integration parameters are straightforward, but optimal parameters related to the grouping are much harder to guess.
Troubleshooting |  Easy to troubleshoot. At any point the actual trajectories and attractors are accessible and why a failure occurs is typically easy to find out. Matching of attractors happens parameter-by-parameter, hence, individual parameter slices can be isolated and analyzed to identify where matching has failed and why (distances between attractors is also provided information). |  Difficult to troubleshoot. It does not find the actual system attractors, so the user must always reason in terms of features. When grouping using clustering (DBSCAN), the grouping process essentially operates as a black box after the features have been computed, making it harder to comprehend failures. Matching of attractors happens at the same time as grouping, making it nearly impossible to understand why an expected matching failed during the continuation process.
Failures |  Algorithm may fail if state space tessellation is not fine enough and a grid cell may contain points from different attractors. Chaotic saddles and other sticky sets generate very long transients (41) and the algorithm can interpret them as attractors. Additionally, the algorithm is sensitive to the time step $\Delta t$ and the used integrator. For limit cycles it is often the case that a non-adaptive integrator needs to be used. |  Sticky sets that are formally not attractors will be interpreted as attractors. Additionally, for ill-defined features, any group of trajectories could be interpreted as an attractor, which is not desirable in the context of this paper. Clustering via DBSCAN may fail unexpectedly, or, finding the optimal radius for the clustering may yield incorrect results.
Table 2: Comparison the two main methods of finding and continuing attractors.

###  2.3 Future extensions

The modularity of our framework allows it to be easily extended by researchers
for their specific use cases. For instance, the default information about the
attractors that is “tracked” during the RAFM continuation is their state space
location. In the future, we plan to further enhance the framework to also
(optionally) track dynamical characteristics of the attractors such as the
Lyapunov spectra or fractal dimensions (19) (computing these quantities after-
the-fact is trivial as illustrated in our example code of Listing. 1).
Utilizing this information will allow us to better label the attractors
(instead of using the integers as we currently do), better match attractors
during the continuation, or even automatically detect bifurcations or changes
in the nature of the attractor, such as transitioning from a periodic
attractor (no positive Lyapunov exponents) to a chaotic one (at least one
positive Lyapunov exponent).

### A. Featurizing methods for finding attractors

We group together two similar methods that have been recently proposed in the
literature for finding attractors. One is called Monte Carlo basin bifurcation
analysis (MCBB) (25), the other is basin stability analysis (bSTAB) (26). Both
methods work by identifying attractors as clusters of user-defined-features of
trajectories. Their first step is to integrate $N$ randomly chosen initial
conditions inside a certain box in state space. The integration is done for
some time $T$, after a transient $T_{\mathrm{tr}}$. Both $T$ needs to be
sufficiently large so that the trajectories correspond to the systems’ long-
term behavior and also $T_{\mathrm{tr}}$ needs to be large to avoid the
transient regime. Each trajectory $\vec{x}(t)$ is then transformed into a
vector of $K$ features $\vec{F}$, specified by some featurizer function
$\vec{f}$ such that $\vec{f}(\vec{x}(t))=\vec{F}$ that has to be defined by
the user. Each vector of features $\vec{F}$ describes a point in the
$K$-dimensional space of features $f_{1}\times f_{2}\times\cdots\times f_{K}$.
The key idea is that features belonging to the same attractor cluster together
in state space, so that each attractor forms a distinct cluster (a cloud of
points) in feature space. The final step in the method is to therefore cluster
the features. The clustering algorithm chosen for this is the Density-based
Spatial Clustering of Applications with Noise (DBSCAN) (27). It first
classifies two points as neighbors if their distance if smaller than a radius
$\epsilon$. Then, it clusters together points with many neighbors (equal or
more than a parameter minPts), and leave as outliers points with too few
neighbors (less than minPts). The radius $\epsilon$ is a crucial parameter for
the algorithm, and often needs fine tuning for proper clustering. The methods
by (25) and (26) use two different ways to identify a value for $\epsilon$.
Authors in (25) use a method that looks at the ordered distance of the $k$
nearest neighbors to each point in the dataset, and finds the first knee (high
derivative point) (27, 42). Authors in (26) iteratively search for the
$\epsilon$ that maximizes a criterion of clustering quality. To calculate this
criterion, they evaluate the silhouette of the each cluster, which measures
how similar each point is to the cluster it currently belongs to, compared to
the other clusters, and ranges from $-1$ (worst matching) to $+1$ (ideal
matching). This leads to one silhouette value per feature; the authors then
take the minimum value as the representative for the clustering quality for
each radius. The chosen radius is thus the value that leads to the highest
minimum silhouette. In both methods, the clusters found by DBSCAN are
considered then as attractors.

### B. Recurrences-based method for finding attractors

The inputs to this method are a dynamical system, a state space box that may
contain the attractors (although initially it may be arbitrarily large), and a
tessellation of the given box into cells. An initial condition of the system
is evolved step-by-step with time step $\Delta t$. At each step, the location
of the trajectory in state space is mapped to its cell, and that cell is
labelled as visited. If the dynamical system has attractors, and they satisfy
the Poincaré recurrence theorem (19, Chap. 9), the trajectory is guaranteed to
revisit cells it has visited before. Once a pre-decided number of recurrences
$n_{f}$ have been accumulated consecutively (i.e., enough previously-visited
cells are visited again), the method claims to have found an attractor. It
then proceeds to locate the attractor accurately, by collecting a pre-decided
number $n_{l}$ of points on the attractor (Figure 1 of Ref. (28) and panel (3)
of Fig. 1). A finite state machine formulation keeps track of coexisting
attractors, so that each attractor is unique. It also keeps track of
divergence to infinity by counting steps $n_{d}$ outside the box, ensures
algorithm termination by setting a total $n_{m}$ of max amount of $\Delta t$
iterations, and makes convergence faster by utilizing information already
encoded in the grid: if the trajectory visits consecutively a relatively small
number $n_{r}$ of cells already labelled as an attractor, convergence is
already eagerly decided. I.e., converging to an already found attractor is
much faster than finding that attractor for the first time. Hence, $\Delta
t,n_{f},n_{l},n_{d},n_{m},n_{r}$ are the meta-parameters of the algorithm and
have sensible default values that work in most cases. More information on the
method can be found in Ref. (28). Notice that the recurrence method is
different at a fundamental level from GAIO(43) and other cell mapping
techniques(44). We expand more on this in the Supplementary Information. Also
note that the method is not perfect; it may identify two attractors when only
one exists due to e.g., choosing too low convergence criteria $n_{l},n_{f}$,
or due to a commensurate period of attractor and integrator time step.
However, once again the importance of finding the “actual” attractors becomes
apparent: further analysis, by e.g., plotting the attractors, immediately
highlights such a failure and how to deal with it, and we also provide several
tips in the documentation of our method in the code implementation (45).

### C. Improvements to the recurrences method

A large drawback of the recurrences method was that it scaled poorly with the
dimension $D$ of the dynamical system. If an $\epsilon$-sized tessellation of
the state space is chosen, then memory allocated scaled as $1/\epsilon^{D}$.
We now use sparse arrays to store accessed grid locations. This changes the
memory scaling to $1/\epsilon^{\Delta}$, with $\Delta$ the capacity dimension
(19) of the attractor, which is typically much smaller than $D$ (and is only
$1$ for limit cycles).

### D. Improvements to the featurizing methods

First, we have changed the criterion used for finding the optimal radius in
the clustering method. We have found the knee method consistently more
unreliable than the iterative search. We have also found that the mean,
instead of the minimum, silhouette value as the measure of clustering quality
leads to better clustering. For instance, this lead to correct clustering in
the Lorenz86 system, whereas the minimum value criterion did not. Furthermore,
our method searches for the optimal-radius radius with a bisection method,
instead of the linear method used by authors in (26). This significantly
speeds up the code. Another simple modification we introduced is to rescale
the features in each dimension ($f_{1},f_{2},\cdots,f_{K}$) into the same
interval, for instance $[0,1]$. We noticed that the clustering method performs
poorly if the features span different ranges of values, and this simple
modification proved to be a very powerful solution. Third, we allow the
integration of all initial conditions to be done in parallel, using several
computer cores, which speeds up the solution. Lastly, grouping in our
framework can also happen based on a histogram in feature space.

### E. Code implementation

The code implementation of our framework is part of the DynamicalSystems.jl
library (29) as the Attractors.jl package (45). The code is open source code
for the Julia language, has been developed following best practices in
scientific code (46), is tested extensively, and is accompanied by a high
quality documentation. An example code snippet is shown in Listing 1.

Besides the quality of the implementation, three more features of the code are
noteworthy. First, that it is part of DynamicalSystems.jl instead of an
isolated piece of code. This integration makes the simplicity and high-
levelness of Listing 1 possible, and makes the input for the code easy to set
up. Moreover, the direct output of the code can be used with the rest of the
library to further analyze attractors in terms of e.g., Lyapunov exponents or
fractal dimensions. Indeed, in the provided code example Listing 1 we compute
the Lyapunov spectra of all found attractors, across all parameter values, in
only two additional lines of code. Second, utilizing the Julia language’s
multiple dispatch system (47), the code is extendable. It establishes one
interface for how to map initial conditions to attractors, and one for how to
group features, both of which can be extended, and yet readily be usable by
the rest of the library such as the continuation methods. Third, a lot of
attention has been put into user experience, by establishing a short learning
curve via a minimal user interface, by carefully considering how to provide
the output in an intuitive format, as well as providing easy-to-use plotting
functions that utilize the code output. More overview and information on the
code or its design can be found in its online documentation or source code
(45).

### F. Code for this article

The code we used to create the figures of this article is fully reproducible
and available online (48).

Listing 1: Julia code snippet showcasing the usage of the DynamicalSystems.jl
implementation of our framework. The code produces panel (a) of Fig. 2. The
main output of the code are two vectors, containing the basins fractions and
attractors at each parameter value respectively. The fractions and attractors
are formulated as dictionaries, mapping attractor labels (the integers) to
basin fractions and sets of points on the attractor, respectively. At its end,
the code snippet computes the Lyapunov spectra of all found attractors, by
using the first point on each attractor as initial condition for the
computation of the Lyapunov spectrum.

⬇

using DynamicalSystems # our framework implementation

using OrdinaryDiffEq # high-accuracy ODE solvers

# create Lorenz84 within DynamicalSystems.jl

function lorenz84_rule(u, p, t)

F, G, a, b = p

x, y, z = u

dx = -y^2 -z^2 -a*x + a*F

dy = x*y - y - b*x*z + G

dz = b*x*y + x*z - z

return SVector(dx, dy, dz)

end

u0 = ones(3) # init. state

p0 = [6.886, 1.347, 0.255, 4.0] # init. parameters

# ODE solver:

diffeq = (alg = Vern9(), reltol = 1e-9, abstol = 1e-9)

# Main object of the library: a ‘DynamicalSystem‘

ds = CoupledODEs(lorenz84_rule, u0, p0; diffeq)

# Provide state space box tessellation to search in

xg = yg = zg = range(-3, 3; length = 600)

grid = (xg, yg, zg)

# initialize recurrences-based algorithm

# and choose its metaparameters

mapper = AttractorsViaRecurrences(ds, grid;

mx_chk_fnd_att = 1000, mx_chk_loc_att = 2000,

mx_chk_lost = 100, mx_chk_safety = 1e8,

Dt = 0.05, force_non_adaptive = true,

)

# find and continue attractors across a given

# parameter range for the ‘pidx‘-th parameter

prange = range(1.34, 1.37; length = 101)

pidx = 2 # index of parameter

sampler = statespace_sampler(grid)[1]

rsc = RAFM(mapper)

# main output:

fractions_curves, attractors_info = continuation(

rsc, prange, pidx, sampler

)

# Estimate Lyapunov spectra for all attractors

# by looping over the parameter range

lyapunovs_curves = map(eachindex(prange)) do index

set_parameter!(ds, pidx, prange[index])

attractor_dict = attractors_info[index]

exponents = Dict(

k => lyapunovspectrum(ds, 10000; u0 = A[1])

for (k, A) in attractor_dict

)

end

The authors would like to thank Ulrike Feudel, Peter Ashwin, Ulrich Parlitz,
and Harry Dankowicz for helpful discussions. G.D. was supported by the Royal
Society International Newton Fellowship. K.L.R. was supported by the German
Academic Exchange Service (DAAD). A.W. was supported by the Spanish State
Research Agency (AEI) and the European Regional Development Fund.

## References

  * (1) U Feudel, AN Pisarchik, K Showalter, Multistability and tipping: From mathematics and physics to climate and brain - Minireview and preface to the focus issue.  Chaos 28, 33501 (2018).
  * (2) F Hellmann, et al., Network-induced multistability through lossy coupling and exotic solitary states.  Nature Communications 11, 592 (2020).
  * (3) H Kim, SH Lee, J Davidsen, SW Son, Multistability and variations in basin of attraction in power-grid systems.  New Journal of Physics 20, 113006 (2018).
  * (4) L Halekotte, A Vanselow, U Feudel, Transient chaos enforces uncertainty in the British power grid.  Journal of Physics: Complexity 2, 035015 (2021).
  * (5) J Marotzke, P Welander, J Willebrand, Instability and multiple steady states in a meridional-plane model of the thermohaline circulation.  Tellus A 40, 162–172 (2016).
  * (6) TM Lenton, Environmental Tipping Points.  Annual Review of Environment and Resources 38, 1–29 (2013).
  * (7) M Hirota, M Holmgren, EHV Nes, M Scheffer, Global Resilience of Tropical Forest and Savanna to Critical Transitions.  Science 334, 232–235 (2011).
  * (8) V Dakos, et al., Ecosystem tipping points in an evolving world.  Nature Ecology & Evolution 3, 355–362 (2019).
  * (9) JL Schwartz, N Grimault, JM Hupé, BCJ Moore, D Pressnitzer, Multistability in perception: binding sensory modalities, an overview.  Philosophical Transactions of the Royal Society B: Biological Sciences 367, 896–905 (2012).
  * (10) JAS Kelso, Multistability and metastability: Understanding dynamic coordination in the brain.  Philosophical Transactions of the Royal Society B: Biological Sciences 367, 906–918 (2012).
  * (11) A Kleinschmidt, P Sterzer, G Rees, Variability of perceptual multistability: from brain state to individual trait.  Philosophical Transactions of the Royal Society B: Biological Sciences 367, 988–1000 (2012).
  * (12) R Zhu, JMd Rio-Salgado, J Garcia-Ojalvo, MB Elowitz, Synthetic multistability in mammalian cells.  Science 375, eabg9765 (2022).
  * (13) T Khazaei, et al., Metabolic multistability and hysteresis in a model aerobe-anaerobe microbiome community.  Science Advances 6, eaba0353 (2020).
  * (14) C Geiß, E Salas, J Guevara-Coto, A Régnier-Vigouroux, RA Mora-Rodríguez, Multistability in Macrophage Activation Pathways and Metabolic Implications.  Cells 11, 404 (2022).
  * (15) K Padiyar, Power System Dynamics: Stability and Control.  (Wiley), (1999).
  * (16) J Lohmann, PD Ditlevsen, Risk of tipping the overturning circulation due to increasing rates of ice melt.  Proc. Natl. Acad. Sci. U. S. A. 118, e2017989118 (2021).
  * (17) P Ashwin, S Wieczorek, R Vitolo, P Cox, Tipping points in open systems: Bifurcation, noise-induced and rate-dependent examples in the climate system.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370, 1166–1184 (2012).
  * (18) L Halekotte, U Feudel, Minimal fatal shocks in multistable complex networks.  Scientific Reports 10, 11783 (2020).
  * (19) G Datseris, U Parlitz, Nonlinear dynamics, Undergraduate Lecture Notes in Physics.  (Springer Nature, Cham, Switzerland), 1 edition, (2022).
  * (20) EJ Doedel, Auto: A program for the automatic bifurcation analysis of autonomous systems.  Congr. Numer 30, 25–93 (1981).
  * (21) A Dhooge, W Govaerts, YA Kuznetsov, Matcont: a matlab package for numerical bifurcation analysis of odes.  ACM Transactions on Mathematical Software (TOMS) 29, 141–164 (2003).
  * (22) H Dankowicz, Recipes for continuation, Computational science & engineering.  (Society for Industrial and Applied Mathematics, Philadelphia, Pa.), (2013).
  * (23) R Veltz, BifurcationKit.jl (2020).
  * (24) PJ Menck, J Heitzig, N Marwan, J Kurths, How basin stability complements the linear-stability paradigm.  Nature Physics 9, 89–92 (2013).
  * (25) M Gelbrecht, J Kurths, F Hellmann, Monte carlo basin bifurcation analysis.  New J. Phys. 22, 033032 (2020).
  * (26) M Stender, N Hoffmann, bSTAB: an open-source software for computing the basin stability of multi-stable dynamical systems.  Nonlinear Dynamics (2021).
  * (27) M Ester, HP Kriegel, J Sander, X Xiaowei, A density-based algorithm for discovering clusters in large spatial databases with noise.  Proceedings of the Second International Conference on Knowledge Discovery and Data Mining (1996).
  * (28) G Datseris, A Wagemakers, Effortless estimation of basins of attraction.  Chaos 32 (2022).
  * (29) G Datseris, Dynamicalsystems.jl: A julia software library for chaos and nonlinear dynamics.  Journal of Open Source Software 3, 598 (2018).
  * (30) EN Lorenz, Irregularity: A fundamental property of the atmosphere.  Tellus A: Dynamic Meteorology and Oceanography 36, 98–110 (1984).
  * (31) M Gelbrecht, V Lucarini, N Boers, J Kurths, Analysis of a bistable climate toy model with physics-based machine learning methods.  European Physical Journal: Special Topics 123 (2021).
  * (32) S Huang, Multistability and multicellularity: cell fates as high-dimensional attractors of gene regulatory networks in Computational Systems Biology.  (Elsevier), pp. 293–326 (2006).
  * (33) J Moehlis, H Faisst, B Eckhardt, A low-dimensional model for turbulent shear flows.  New Journal of Physics 6, 56 (2004).
  * (34) J Huisman, FJ Weissing, Fundamental Unpredictability in Multispecies Competition.  The American Naturalist 157, 488–494 (2001).
  * (35) F Hausdorff, Grundzuge Der Mengenlehre, AMS Chelsea Publishing.  (American Mathematical Society, Providence, RI), (1949).
  * (36) P Schultz, PJ Menck, J Heitzig, J Kurths, Potentials and limits to basin stability estimation.  New Journal of Physics 19 (2017).
  * (37) JM Cortes, et al., Short-term synaptic plasticity in the deterministic tsodyks–markram model leads to unpredictable network dynamics.  Proceedings of the National Academy of Sciences 110, 16610–16615 (2013).
  * (38) R Veltz, BifurcationKit.jl (2020).
  * (39) JH Verner, Numerically optimal runge–kutta pairs with interpolants.  Numerical Algorithms 53, 383–396 (2010).
  * (40) C Rackauckas, Q Nie, Differentialequations.jl – a performant and feature-rich ecosystem for solving differential equations in julia.  The Journal of Open Research Software 5 (2017) Exported from https://app.dimensions.ai on 2019/05/05.
  * (41) YC Lai, T Tél, Transient Chaos: Complex Dynamics on Finite-TIme Scales, Applied Mathematical Sciences.  (Springer New York, NY), (2009).
  * (42) M Hahsler, M Piekenbrock, D Doran, dbscan: Fast density-based clustering with r.  Journal of Statistical Software 91, 1–30 (2019).
  * (43) R Gerlach, A Ziessler, B Eckhardt, M Dellnitz, A set-oriented path following method for the approximation of parameter dependent attractors.  SIAM Journal on Applied Dynamical Systems 19, 705–723 (2020).
  * (44) JQ Sun, FR Xiong, O Schütze, C Hernández, Cell mapping methods.  (Springer), (2018).
  * (45) G Datseris, KL Rossi, A Wagemakers, Juliadynamics/attractors.jl: v1.2.5 (2023).
  * (46) G Datseris, Good scientific code workshop (2023).
  * (47) J Bezanson, A Edelman, S Karpinski, VB Shah, Julia: A fresh approach to numerical computing.  SIAM review 59, 65–98 (2017).
  * (48) G Datseris, Datseris/frameworkglobalstability: v0.1 (2023).


# A generalized expansion method for computing Laplace–Beltrami eigenfunctions
on manifolds

Jackson C. Turner  Department of Applied Physics and Applied Mathematics,
Columbia University, New York City 10027<EMAIL_ADDRESS>,  Elena
Cherkaev  Department of Mathematics, University of Utah, Salt Lake City 84112
<EMAIL_ADDRESS>and  Dong Wang  School of Science and Engineering, The
Chinese University of Hong Kong, Shenzhen, Guangdong, 518172, China
Shenzhen International Center for Industrial and Applied Mathematics, Shenzhen
Research Institute of Big Data, Guangdong, 518172, China<EMAIL_ADDRESS>
###### Abstract.

Eigendecomposition of the Laplace–Beltrami operator is instrumental for a
variety of applications from physics to data science. We develop a numerical
method of computation of the eigenvalues and eigenfunctions of the
Laplace–Beltrami operator on a smooth bounded domain based on the relaxation
to the Schrödinger operator with finite potential on a Riemannian manifold and
projection in a special basis. We prove spectral exactness of the method and
provide examples of calculated results and applications, particularly, in
quantum billiards on manifolds.

###### Key words and phrases:

Laplace operator, Laplace–Beltrami operator, fictitious domain methods,
expansion method, quantum billiards

###### 2020 Mathematics Subject Classification:

35J05, 65N85, 47A70, 65N25,

##  1\. Introduction

The Laplace–Beltrami operator plays an important role in the differential
equations that describe many physical systems. These include, for example,
vibrating membranes, fluid flow, heat flow, and solutions to the Schrödinger
equation. Another example is that of spectral partitions—collections of $k$
pairwise disjoint open subsets such that the sum of their first
Laplace–Beltrami eigenvalues is minimal [13, 14, 38, 39]. This has a wide
class of applications including data classification [32], interacting agents
[18, 16, 19], and so on. In all the above applications, the fundamental
question is how to efficiently compute the eigenvalues of the Laplace–Beltrami
operator in an arbitrary domain with a proper boundary condition. Also, the
Laplace–Beltrami operator is crucial to understanding systems described by
nonlinear Schrödinger equations, such as the propagation of Langmuir waves in
an ionized plasma [24, 7], the single-particle ground-state wavefunction in a
Bose–Einstein condensate [7], the slowly-varying envelope of light waves in
Kerr media [20], and water surface wave packets [41].

The Laplace–Beltrami operator of a scalar function $f$ on a Riemannian
manifold $(\mathcal{M},g)$ is defined as the surface divergence of the vector
field gradient of $f$,

(1.1) $\Delta_{g}f=\nabla\cdot\nabla f.$

The divergence of a vector field $X$ with metric $g$ is (in Einstein notation)

(1.2) $\nabla\cdot
X=\frac{1}{\sqrt{|g|}}\partial_{i}\mathopen{}\mathclose{{}\left(\sqrt{|g|}X^{i}}\right),$

and the gradient of a scalar function $f$ is

(1.3) $(\operatorname{grad}f)^{i}=\partial^{i}f=g^{ij}\partial_{j}f.$

From above, we obtain the Laplace–Beltrami operator acting on functions over
$(\mathcal{M},g)$,

(1.4)
$\Delta_{g}f=\frac{1}{\sqrt{|g|}}\partial_{i}\mathopen{}\mathclose{{}\left(\sqrt{|g|}g^{ij}\partial_{j}f}\right).$

In general, the Helmholtz equation (i.e. Laplace–Beltrami eigenvalue problem)
with Dirichlet boundary conditions on $\Omega\subset(\mathcal{M},g)$ is

(1.5) $\begin{cases}-\Delta_{g}u(x)=\lambda u(x),&x\in\Omega\\\
u(x)=0,&x\in\partial\Omega.\end{cases}$

In this paper, we develop a numerical method to find the eigenvalues and
eigenfunctions of the Laplace–Beltrami operator with Dirichlet and periodic
boundary conditions for arbitrary domains on various surfaces. The idea is
highly motivated by the Schrödinger operator and is based off the method given
in [27]. By using the Schrödinger operator relaxation, we relax the eigenvalue
problem on an arbitrary domain into the eigenvalue problem for a Schrödinger
operator on a regular domain which is convenient for the numerical
discretization.

In [34], comparable methods on manifolds using linear and cubic FEM operators
and discrete geometric Laplacians are explored, and [17] provides a method for
hyperbolic domains. There is extensive literature on the Laplacian for planar
regions [26, 29, 12, 8, 1]. In methods for solving nonlinear Schrödinger
equations, finite difference discretizations of the Laplace operator are often
used [5, 6, 16]. It is likely many of these methods above can be extended to
the Laplace–Beltrami operator on manifolds. The method we present in this
paper has some immediate advantages over the finite difference method—since
the boundary of domains are characterized by a potential function (see Theorem
7), no creation of a complicated mesh is needed, allowing for more generic
domains and producing smooth solutions. Also, the method has promise to be
quite robust in discretizing the operator on domains with corners (as in Table
I), especially in applications when computation of many eigenvalues is
required, whereas the finite difference method is notoriously inefficient on
such domains.

The rest of the paper is organized as follows. In Section 2, we recall the
Schrödinger operator and introduce the generalized expansion method. We
discuss and prove the convergence and accuracy of the relaxation and
approximation in Section 3 and show extensive numerical experiments in Section
4. We investigate applications to spherical domains, periodic domains, and
billiard problems in Section 5 and draw some conclusion in Section 6.

##  2\. Generalized Expansion Method

The time-independent Schrödinger equation on a Riemannian manifold
$\mathcal{M}$ with metric $g$, potential ${V}(x)$, and energy levels $E_{n}$
is

(2.1)
$\hat{H}\psi_{n}(x)=\mathopen{}\mathclose{{}\left[-\Delta_{g}+{V}(x)}\right]\psi_{n}(x)=E_{n}\psi_{n}(x),\\\
$

where $\Delta_{g}$ is the Laplace–Beltrami operator on $(\mathcal{M},g)$ as in
(1.4). The Schrödinger equation is an eigenvalue problem for the Schrödinger
operator $\hat{H}=-\Delta_{g}+V(x)$. When

(2.2)
${V}(x)=\mathopen{}\mathclose{{}\left\\{\begin{array}[]{ll}0&x\in\Omega\\\
\infty&x\not\in\Omega,\end{array}}\right.$

the eigenvalue problem for $\hat{H}$ is equivalent to (1.5). Eigenfunctions
are normalized by setting

(2.3) $\int_{\Omega}|\psi_{n}(x)|^{2}dx=1,$

where $|\psi_{n}(x)|^{2}dx$ is a probability density.

In [27], a method is given to solve (1.5) with $g=I_{2}$ on any bounded smooth
$\Omega\subset\mathbb{R}^{2}$ by embedding it in a rectangle, as in Figure 1.
In order to evaluate the Laplace–Beltrami eigenvalues for $\Omega$ on a 2-D
surface, we generalize this method when considering $\Omega$ as a smooth
subset of a bounded manifold $S=(\mathcal{M},g)$ using a complete set of
orthonormal eigenfunctions
$\mathcal{F}^{S}_{\infty}=\\{\phi_{n}\\}_{n=1}^{\infty}$ on $S$ with
corresponding eigenvalues $\\{\lambda_{n}(S)\\}_{n=1}^{\infty}$ (with
$\lambda_{1}<\lambda_{2}\leq\lambda_{3}\leq\dots$). Here, we assume all
$\phi_{n}\in\mathcal{F}^{S}_{\infty}$ have Dirichlet boundary conditions on
$\partial S$, but for cases when $|\partial\Omega\cap\partial S|>0$, one may
choose to use other boundary conditions to obtain solutions of (1.5) with
$\mathopen{}\mathclose{{}\left.u}\right|_{x\in\partial\Omega\cap\partial
S}\not\equiv 0$, as in the periodic examples in Section 5.2.

$S$$\Omega$ Figure 1. With the expansion method, a bounded domain $\Omega$ is
embedded onto a rectangle $S$ with euclidean geometry and $\hat{H}$ is
projected onto $\mathcal{F}^{S}_{N}$.

In this method, we use $\mathcal{F}^{S}_{N}=\\{\phi_{n}\\}_{n=1}^{N}$ as a
basis on which we expand the operator $\hat{H}$ and seek solutions of (1.5),
or the equivalent equation involving the eigenvalue problem of the Schrödinger
operator $\hat{H}$,

(2.4)
$\hat{H}\psi(x)=\mathopen{}\mathclose{{}\left[-\Delta_{g}+\tilde{V}(x)}\right]\psi(x)=\lambda\psi(x),$

with $\tilde{V}(x)$ defined as

(2.5)
$\tilde{V}(x)=\mathopen{}\mathclose{{}\left\\{\begin{array}[]{ll}0&x\in\Omega\\\
\infty&x\in S\setminus\Omega.\end{array}}\right.$

We approximate $\tilde{V}(x)$ as

(2.6)
$V(x)=V_{0}\chi_{S\setminus\Omega}(x)=\mathopen{}\mathclose{{}\left\\{\begin{array}[]{ll}0&x\in\Omega\\\
V_{0}\gg 1&x\in S\setminus\Omega.\end{array}}\right.$

This allows us to discretize the operator $\hat{H}\approx H_{N}$,

(2.7)
$H_{N_{nm}}=\langle\phi_{n},\hat{H}\phi_{m}\rangle=\lambda_{n}^{S}\delta_{nm}+\int_{S}V(x)\phi^{*}_{n}\phi_{m}dx,$

where we truncate $n,m\leq N$. The eigenvalues and eigenvectors of $H_{N}$
approximate those of $\hat{H}$.

##  3\. Convergence analysis

In this section, we provide a rigorous proof on the convergence of the
proposed method in the sense of $V_{0}\rightarrow\infty$ and
$N\rightarrow\infty$. To keep self-consistency of the paper, we first recall
some definitions and preliminary results from [25, 36, 22].

######  Definition 1.

Suppose $A_{n}$ and $A$ are self-adjoint operators. We say that $A_{n}$
converges to $A$ in the strong resolvent sense, if

(3.1) $\|(R_{A_{n}}(z)-R_{A}(z))\phi\|\to
0,\quad\forall\phi\in\mathfrak{D}(A_{n})$

for some
$z\in\Gamma=\mathbb{C}\setminus\Sigma,\Sigma=\sigma(A)\cup\bigcup_{n}\sigma\mathopen{}\mathclose{{}\left(A_{n}}\right)$
where the function $R_{A}(z)$ is the resolvent of $A$.

######  Definition 2.

A subset $\mathfrak{D}_{0}\subseteq\mathfrak{D}(A)$ is a core of $A$ when
$\\{(x,Ax):x\in\mathfrak{D}_{0}\\}$ is dense in
$\\{(x,Ax):x\in\mathfrak{D}(A)\\}$.

######  Lemma 3.

(6.36 of [36]): Let $A_{n},A$ be self-adjoint operators. Then $A_{n}$
converges to $A$ in the strong resolvent sense if there is a core
$\mathfrak{D}_{0}$ of $A$ such that for any $\psi\in\mathfrak{D}_{0}$ we have
$P_{n}\psi\in\mathfrak{D}\mathopen{}\mathclose{{}\left(A_{n}}\right)$ for $n$
sufficiently large and $A_{n}P_{n}\psi\rightarrow A\psi$.

######  Theorem 4.

(6.38 of [36]). Let $A_{n}$ and $A$ be self-adjoint operators. If $A_{n}$
converges to $A$ in the strong resolvent sense, we have
$\sigma(A)\subseteq\lim_{n\rightarrow\infty}\sigma\mathopen{}\mathclose{{}\left(A_{n}}\right)$.

######  Theorem 5.

(2.2.3 of [25]). Let $L_{n}$ be a sequence of uniformly elliptic operators
defined on an open set $D$ by

(3.2) $L_{n}u:=-\sum_{i,j=1}^{N}\frac{\partial}{\partial
x_{i}}\mathopen{}\mathclose{{}\left(a_{ij}^{n}(x)\frac{\partial u}{\partial
x_{j}}}\right)+a_{0}^{n}(x)u.$

We assume that, for fixed $i,j$, the sequence $a_{i,j}^{n}$ is bounded in
$L^{\infty}$ and converge almost everywhere to a function $a_{i,j}$; we also
assume that the sequence $a_{0}^{n}$ is bounded in $L^{\infty}$ and converges
weakly-* in $L^{\infty}$ to a function $a_{0}$. Let $L$ be the (elliptic)
operator defined on $D$ as in (3.2) by the functions $a_{i,j}$ and $a_{0}$.
Then each eigenvalue of $L_{n}$ converges to the corresponding eigenvalue of
$L$.

######  Theorem 6.

(9.29 of [22]) Let $u\in W^{2,p}(S)\cap C^{0}(\bar{S})$ satisfy $Lu=f$ in
$S,u=\varphi$ on $\partial S$ where $f\in L^{p}(S),\varphi\in
C^{\beta}(\bar{S})$ for some $\beta>0$, and suppose that $\partial S$
satisfies a uniform exterior cone condition. Then $u\in C^{\alpha}(\bar{S})$
for some $\alpha>0$.

Now, we prove spectral exactness of the expansion method in $V_{0}$ and $N$ in
Theorems 7 and 9. We give an example of calculating solutions and the rate of
convergence in $V_{0}$ of eigenvalues for the relaxed problem on an interval
in Example 8. We provide intuition for efficient implementation of the
expansion method in Remark 10.

######  Theorem 7.

The eigenvalues of the Schrödinger operator

(3.3) $\hat{H}(V_{0})=-\Delta_{g}+V(x),\qquad
V(x)=\mathopen{}\mathclose{{}\left\\{\begin{array}[]{ll}0,&x\in\Omega\\\
V_{0},&x\not\in\Omega\end{array}}\right.$

acting on a bounded Riemannan manifold $(\mathcal{M},g)$, with $\Omega$ smooth
in $(\mathcal{M},g)$ converge monotonically to the eigenvalues of
$-\Delta_{g}$ with Dirichlet boundary conditions on $\Omega$ as
$V_{0}\rightarrow\infty$.

###### Proof.

By considering the volume form on the manifold, we have the inner product:

(3.4) $\langle
f_{1},f_{2}\rangle_{g}=\int_{\mathcal{M}}\overline{f_{1}(x)}f_{2}(x)(\operatorname{det}g)^{1/2}dx_{1}\cdots
dx_{n}$

Now, from the Rayleigh quotient of an elliptic linear operator $\mathcal{L}$
on a Riemannian manifold $(\mathcal{M},g)$,

(3.5) $R(\mathcal{L},u)=\frac{\langle u,\mathcal{L}u\rangle_{g}}{\langle
u,u\rangle_{g}},$

we have

(3.6) $R(\hat{H}(V_{0}),u)=\frac{\langle u,-\Delta_{g}u\rangle_{g}}{\langle
u,u\rangle_{g}}+\frac{\langle u,Vu\rangle_{g}}{\langle u,u\rangle_{g}}.$

By the Courant-Fischer formula,

(3.7)
$\lambda_{k}(\hat{H}(V_{0}))=\inf_{\mathcal{Q}\in\mathcal{Q}_{k}}\sup_{u\in\mathcal{Q}}\mathopen{}\mathclose{{}\left(R(\hat{H},u)}\right)$

with $\mathcal{Q}_{k}$ being the family of subspaces of
$H_{0}^{1}(\mathcal{M})$ of dimension $k$, we obtain the following inequality,

(3.8)
$\lambda_{k}(\hat{H}(V_{0}^{*}))\geq\lambda_{k}(\hat{H}(V_{0})),\quad\text{for
}V_{0}^{*}>V_{0},$

giving us monotonicity. We also have for $u\in H^{1}_{0}(\mathcal{M})$ and as
$V_{0}\rightarrow\infty$, $R(\hat{H}(V_{0}),u)<\infty$ if and only if
$\operatorname{supp}(u)\subseteq\Omega$ almost everywhere. Hence, we have

(3.9)
$\lim_{V_{0}\rightarrow\infty}\lambda_{k}(\hat{H}(V_{0}))=\inf_{\mathcal{Q}^{*}\in\mathcal{Q}^{*}_{k}}\sup_{u\in\mathcal{Q}^{*}}\mathopen{}\mathclose{{}\left(R(\hat{H},u)}\right)$

with $\mathcal{Q}^{*}_{k}$ being the family of subspaces of
$H_{0}^{1}(\Omega)$ of dimension $k$. This is precisely the Courant-Fischer
definition of the eigenvalues of the Laplace–Beltrami operator with Dirichlet
boundary conditions on $\Omega$. ∎

######  Example 8.

Consider the regions $S=(0,2)$ and $\Omega=(0,1)$. The eigenvalues of the
Helmholtz equation on $\Omega$,

(3.10) $-v_{k}^{\prime\prime}=\mu_{k}v_{k},\qquad v_{k}(0)=v_{k}(1)=0,$

can be approximated by the eigenvalues of the Schrödinger operator
$\hat{H}_{V_{0}}=-\partial_{x}^{2}+V(x;V_{0})$ on $S$ with
$V(x;V_{0})=V_{0}\chi_{S\setminus\Omega}(x)$ and large $V_{0}\gg 1$, with

(3.11) $\hat{H}_{V_{0}}u_{k}=\lambda_{k}(V_{0})u_{k},\qquad u(0)=u(2)=0,$

and rate of convergence

(3.12) $|\lambda_{k}(V_{0})-\mu_{k}|\sim\frac{1}{\sqrt{V_{0}}}$

up to a constant, as $V_{0}\uparrow\infty.$

###### Proof.

We have (3.10) has eigenvalues $\mu_{k}=\pi^{2}k^{2}$ for
$k\in\mathbb{N}\setminus\\{0\\}$. Solutions to (3.11) are of the form:

(3.13) $u_{k}=\begin{cases}\sin\sqrt{\lambda_{k}}x,&x\in(0,1]\\\
\frac{\sin\sqrt{\lambda_{k}}}{\sinh\sqrt{V_{0}-\lambda_{k}}}\sinh(\sqrt{V_{0}-\lambda_{k}}(2-x)),&x\in(1,2).\end{cases}$

By setting $u_{k}\in C^{1}(0,2)$, we arrive at

(3.14)
$\displaystyle\sqrt{\lambda_{k}}\cot\sqrt{\lambda_{k}}=-\sqrt{V_{0}-\lambda_{k}}\coth\sqrt{V_{0}-\lambda_{k}}\implies\frac{\sqrt{\lambda_{k}}}{\sin\sqrt{\lambda}_{k}}\sim\sqrt{V_{0}},$

as $V_{0}\uparrow\infty$. We then have by Taylor expansion of $\sin(\cdot)$
about $\sqrt{\mu_{k}}$,

(3.15)
$\displaystyle\sin\sqrt{\lambda_{k}}=\sin(\sqrt{\lambda_{k}}-\sqrt{\mu_{k}}+\sqrt{\mu_{k}})\approx\pm(\sqrt{\lambda_{k}}-\sqrt{\mu_{k}})\sim\frac{\sqrt{\lambda_{k}}}{\sqrt{V_{0}}}\sim\frac{1}{\sqrt{V_{0}}}$

since $\sin\sqrt{\mu_{k}}=0$. The last relation above arises from
$\sqrt{\lambda_{k}}\sim\sqrt{\mu_{k}}$, a constant. Hence altogether,

(3.16)
$|\lambda_{k}-\mu_{k}|=|\sqrt{\lambda_{k}}-\sqrt{\mu_{k}}|\cdot|\sqrt{\lambda_{k}}+\sqrt{\mu_{k}}|\sim\frac{1}{\sqrt{V_{0}}}.$

∎

######  Theorem 9.

Given a complete orthonormal basis $\mathcal{F}^{S}_{\infty}\subset
H^{1}_{0}(S)$ of Laplace–Beltrami eigenfunctions on a bounded smooth domain
$S=(\mathcal{M},g)$, the Dirichlet eigenvalues of the $N$-dimensional operator
$H_{N}$ where $H_{N_{ij}}=\langle\phi_{i},\hat{H}\phi_{j}\rangle$ for $i,j\leq
N$ converge to those of $\hat{H}$ as $N\rightarrow\infty$ where

(3.17) $\hat{H}=-\Delta_{g}+V(x),\qquad V(x)\in L^{\infty}(S,\mathbb{R}).$

The results in this theorem hold for Neumann and periodic boundary conditions
as well, using the appropriate basis and Sobolev space.

###### Proof.

We first note

(3.18)
$H_{N}:\mathcal{F}^{S}_{N}\to\mathcal{F}^{S}_{N}\qquad\hat{H}:H^{1}_{0}(S)\to
L^{2}(S),$

where $\mathcal{F}^{S}_{N}=\operatorname{span}\\{\phi_{j}\\}_{j=1}^{N}$. We
also have, by the definition of $H_{N}$, with
$u=\sum_{j=1}^{\infty}c_{j}\phi_{j}$:

(3.19) $\displaystyle H_{N}u$
$\displaystyle:=\sum_{i=1}^{N}\mathopen{}\mathclose{{}\left(\sum_{k=1}^{N}\langle\phi_{i},\hat{H}\phi_{k}\rangle
c_{k}}\right)\phi_{i}=\sum_{j=1}^{N}\langle\phi_{j},\hat{H}u\rangle\phi_{j}=P_{N}\hat{H}u$
$\displaystyle P_{N}v$
$\displaystyle:=\sum_{j=1}^{N}\langle\phi_{j},v\rangle\phi_{j}$

Furthermore, by convention we may extend the domain of $H_{N}$ to
$H^{1}_{0}(S)$ by the following extension, which we will now use in the proof:

(3.20) $H_{N}u:=P_{N}\hat{H}P_{N}u.$

Without loss of generality, because $V(x)$ is bounded we can assume further
that $V(x)\in[0,\infty)$ since the resulting spectra are merely shifted by a
constant $\alpha$ when adding $\alpha$ to $V(x)$, hence we have
$\operatorname{ker}(\hat{H}^{*})=\emptyset$ by injectivity and self-
adjointness, giving us $\overline{\operatorname{ran}(\hat{H})}=L^{2}(S)$. We
consider the set $\mathfrak{D}_{0}=\mathcal{F}^{S}_{\infty}\cap H^{1}_{0}(S)$,
and we have $\overline{\mathfrak{D}_{0}}=L^{2}(S)=\overline{H^{1}_{0}(S)}$,
and
$\overline{H\mathfrak{D}_{0}}=\overline{\operatorname{span}\\{H\phi_{j}\\}_{j=1}^{\infty}}\cap\overline{L^{2}(S)}=L^{2}(S)$,
by $V\in L^{\infty}(S)$. Hence, the graph
$\\{(x,\hat{H}x):x\in\mathfrak{D}_{0}\\}$ is dense in
$\\{(x,\hat{H}x):x\in\mathfrak{D}(\hat{H})\\}$, therefore $\mathfrak{D}_{0}$
is a core of $\hat{H}$. Furthermore, we have for all $\psi\in\mathfrak{D}_{0}$
that $P_{N}\psi\in\operatorname{span}\\{\phi_{j}\\}_{j=1}^{N}$ and
$H_{N}\psi=P_{N}\hat{H}P_{N}\psi\to\hat{H}P_{N}\psi$, so by Lemma 3 we have
strong convergence in the resolvent sense, and hence the conditions for
Theorem 4 are met, and we have

(3.21) $\lim_{N\to\infty}\sigma(H_{N})\supseteq\sigma(\hat{H}).$

Now, it is well-known these operators have purely point spectra. We consider
$V_{N}=P_{N}V$ and consider $\lambda_{N}\to\lambda_{*}$, some converging
sequence of eigenvalues of $H_{N}$ with corresponding eigenvectors $u_{N}$. We
have $P_{N}u_{N}=u_{N}$ and $P_{N}\Delta_{g}u_{N}=\Delta_{g}u_{N}$, giving us

$\displaystyle H_{N}u_{N}=(-\Delta_{g}+V_{N})u_{N}=\lambda_{N}u_{N}$
$\displaystyle\implies\lambda_{N}\in\sigma(-\Delta_{g}+V_{N}),$
$\displaystyle\sigma(-\Delta_{g}+V_{N})\longrightarrow\sigma(\hat{H})$
$\displaystyle\implies\lambda_{*}\in\sigma(\hat{H})$

by Theorem 5, giving us

(3.22) $\lim_{N\to\infty}\sigma(H_{N})\subseteq\sigma(\hat{H}).$

Altogether, we have the desired result,

(3.23) $\lim_{N\to\infty}\sigma(H_{N})=\sigma(\hat{H}).$

Similar proofs can be made for the Neumann and periodic boundary cases. ∎

######  Remark 10.

For fixed $\Omega,\ V_{0},$ and $N$, an efficient implementation of this
method is to seek an integrable domain $S\supset\Omega$ to minimize the
following value involving the $L^{2}$-induced norm of the difference of the
operators acting on the finite-dimensional space $\mathcal{F}^{S}_{N}$:

(3.24)
$\tau_{\Omega}(S):=V_{0}^{-1}\cdot\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left.[\hat{H}-H_{N}(S)]}\right|_{\mathcal{F}^{S}_{N}}}\right\|_{\mathcal{B}(L^{2}(\Omega),L^{2}(\Omega))}=V_{0}^{-1}\cdot\mathopen{}\mathclose{{}\left\|V-\sum_{j=1}^{N}\langle\phi_{j},V\rangle\phi_{j}}\right\|_{L^{\infty}(\Omega)}$

This procedure is equivalent to fitting a domain $\Omega$ properly into a
solvable set $S$ so that the potential $V$ that characterizes $\Omega$ is
well-approximated by $P_{N}V$.

##  4\. Numerical Accuracy

Note that contrary to Theorem 7, in practice, the method loses accuracy if
$V_{0}$ is chosen to be too large, due to floating point round-off errors. We
provide a numerical example here. One can identify the shape of a triangle
given the spectrum of the solution to the Helmholtz equation [23], and a
formula can be derived for the Laplacian eigenvalues $\lambda_{n}$ of an
equilateral triangle with Dirichlet boundary conditions [11], for positive
integers $p$ and $q$,

(4.1)
$\lambda_{n}\equiv\lambda_{pq}=\Big{(}\frac{4\pi}{3}\Big{)}^{2}(p^{2}+q^{2}-pq),\
1\leq q\leq p/2,$

where $\lambda_{n}$ is a multiple eigenvalue with multiplicity 2 if $p\neq
2q$. $\lambda_{n}$ is in units ${1}/{a^{2}}$ and $a$ is the side length of the
triangle. In Figure 2, we compare these known eigenvalues with those computed
using the expansion method with varying $V_{0}$. Although we showed monotonic
convergence as $V_{0}\to\infty$ in 3, a properly chosen value would be at
about $V_{0}\approx 2.6\times 10^{6}$, depending on how many eigenvalues one
wishes to compute and the chosen domain.

Figure 2. Average relative error for the first $n$ eigenvalues for varying
values of $V_{0}^{-1}$.

We now solve for the L-shaped domain modes numerically using the expansion
method. Here, we use $\mathcal{F}^{S}_{N}=\\{\phi_{n}\\}_{n=1}^{N}$ on
$S=(0,2)^{2}$ and $\Omega=S\setminus[1,2)^{2}$ with

(4.2)
$\phi_{n}=\phi_{(n_{1},n_{2})}=\sqrt{\frac{2}{a_{1}}}\sin\mathopen{}\mathclose{{}\left(\frac{\pi}{a_{1}}n_{1}x_{1}}\right)\sqrt{\frac{2}{a_{2}}}\sin\mathopen{}\mathclose{{}\left(\frac{\pi}{a_{2}}n_{2}x_{2}}\right)=\sin\mathopen{}\mathclose{{}\left(\frac{\pi}{2}n_{1}x_{1}}\right)\sin\mathopen{}\mathclose{{}\left(\frac{\pi}{2}n_{2}x_{2}}\right),$

just as in [27]. We provide examples of the computed eigenmodes in Figure 3.
This L-shaped domain is common in the literature, as it is a simple
construction of a domain with no closed form solution [37, 40, 21]. In Table
I, we provide the computed eigenvalues corresponding to the provided
eigenmodes, along with those computed using a second-order finite difference
operator (see [30]), with uniform grid spacings $h$ for both the $x$ and $y$
axes:

$\displaystyle H^{\operatorname{FD}}$ $\displaystyle=-(I\otimes
D^{2}_{x}+D_{y}^{2}\otimes I)+V_{\operatorname{d}}$ $\displaystyle D_{x}^{2}$
$\displaystyle=D_{y}^{2}=\frac{1}{h^{2}}\mathopen{}\mathclose{{}\left(\begin{array}[]{cccc}-2&1&&\\\
1&\ddots&\ddots&\\\ &\ddots&\ddots&1\\\ &&1&-2\end{array}}\right)$
$\displaystyle V_{\operatorname{d}}$
$\displaystyle=\mathopen{}\mathclose{{}\left(\begin{array}[]{cccccccc}v_{11}&&&&&&&\\\
&v_{21}&&&&&&\\\ &&\ddots&&&&&\\\ &&&v_{N1}&&&&\\\ &&&&v_{22}&&&\\\
&&&&&\ddots&&\\\ &&&&&&&v_{NN}\end{array}}\right)$ $\displaystyle v_{nm}$
$\displaystyle=V(x_{n},y_{m}).$

The results in Table I of the computation of Laplacian eigenvalues for the
L-shaped domain using the expansion method and FD are compared with the known
values from [37], which were remarkably computed with up to 8 digits of
accuracy. This comparison provides numerical validation of the expansion
method.

Figure 3. Various states $|\psi|$ for the L-shaped region on a plane, computed
using the expansion method.

$n$ | FD Method | Exp. Method | Known Values
---|---|---|---
1 | 9.33328 | 9.63359 | 9.63972
2 | 14.8927 | 15.1964 | 15.1973
3 | 19.4634 | 19.7385 | 19.7392
4 | 29.2480 | 29.5209 | 29.5215
5 | 31.2302 | 31.8982 | 31.9126
6 | 40.3540 | 41.4629 | 41.4745
20 | 98.9878 | 101.585 | 101.605
50 | 239.487 | 250.777 | 250.785
104 | 462.102 | 493.543 | 493.480

$n$ | FD Method | Exp. Method | Known Values
---|---|---|---
1 | 9.69329 | 10.1213 | 9.63972
2 | 14.9296 | 15.7156 | 15.1973
3 | 19.4097 | 20.1868 | 19.7392
4 | 28.5903 | 29.8785 | 29.5215
5 | 31.171 | 32.8286 | 31.9126
6 | 39.4412 | 42.9975 | 41.4745
20 | 91.6147 | 105.747 | 101.605
50 | 179.955 | 255.717 | 250.785
104 | 306.101 | 514.121 | 493.480

Table I. Computed Laplacian eigenvalues of the L-shaped region. Both schemes
used a discretized 2-D Schrödinger operator in $\mathbb{R}^{N\times N}$ and
$V_{0}=2.1\times 10^{5}$. $N=2500$ (left) and $N=225$ (right). The known
values were taken from [37].

##  5\. Applications

###  5.1. Spherical domains.

The expansion method can be used on a variety of manifolds, and for example,
on a spherical surface. The eigenfunctions of the Laplace–Beltrami operator on
a sphere are the spherical harmonics $Y$, which are solutions to

(5.1)
$\frac{1}{\sqrt{|g|}}\partial_{i}\mathopen{}\mathclose{{}\left(\sqrt{|g|}g^{ij}\partial_{j}Y(\theta,\phi)}\right)=\lambda
Y(\theta,\phi)$

where

(5.2) $g=\mathopen{}\mathclose{{}\left[\begin{array}[]{cc}1&0\\\
0&\sin^{2}\theta\end{array}}\right].$

Spherical harmonics $Y_{m\ell}$ provide a set of orthonormal functions and
thus can be used as a basis. These functions are defined over the indices $m$
(integers) and $\ell$ (non-negative integers), where $Y_{m\ell}$ is defined
for $|m|\leq\ell$. These functions are known explicitly ($P^{m}_{\ell}$
denoting associated Legendre polynomials),

(5.3)
$Y_{m\ell}=\mathopen{}\mathclose{{}\left\\{\begin{array}[]{ll}(-1)^{m}\sqrt{2}\sqrt{\frac{2\ell+1}{4\pi}\frac{(\ell-|m|)!}{(\ell+|m|)!}}P_{\ell}^{|m|}(\cos\theta)\sin(|m|\varphi)&\text{
if }m<0\\\ \sqrt{\frac{2\ell+1}{4\pi}}P_{\ell}^{m}(\cos\theta)&\text{ if
}m=0\\\
(-1)^{m}\sqrt{2}\sqrt{\frac{2\ell+1}{4\pi}\frac{(\ell-m)!}{(\ell+m)!}}P_{\ell}^{m}(\cos\theta)\cos(m\varphi)&\text{
if }m>0.\end{array}}\right.$

The eigenfunctions of the Laplace–Beltrami operator with Dirichlet boundary
conditions for some smooth region on a sphere can be expressed in $L^{2}(S)$
as linear combinations of spherical harmonics,

(5.4) $\psi_{j}(\phi,\theta)=\sum_{\begin{subarray}{c}m,\
\ell\end{subarray}}c^{(j)}_{m\ell}Y_{\ell}^{m}(\phi,\theta).$

As in (2.7), the matrix representation of the Schrödinger operator $\hat{H}$
in the space composed of the basis functions is given by

(5.5) $H_{ij}=\int_{S}Y_{i}(\phi,\theta)\hat{H}Y_{j}(\phi,\theta)ds,$

where $i$ and $j$ each represent an index pair $(m,\ell)$. By substituting
$\mathopen{}\mathclose{{}\left[-\Delta_{g}+V(\phi,\theta)}\right]$ for
$\hat{H}$ and making a change of variables, we obtain $ds=r\sin(\theta)d\phi
d\theta$ on the unit sphere and the discretized Hamiltonian,

(5.6) $\displaystyle H_{N_{ij}}$ $\displaystyle=\langle
Y_{i},\hat{H}Y_{j}\rangle$
$\displaystyle=\ell(\ell+1)\delta_{ij}+V_{0}\int_{S}Y_{i}(\phi,\theta)Y_{j}(\phi,\theta)\sin{(\theta)}d\phi
d\theta,$

for a large value $V_{0}\gg 1$. We can then calculate the matrix $H_{N}$ and
its eigenpairs numerically. We then expand the eigenvectors back into the
spherical harmonic basis.

Using the presented method, we have calculated and plotted the first twelve
states for the half-sphere, octant, and spherical square alongside their
planar analogs (in grayscale) to illustrate the utility of the expansion
method. We have plotted the absolute value to distinguish the nodal lines.

Figure 4. The first 12 modes of a half-sphere (above) compared to those of a
planar disk (below). Notice the ordering of the degenerate modes does not
match.

Figure 5. The first 12 modes of the spherical octant (above) compared with
those of an equilateral triangle (below). Notice the ordering of the
degenerate modes does not match.

Figure 6. The first 12 modes of the hexant or spherical square (above) and
those of a planar square (below).

###  5.2. Periodic domains.

The Schrödinger equation (2.1) in periodic domains has important application
to solving for Bloch states of electrons in a crystalline solid [28]. Here we
present an extension of the expansion method to computing eigenvalues and
eigenfunctions of the Schrodinger and Helmholtz equations for periodic domain.

######  Definition 11.

A $d$-dimensional lattice $\Gamma_{B}$ is the set $\\{\sum
x_{i}{b}_{i}:x_{i}\in\mathbb{Z}\\}$ with vectors
$\\{{b}_{i}\\}\in\mathbb{R}^{d}$. Consequently, $\Gamma_{B}=\mathbb{Z}^{d}B$
for $B\in GL(d,\mathbb{R})$, the group of $d\times d$ invertible real
matrices.

######  Definition 12.

The dual of a lattice $\Gamma_{B}$ is
$\Gamma_{B}^{*}=\\{{x}\in\mathbb{R}^{d}:\langle{x},{y}\rangle\in\mathbb{Z},\forall{y}\in\Gamma_{B}\\}$.
Consequently, $\Gamma_{B}^{*}=\mathbb{Z}^{d}B^{-T}$ for $B\in
GL(d,\mathbb{R})$.

######  Definition 13.

A $d$-dimensional flat torus $T_{B}$ is defined as the quotient space
$T_{B}={\mathbb{R}^{d}}/{\Gamma_{B}}$ for $B\in GL(d,\mathbb{R})$.

Figure 7. Fundamental region of a lattice in $\mathbb{R}^{2}$ (left). A
$2$-dimensional flat torus is formed by periodic boundary conditions over
opposite edges of the fundamental region (middle). An example of a periodic
domain—the dark regions indicate holes with Dirichlet boundary conditions
where $V_{0}\gg 1$ (right).

For Euclidean space $\mathbb{R}^{d}$, eigenfunctions for $-\Delta$ on
$\Omega=T_{B}$ are of the form $\phi_{n}({x})=\exp(2\pi
i\langle{x},{w}\rangle)$ for $w\in\Gamma^{*}_{B}$, and the eigenvalues are
$-\Delta\phi_{n}({x})/\phi_{n}({x})=4\pi^{2}\|{w}\|^{2}$. We extend the
expansion method to these domains using these eigenpairs and the potential
$V(x)=V_{0}\chi_{\Omega^{c}}(x)$. This allows us to compute eigenfunctions of
domains with mixed Dirichlet and periodic boundaries. In Figure 7, we show an
example of a periodic domain with a hole removed in each cell.

Just as in previous sections, the expansion method for a flat torus is given
by the discretization of $\hat{H}$ where

(5.7) $\begin{array}[]{ll}H_{nm}&=\langle\phi_{n},\hat{H}\phi_{m}\rangle\\\
&=\lambda_{n}(T_{B})\delta_{nm}+V_{0}{\int}_{T_{B}\setminus\Omega}\phi_{n}^{*}\phi_{m}dx\end{array}$

and $\lambda_{n}(T_{B})$ denotes the eigenvalues of the flat torus $T_{B}$
itself. The eigenvalues and eigenvectors of $H_{nm}$ are approximations of the
eigenvalues and eigenfunctions of $\hat{H}$ (in the basis $\\{\phi_{n}\\}$).
Figure 8 displays computed eigenmodes on a periodic domain with holes (domain
shown in Figure 13).

Figure 8. The absolute value of the first 12 computed eigenfunctions for a
periodic domain (same domain shown in Figure 13). We impose periodic
conditions on the boundary of the square cells and Dirichlet conditions along
the boundary of the omitted shape.

###  5.3. Spectral Clustering and Billiard Trajectories

Here, we give several examples of implementing the generalized expansion
method to heuristically explore quantum signatures of chaos in classical
billiards from computed eigenvalue statistics. This is similar to what is done
in [27, 2], but on manifolds. We use the following conjectures as foundations
for the heuristic.

######  Conjecture 14.

(Berry-Tabor) The spectral value spacings of generic integrable systems
coincide with those of uncorrelated random numbers from a Poisson process.

######  Conjecture 15.

(Bohigas-Gianonni-Schmit) The spectral value spacings of generic classically
chaotic systems coincide with those of random matrices from the Gaussian
Ensembles.

Figure 9. These billiards give simple examples of chaotic trajectories (left)
and regular trajectories (right). The chaotic billiard here is the well-known
Sinai Billiard.

It is understood that generic systems have zero probability to have
symmetries, and Hamiltonians with symmetries have degenerate states that are
not relevant to the discussion of billiard dynamics [10]. Hence, in this
discussion and in common practice, we desymmetrize the domains as much as
possible before computing eigenvalues in order to ignore the symmetric modes,
and we are left with billiards which we assume are sufficiently generic.
Further insight into these conjectures and their relation to random matrix
theory can be found in [3, 35].

Now, by considering the normalized distribution of the computed first $n$
spacings between consecutive Laplace–Beltrami eigenvalues of some
desymmetrized region $\Omega$ with Dirichlet boundary conditions, we apply the
Conjectures 14 and 15 as a heuristic to verify the trajectory type (see Figure
9) of the following regions by comparing these distributions to the Poisson
distribution $P_{0}(s)=e^{-s}$ and GOE distribution
$P_{\operatorname{GOE}}(s)=\frac{1}{2}\pi se^{-\pi s^{2}/4}$.

#### Planar billiards.

As the Sinai Billiard is a well-known classically chaotic billiard, we have
used this domain as an example to perform the generalized expansion method and
compare the resulting eigenvalue spacing distribution to the expected GOE
distribution. As the equilateral triangle is a known integrable system, we
expect its eigenvalue spacing distribution to coincide with a Poisson
distribution. Using the method, we indeed arrive at these results and show
them in Figure 11. We perform, when possible, a terminating sequence of
desymmetrizations on the domains, as shown in 10.

Figure 10. Symmetries of the Sinai billiard (left), desymmetrized domain
embedded in a square (middle), and nodal lines of 40th eigenstate (right).

Figure 11. The eigenvalue spacings of an equilateral triangle (left) compared
to those of the Sinai Billiard (right).

#### Spherical billiards.

Here, we consider spherical domains: octant (spherical equilateral triangle)
and octant with a hole removed as shown in Figure 12. The former does not have
a terminating series of desymmetrizations (so we leave it as is). However, we
indeed can perform a desymmetrization on the latter, as shown. The eigenvalue
distributions are compared in Figure 12. These eigenvalue statistics suggest
the domains are regular and chaotic, respectively.

Figure 12.  Spherical octant (top left), octant with hole (top middle), and
desymmetrized octant with hole (top right). Their respective eigenvalue
spacings (bottom).

#### Periodic billiards.

In Figure 13 we illustrate and compute eigenvalue spacings for two periodic
domains after necessary desymmetrizing. They both have eigenvalues with low
clustering and appear to take on a Poisson distribution, indicating chaotic
trajectories.

Figure 13. Four cells of two periodic domains (top) and their respective
histograms of the normalized Eigenvalue Spacing Distributions for the regions
obtained by desymmetrizing (bottom). The grey regions indicate $\Omega$, and
we impose Dirichlet boundary conditions where they meet the white regions. The
first twelve eigenfunctions for the asymmetric domain are shown in Figure 8.

##  6\. Conclusion

The Laplace–Beltrami operator is crucial to describing many physical phenomena
on manifolds, and calculation of its eigenvalues is important to many
applications involving non-Euclidean media. The generalized expansion method
described in this paper provides a straightforward approach to discretize the
Laplace–Beltrami operator to approximate its eigenmodes and eigenvalues. We
provided proofs for its spectral convergence (Theorems 7 and 9) along with
various analytic and numerical examples, including its application to studying
billiard problems on surfaces. Notable applications for this method exist in
nonlinear systems such as the study of Kerr media and Bose-Einstein
condensates where one may approximate solutions by iterating on the ground
state solution of the Schrödinger equation [15, 5, 6]. Additionally, many
applications exist in condensed matter physics. For example, as demonstrated
in Section 5.2, this method can be used to solve for Bloch states of periodic
domains defined on a lattice. Other applications include theories of 2D
materials [4, 33], superconductors [9], and types of soft matter such as
membranes [31].

### Acknowledgement

The authors acknowledge the support from Department of Mathematics at
University of Utah where this project was initialized. Elena Cherkaev
acknowledges support from the U.S. National Science Foundation through grants
DMS-1715680 and DMS-2111117. Dong Wang acknowledges the support from National
Natural Science Foundation of China (NSFC) grant 12101524 and the University
Development Fund from The Chinese University of Hong Kong, Shenzhen
(UDF01001803).

## References

  * [1] Paolo Amore  “Solving the Helmholtz equation for membranes of arbitrary shape: numerical results”  In _Journal of Physics A: Mathematical and Theoretical_ 41.26 IOP Publishing, 2008, pp. 265206
  * [2] Paolo Amore  “Spectroscopy of drums and quantum billiards: Perturbative and nonperturbative results”  In _Journal of mathematical physics_ 51.5 American Institute of Physics, 2010, pp. 052105
  * [3] YY Atas, E Bogomolny, O Giraud and G Roux  “Distribution of the ratio of consecutive level spacings in random matrix ensembles”  In _Physical review letters_ 110.8 APS, 2013, pp. 084101
  * [4] Phaedon Avouris, Tony F Heinz and Tony Low  “2D Materials”  Cambridge University Press, 2017
  * [5] W. Bao  “Ground states and dynamics of multicomponent Bose–Einstein condensates”  In _Multiscale Modeling & Simulation_ 2.2, 2004, pp. 210–236  DOI: 10.1137/030600209
  * [6] W. Bao and Q. Du  “Computing the ground state solution of Bose–Einstein condensates by a normalized gradient flow”  In _SIAM Journal on Scientific Computing_ 25.5, 2004, pp. 1674–1697  DOI: 10.1137/s1064827503422956
  * [7] Weizhu Bao  “The nonlinear Schrödinger equation and applications in Bose-Einstein condensation and plasma physics”  In _Dynamics in models of coarsening, coagulation, condensation and quantization_ 9 World Scientific River Edge, NJ, USA, 2007, pp. 141–240
  * [8] Alexander Harvey Barnett  “Dissipation in Deforming Chaotic Billiards”  Harvard University, 2000
  * [9] Oleg L Berman, Yurii E Lozovik, Sergey L Eiderman and Rob D Coalson  “Superconducting photonic crystals: Numerical calculations of the band structure”  In _Physical Review B_ 74.9 APS, 2006, pp. 092505
  * [10] Michael V Berry  “Semiclassical mechanics of regular and irregular motion”  In _Les Houches lecture series_ 36 North-Holland Amsterdam, 1983, pp. 171–271
  * [11] Michael Victor Berry and Mark Wilkinson  “Diabolical points in the spectra of triangles”  In _Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences_ 392.1802 The Royal Society London, 1984, pp. 15–43
  * [12] Timo Betcke  “The generalized singular value decomposition and the method of particular solutions”  In _SIAM Journal on Scientific Computing_ 30.3 SIAM, 2008, pp. 1278–1295
  * [13] Beniamin Bogosel  “Efficient algorithm for optimizing spectral partitions”  In _Applied Mathematics and Computation_ 333 Elsevier BV, 2018, pp. 61–75  DOI: 10.1016/j.amc.2018.03.087
  * [14] B. Bourdin, D. Bucur and É. Oudet  “Optimal partitions for eigenvalues”  In _SIAM Journal on Scientific Computing_ 31.6, 2010, pp. 4100–4114  DOI: 10.1137/090747087
  * [15] Jared C Bronski, Lincoln D Carr, Bernard Deconinck and J Nathan Kutz  “Bose-Einstein condensates in standing waves: The cubic nonlinear Schrödinger equation with a periodic potential”  In _Physical Review Letters_ 86.8 APS, 2001, pp. 1402
  * [16] S.-M. Chang, C.-S. Lin, T.-C. Lin and W.-W. Lin  “Segregated nodal domains of two-dimensional multispecies Bose–Einstein condensates”  In _Physica D: Nonlinear Phenomena_ 196.3, 2004, pp. 341–361  DOI: 10.1016/j.physd.2004.06.002
  * [17] Howard S Cohl and Ernie G Kalnins  “Fundamental solution of the Laplacian in the hyperboloid model of hyperbolic geometry”  In _arXiv preprint arXiv:1201.4406_ , 2012
  * [18] M. Conti, S. Terracini and G. Verzini  “An optimal partition problem related to nonlinear eigenvalues”  In _Journal of Functional Analysis_ 198.1, 2003, pp. 160–196  DOI: 10.1016/s0022-1236(02)00105-2
  * [19] O Cybulski and R Holyst  “Three-dimensional space partition based on the first Laplacian eigenvalues in cells”  In _Physical Review E_ 77.5, 2008, pp. 56101  DOI: 10.1103/physreve.77.056101
  * [20] Gadi Fibich  “The Nonlinear Schrödinger Equation”  Springer, 2015
  * [21] L Fox, P Henrici and C Moler  “Approximations and bounds for eigenvalues of elliptic operators”  In _SIAM Journal on Numerical Analysis_ 4.1 SIAM, 1967, pp. 89–102
  * [22] David Gilbarg, Neil S Trudinger, David Gilbarg and NS Trudinger  “Elliptic Partial Differential Equations of Second Order”  Springer, 1977
  * [23] Daniel Grieser and Svenja Maronna  “Hearing the shape of a triangle”  In _Notices of the AMS_ 60.11, 2013, pp. 1440–1447
  * [24] Boling Guo, Zaihui Gan, Linghai Kong and Jingjun Zhang  “The Zakharov System and Its Soliton Solutions”  Springer, 2016
  * [25] Antoine Henrot  “Extremum Problems for Eigenvalues of Elliptic Operators”  Springer Science & Business Media, 2006
  * [26] Jeremy G Hoskins, Vladimir Rokhlin and Kirill Serkh  “On the numerical solution of elliptic partial differential equations on polygonal domains”  In _SIAM Journal on Scientific Computing_ 41.4 SIAM, 2019, pp. A2552–A2578
  * [27] D Kauffman, I Kosztin and K Schulten  “Expansion method for stationary states of quantum billiards”  In _Am. J. Phys_ 67, 1999, pp. 133–141
  * [28] Charles Kittel, Paul McEuen and Paul McEuen  “Introduction to Solid State Physics”  Wiley New York, 1996
  * [29] James R Kuttler and Vincent G Sigillito  “Eigenvalues of the Laplacian in two dimensions”  In _Siam Review_ 26.2 SIAM, 1984, pp. 163–193
  * [30] Randall J LeVeque  “Numerical Methods for Conservation Laws”  Springer, 1992
  * [31] Jannik C Meyer, Andre K Geim, Mikhail I Katsnelson, Konstantin S Novoselov, Tim J Booth and Siegmar Roth  “The structure of suspended graphene sheets”  In _Nature_ 446.7131 Nature Publishing Group, 2007, pp. 60–63
  * [32] B. Osting, C.. White and É. Oudet  “Minimal Dirichlet energy partitions for graphs”  In _SIAM J. Scientific Computing_ 36.4, 2014, pp. A1635–A1651  DOI: 10.1137/130934568
  * [33] JT Paul, AK Singh, Zheng Dong, Houlong Zhuang, BC Revard, B Rijal, M Ashton, A Linscheid, M Blonsky and D Gluhovic  “Computational methods for 2D materials: discovery, property characterization, and application design”  In _Journal of Physics: Condensed Matter_ 29.47 IOP Publishing, 2017, pp. 473001
  * [34] Martin Reuter, Silvia Biasotti, Daniela Giorgi, Giuseppe Patanè and Michela Spagnuolo  “Discrete Laplace–Beltrami operators for shape analysis and segmentation”  In _Computers & Graphics_ 33.3 Elsevier, 2009, pp. 381–390
  * [35] S Harshini Tekur and MS Santhanam  “Symmetry deduction from spectral fluctuations in complex quantum systems”  In _Physical Review Research_ 2.3 APS, 2020, pp. 032063
  * [36] Gerald Teschl  “Mathematical Methods in Quantum Mechanics”  In _Graduate Studies in Mathematics_ 99 American Mathematical Society Providence, RI, USA, 2009, pp. 106
  * [37] Lloyd N Trefethen and Timo Betcke  “Computed eigenmodes of planar regions”  In _Contemporary Mathematics_ 412 Providence, RI: American Mathematical Society, 2006, pp. 297–314
  * [38] Dong Wang  “An efficient unconditionally stable method for dirichlet partitions in arbitrary domains”  In _SIAM Journal on Scientific Computing_ 44.4 SIAM, 2022, pp. A2061–A2088
  * [39] Dong Wang and Braxton Osting  “A diffusion generated method for computing Dirichlet partitions”  In _Journal of Computational and Applied Mathematics_ 351 Elsevier BV, 2019, pp. 302–316  DOI: 10.1016/j.cam.2018.11.015
  * [40] Quan Yuan and Zhiqing He  “Bounds to eigenvalues of the Laplacian on L-shaped domain by variational methods”  In _Journal of computational and applied mathematics_ 233.4 Elsevier, 2009, pp. 1083–1090
  * [41] Vladimir E Zakharov  “Stability of periodic waves of finite amplitude on the surface of a deep fluid”  In _Journal of Applied Mechanics and Technical Physics_ 9.2 Springer, 1968, pp. 190–194