# Complex field reversal dynamics in nanomagnetic systems

Michael Saccone  Los Alamos National Laboratory, Los Alamos, New Mexico 87545,
USA  Francesco Caravelli  Los Alamos National Laboratory, Los Alamos, New
Mexico 87545, USA

###### Abstract

Nanomagnetic materials, built from thin, patterned films of ferromagnetic
materials, began as analogues to frustrated magnetism. Their low energy of
operation and emergent properties make them strong candidates for physics
based devices. A recent model of how nanomagnetic domains flip, the Glauber
mean-field model, is used here to understand how systems of nanomagnets evolve
when opposed by external field. This reversal can be expressed in an
analytical form in the case of one-dimensional chains and trees at zero
temperature, where the cascade of spin flips gives rise to harmonic power
spectra. The same cascades in two and three dimensions form fractal field
reversal clusters whose shape depends on the strength of the field and the
tuning of interactions between nanomagnets.

††preprint: APS/123-QED

##  I  Introduction

To engineer physics that goes beyond that of constituent materials, scientists
manufacture metamaterials using traditional and nanofabrication techniques
Zheludev and Kivshar (2012); Kadic  _et al._ (2019). Some metamaterials pursue
a specific functionality impossible in pure materials, such as subwavelength
focusing and invisibility cloaking Hess  _et al._ (2012), while others seek to
emulate other physical systems with greater control over the parameter space
and improved imaging of the system microstates Schiffer and Nisoli (2021);
Skjærvø  _et al._ (2020). More finds new ways to be different and the gestalt
behaviors of these systems do not perfectly emulate their supposed
counterparts due to disorder Budrikis  _et al._ (2012), dynamical differences
Morley  _et al._ (2017), or entirely new degrees of freedom Gartside  _et al._
(2022) that prove non-negligible.

In particular, the field of artificial spin ice Schiffer and Nisoli (2021);
Skjærvø  _et al._ (2020) has begun to grapple with this concept. Initially a
means of directly imaging patterned, Ising-like nanomagnets with dipolar
interactions that map onto problems in statistical physics and frustrated
magnetism, the field has since grown to encompass device-oriented approaches
to computation and evaluate the collective behavior of nanomagnets beyond
simple, Ising spins Dion  _et al._ (2022); Gartside  _et al._ (2022); Vidamour
_et al._ (2022). Visualizing the nanomagnets in real time revealed that their
fluctuations do not purely correspond to a thermal ensemble, but rather
incorporate the complexities of relaxation pathways Arava  _et al._ (2019);
Farhan  _et al._ (2013), system topology Lao  _et al._ (2018), deviation from
ergodicity Lammert  _et al._ (2010); Lao  _et al._ (2018), and innate material
properties Morley  _et al._ (2017); Drisko  _et al._ (2015). While this may be
modeled by micromagnetic simulations of the LLG equations Velo  _et al._
(2020), this computational approach is prohibitively costly for large systems
and does not provide an analytically tractable means of understanding
nanomagnet behavior. A recently introduced model of mean-field Glauber
dynamics Saccone  _et al._ (2023) represents the major magnetic domains in a
nanomagnetic system as mean-field, continuous variables between -1 and 1 (Fig.
1a) and parameterizes the interaction strength between domains. This reduces
the spatial resolution of the full micromagnetic simulation for an assumed set
of relevant variables. While this hides details such as the curving of
magnetization at the tips of the nanomagnets, the computational overhead is
significantly lowered and the increased symmetry admits insightful analytical
solutions. Within this work, we leverage this model to explain why
nanomagnetic systems deviate from simple Ising systems and predict what
untapped dynamics should emerge under the correct conditions. Pursuant to
prior experimental studies Bingham  _et al._ (2021); Hallén  _et al._ (2022);
Zeissler  _et al._ (2016); Hügli  _et al._ (2012), we begin with avalanches of
spin flips in one-dimensional chains and trees (Fig. 1b), then expand to the
more realistic two-dimensional square (Fig. 1c) and diamond lattice (Fig. 1d)
systems to find complex fractal avalanches, similar to previous numerical
experiments Chern  _et al._ (2014), and harmonic spectra. Avalanche dynamics
have long been considered crucial to information processing in the brain
Shriki  _et al._ (2013) and helpful in generating neuromorphic computation in
nanostructed networks Mallinson  _et al._ (2019); Pike  _et al._ (2020). It is
likely that spin ice avalanches can similarly enhance information processing.

Figure 1: a) The mean-field model sums over the thermal ensemble of spins in a
nanomagnetic domain and represents their average moment as a single variable
between -1 and 1, $m_{i}(t)$. b) An illustration of a tree of coupled
nanomagnets. c) The square lattice, a two-dimensional model system for
frustrated magnetism. In the nearest neighbor approximation, nanomagnets
sharing a vertex interact with one another. d) A unit cell of the diamond
lattice, a three-dimensional system where interaction occurs between tetrapod-
sharing nanomagnets.

##  II Nanomagnet chain and trees

The model explored in this paper is that of mean-field Glauber dynamics. As
derived previously for spin-ice Saccone  _et al._ (2023) as a generalization
of the Contucci model Contucci and Ghirlanda (2007); Krapivsky  _et al._
(2010), each nanomagnetic island is represented by a mean-field spin,
$m_{i}\in[-1,1]$, that evolves over time as governed by

$\dot{\vec{m}}=-\vec{m}+\tanh(\beta Q\vec{m}+\beta\vec{h}(t)).$ (1)

The magnets are subject to an interaction matrix $Q$ and external field
$\vec{h}(t)$, both modulated by an inverse temperature parameter $\beta$. In
this deunitized result, the natural relaxation time of a domain is 1. While
previous work has determined the attractors of this model and their stability
for small systemsSaccone  _et al._ (2023), here we expand the analysis to
extensive systems with a focus on their response to the external field.

In the simplest extensive system, a one-dimensional chain, the elements of the
interaction matrix are

$Q_{ij}=q\delta_{ij}+J(\delta_{i,j+1}+\delta_{i,j-1}),$ (2)

where $q$ is the self-energy of the system which penalizes deviation from
$m_{i}=\pm 1$ and $J$ is the ferromagnetic (+) or antiferromagnetic (-)
interaction strength of the neighboring magnets. A chain is a special case of
a tree (Fig. 1b) with consistent degree 2, and many of the statements made
about magnetization dynamics on a physical chain of nanomagnets can be
extended to general trees. By approximating the governing differential
equations in different regimes, we can gain insight into how the system
responds to external field. For example, the high temperature or paramagnetic
case linearizes the equations of motion and results in a system equivalent to
driven RL circuits (See Supplementary Information).

###  II.1 Low temperature limit

Low-temperature systems approach fully magnetized states that only change when
an island’s coercivity is exceeded by all fields acting on the island. This is
reflected in the hyperbolic tangent function approaching a sign function,
$\text{sgn}(x)$:

$\dot{\vec{m}}=-\vec{m}+\text{sgn}(Q\vec{m}+\vec{h}).$ (3)

Changes in the moments occur in piecewise functions determined by the argument
of the sign function changing from positive to negative, or vice versa. The
attractors of the individual magnets are the binary values $\pm 1$. The
duration of the piecewise periods may be generally difficult to determine but
is solvable for simple scenarios such as a tree interaction structure (See
Supplementary Information). The collective behavior of the system depends on
what sequence of flips occur and how they influence one another. The initial
condition of the system and the spatial pattern of the field all can evoke
different behaviors. One commonly encountered scenario is a magnetization
reversal from saturation, all magnets aligned in one direction, under a
uniform external field. We will explore how systems behave in this deceptively
simple scenario for the remainder of the paper.

For a nanomagnet chain, assume that all spins begin pointing to the left, the
system is ferromagnetic $J>0$, and that the system is subject to a static
field pointing to the right, $h>0$. The field is great enough to overcome the
coercivity of a single spin, but not great enough to flip all spins at once
considering their local field from ferromagnetic interactions is $2J$ when
they all point to the left. Thus, $q<h<2J$, or for a scaled field
$h^{*}=(h-q)/2J$, $0<h^{*}<1$. Let a single nanomagnet at location $i$ start
at time zero with a moment $m_{i}(0)=1$ to nucleate an avalanche. This is an
idealized case of a single nanomagnet flipping before its neighbors due to
disorder in fabrication Budrikis  _et al._ (2012). This magnet reverses its
field on the moments at $i\pm 1$. Those nanomagnets flip their neighbors at
$i\pm 2$, which continue the chain until the flips terminate. Each $m$
increases monotonically once it is associated sign function potential changes
sign as all arguments to the sign function only increase in magnitude
throughout this process, bringing all moments into alignment with the field
and their neighbors. This logic holds on a general tree as well, with
magnetization flips propagating from an initial site outwards with fixed time
delays determined by the strength of the field (See Supplementary Information
for the full derivation). Specifically, we find that the time between moment
flips for the uniform field is

$\tau=-\log(h^{*}).$ (4)

The delay is consistent for a given field, meaning that the strength of the
signal neither grows nor decays, precisely the ”branching fraction” of 1
determined to be optimal in neuron avalanches Shriki  _et al._ (2013). This
time delay slowly grows infinite as the field shrinks to the threshold of
$h^{*}=0$ and becomes zero as the field approaches $h^{*}=1$, implying
infinite speed as the slightest movement of a neighboring magnet will initiate
a flip. It is worth noting that this is a form of superluminal motion, in
which successive events propagate at a rate that appears faster than the speed
of light to an observer Ginzburg (1960). This apparent superluminal motion is
due to the fact that the external field acts on all the spins effectively
independently, and $h^{*}$ corresponds exactly to the value at which the
flipping stops being cooperative. The precise details of delayed vector
potentials between nanomagnets are not accounted for and may become relevant
if precise experiments near $h^{*}=1$ are conducted.

If the first flip occurs at site $j$ during $t=0$, the string propagates $R$
sites away from the initial island, and the first flipped moment is ignored
due to its lack of dynamic contribution, the total magnetization of the chain
is

$M(t)=2\sum_{k=0}^{R}\big{(}1-2e^{-(t-k\tau)}\big{)}\Theta(t-k\tau)-\Theta(k\tau-t)$
(5)

where $\Theta(x)$ is the Heaviside function. For a tree with uniform degree
$d$ (a Cayley tree), the number of flipping magnets grows exponentially with
each layer flipped from the root moment. The magnetization is then

$\displaystyle M(t)$ $\displaystyle=$ $\displaystyle
d\big{[}(1-2e^{-t}\big{)}\Theta(t)-\Theta(-t)\big{]}$ (6) $\displaystyle+$
$\displaystyle
d\sum_{k=1}^{R}(d-1)^{k}\big{(}1-2e^{-(t-k\tau)}\big{)}\Theta(t-k\tau)$
$\displaystyle-\Theta(k\tau-t).$

To connect to studies of magnetic noise spectra Dusad  _et al._ (2019) and
generally elucidate the time-series, we look to the power spectrum,
$P(\omega)=|\mathcal{F}(M(t))|.^{2}$, where $\mathcal{F}(x(t))$ is the Fourier
transform of a function. For the chain, this gives

$P(\omega)=\frac{8}{\pi}\frac{1}{\omega^{4}+\omega^{2}}\frac{1-\cos
R\tau\omega}{1-\cos\tau\omega}.$ (7)

Figure 2: The power spectra $P(f)$ for $R=10$ plotted as the external field
$h^{*}$ is varied for a (a) chain of flips and (b) a Cayley tree with $d=3$.

Holistically, the power spectrum is enveloped by
$\frac{1}{(\omega^{4}+\omega^{2})}$, quickly approaching
$\frac{1}{\omega^{4}}$, which is unitless because the natural relaxation
frequency of a single nanomagnet is set to 1 in the governing differential
equations. A signal like this would be hidden in magnetization noise due to
the $1/\omega^{4}$ envelope, masked by the stronger $1/\omega^{2}$ signals
from Brownian motion at higher temperatures or stronger still subdiffusive
motion of monopoles at lower temperatures Nisoli (2021). The driven motion of
the magnetization motion has the same $1/\omega^{4}$ signature as
superdiffusivity Nisoli (2021). The rest of the function is modulated by time
constants $\tau$ and $R\tau$, giving it a harmonic structure that is audible
when converted to the audio output (see Supplementary Material) and visible
when the spectra are plotted (Fig. 2). The positions of these harmonics depend
on the strength of the external field, as seen in Fig. 2a), and vary most
rapidly around $h^{*}=1$, the transition to collective spin flipping.
Increasing the number of spin flips adds more ripples between primary
harmonics. The form of the Cayley tree power for $d=3$ is similar outside of
more complex harmonic terms:

$P(\omega)=\frac{18}{\pi}\frac{1}{\omega^{4}+\omega^{2}}\frac{1+4^{R+1}-2^{R+2}\cos(R+1)\tau\omega}{5-4\cos\tau\omega}.$
(8)

The scaling envelope is the same, but the power is much larger due to the
exponential growth with increasing $R$, while the harmonics are proportionally
weaker (Fig. 2b)

Note that multiple of these processes may occur at the same or staggered times
and they will be no different aside from the early collision of domain walls
with one another. The specific site chosen to begin the propagation may occur
experimentally due to variance in the coercivity of individual nanomagnets or
spatial fluctuation of the field. When the field is increased gradually, this
is more reasonable to expect than the simultaneous switching of all
nanomagnets. As will be discussed more in the conclusion, this also explains
why experiments show chains of nanomagnets flipping in rapid succession
despite system dynamics still occurring on a longer timescale. Hypothetically,
a thermal flip may encourage a series of ”kinetic” flips, shown here to
require no addition of energy from the environment, which then terminate due
to a higher local coercivity.

Figure 3: The fractal dimension of the field reversal clusters at a)
$\beta=2000$ and b) $\beta=200$ plotted against interaction strength ratio
$J_{1}/J_{2}$ and scaled field $h^{*}$. Increasing the resolution of the
diagram reveals that the jagged structure of the diagram in c) rapidly becomes
smoother as temperature increases in d). e) At low temperature ($\beta=2000$),
intermediate external field ($h^{*}=0.40$), and equivalent interaction
strengths ($J_{1}/J_{2}=1$), each spin’s magnetic moment $m$ in the field
reversal simulation is plotted as a heatmap. Considering the cluster to be all
moments greater than 0, the fractal dimension of this shape is 1.71. f) When
the field is increased to $h^{*}=0.90$ ceteris paribus, the fractal dimension
remains the same while moments tend to flip at lower values of their
neighboring moments. g) At an increased temperature ($\beta=200$),
intermediate field ($h^{*}=0.45$), and on the edge of the interaction strength
ratio that results in fractal domains ($J_{1}/J_{2}=0.9827$), the moments tend
to reverse in single strings but occasionally bifurcate into more fractious
structures, resulting in a fractal dimension of 1.65.

##  III Square Lattice

A particularly captivating feature of spin ice is vertex frustration: the
symmetry of interaction energies where spins meet at a vertex means they may
take on several orientations while maintaining the same local energy.
Artificial square ice connects vertices on a square grid with nanomagnetic
edges Möller and Moessner (2006). The ground state, when all interactions are
symmetric and local, arranges spins into two-in/two-out configurations (Fig.
1c). The magnets in this pattern may be labelled by $i_{x}$ and $i_{y}$, each
with the range $[0,\sqrt{N}-1]$, that correspond to $i=i_{x}+\sqrt{N}i_{y}$,
making the interaction matrix

$Q_{ij}=q\delta_{ij}+J_{1}(\delta_{i_{x},j_{x}+1}\delta_{i_{y},j_{y}}+\delta_{i_{x},j_{x}-1}\delta_{i_{y},j_{y}}-\delta_{i_{x},j_{x}}\delta_{i_{y},j_{y}+1}\\\
-\delta_{i_{x},j_{x}}\delta_{i_{y},j_{y}-1})+J_{2}(\delta_{i_{x},j_{x}-1}\delta_{i_{y},j_{y}-1}+\delta_{i_{x},j_{x}+1}\delta_{i_{y},j_{y}+1}),$
(9)

when $i_{x}+i_{y}$ is even and

$Q_{ij}=q\delta_{ij}+J_{1}(\delta_{i_{x},j_{x}+1}\delta_{i_{y},j_{y}}+\delta_{i_{x},j_{x}-1}\delta_{i_{y},j_{y}}-\delta_{i_{x},j_{x}}\delta_{i_{y},j_{y}+1}\\\
-\delta_{i_{x},j_{x}}\delta_{i_{y},j_{y}-1})+J_{2}(\delta_{i_{x},j_{x}+1}\delta_{i_{y},j_{y}-1}+\delta_{i_{x},j_{x}-1}\delta_{i_{y},j_{y}+1}),$
(10)

when $i_{x}+i_{y}$ is odd. $J_{1}$ is the typically stronger interaction
energy of perpendicular magnets and $J_{2}$ is the typically weaker energy of
collinear magnets. Various experiments use fabrication techniques to modify
the ratio of $J_{1}/J_{2}$, typically reducing $J_{1}$ to restore symmetry to
the model Perrin  _et al._ (2016); Östman  _et al._ (2018); Farhan  _et al._
(2019). Despite the additional interactions, several modes of the system’s
field response are largely the same. In the high-temperature limit, the
homogeneous solution decays and the particular solution still responds site by
site to the external field. All previous conclusions about field response
being RL circuit-like hold. The collective modes in the low-temperature
approximation also do not occur below $h=q$ and flip the whole system
simultaneously above $h=q+2J$, encouraging the use of the same scaled field
$h^{*}$.

The most substantial new response lies in sequential spin flips. There are
parameter regimes in which a single string of spins flips, but the underlying
spin texture and types of allowed transitions influence the sequence of flips.
We simulate the field reversal of a saturated square ice with periodic
boundary conditions and a single, central nanomagnet with no coercivity to
seed a domain flip. Many transitions not anticipated through single spin flips
occur, such as direct type I to type II conversion when the delay between two
flips is low and the ratio of $J_{1}/J_{2}$ is close to 1. The process cannot
be simplified to evenly delayed spin flips but the cluster grows at its edge
as new transitions become possible. Flipped spins beget more flipped spins,
analogous to diffusion-limited aggregation (DLA) Daccord  _et al._ (1986);
Mathiesen  _et al._ (2006). DLA is a process of cluster growth in which
randomly walking particles are attached to an existing cluster when they
contact it. A relatively simple ruleset results in fractal structures with
characteristic fractal dimensions around 1.7 Daccord  _et al._ (1986);
Mathiesen  _et al._ (2006) when embedded in two dimensions and around 2.5 when
embedded in three dimensions Tokuyama and Kawasaki (1984). That is, the
fractal cluster only occupies a fraction of the whole space in which it
exists, but is not constrained to an integer lower dimension such as a line or
plane. To compare magnetic reversal clusters to DLA, the fractal nature was
quantified with a network approach to fractal dimension. This dimension
compares the accumulation of points on subsets of the network, $m$, which are
all a distance $r$ traversed from arbitrarily chosen central points. The
growth rate is assumed to be a power law, as is the case with area or volume,
$m=Ar^{d}$, where $A$ is the rate of scaling and $d$ is the dimension of the
cluster. $d$ is then extracted from the standard least squares fit of the
relationship between $m$ and $r$.

Figure 4:  a) The fractal dimension of the field reversal cluster on a
pyrochlore lattice plotted versus temperature $T$ and scaled field $h^{*}$. b)
The power spectra of the magnetization time-series at zero temperature for
increasing values of the field. c) Visualization of the moments on the
pyrochlore lattice after field reversal with $h^{*}=0.1$ and $T=0$ as seen
from slightly off the (011) direction and the d) (001) direction. The area of
the dots scales with the deviation of the moments from their initial value of
-1 and the color linearly scales between blue and red as the moments reverse.

In Fig. 3(a-b) we plot a heatmap of the fractal dimension against the ratio of
$J_{1}/J_{2}$ and the external field strength at $\beta=2000$ (Fig. 3a) and
200 (Fig. 3a). A complex phase diagram is revealed, showing distinct regions
where $f=$ 1, 2, and values between 1.6 to 1.8. Though the precise transitions
that govern the aggregation of a fractal domain vary wildly throughout the
phase space and accordingly change the exact value of the fractal dimension
(see Supplementary Material), the non-integer dimensions observed correspond
to additional work on diffusion-limited aggregation which can fluctuate based
on the precise assumptions of the simulation or experimental scenarios. One-
dimensional reversal occurs along chains of collinear moments when
$J_{1}/J_{2}$ is less than a critical, field-dependent value, and along
perpendicular moments when $J_{1}/J_{2}$ is above another value, the direction
of propagation reflects stronger interactions. For finite temperatures, there
is a region of high $h^{*}$ less than 1 where all spins flip simultaneously,
making a cluster of dimension 2. The exact shape of the boundary between these
phases is detailed and, as zooming in on the diagram reveals in Fig. 3c-d),
loses detail as the temperature is increased. The polygonal structure at
$\beta=2000$ (Fig. 3c) rapidly becomes smoother as temperature increases to
$\beta=200$ (Fig. 3d). Although the details of this phase may be difficult to
observe experimentally, the diagram suggests a stark transition between one-
dimensional reversal and fractal reversal, then fractal to total reversal at
critical fields when $J_{1}\neq J_{2}$, simply observable with precise field
control. The variety of shapes observable is broad, particularly near the
phase boundary. While interaction energies are equivalent ($J_{1}/J_{2}=1$),
the temperature is relatively low ($\beta=2000$), and the field is moderate
($h^{*}=0.40$), the cluster is highly symmetric with a fractal dimension of
$f=1.73$ (Fig. 3e). Increasing the field to $h^{*}=0.90$ (Fig. 3f) makes the
moments flip earlier in the motion of their neighbors, resulting in a similar
shape and fractal dimension of $f=1.71$. At a higher temperature ($\beta=200$)
and asymmetric interaction energies ($J_{1}/J_{2}=0.9827$), the cluster forms
through strings of reversals that occasionally bifurcate into branching
structures before recombining again into strings (Fig. 3g). This reduces the
fractal dimension to 1.65. Altogether, the strength of the field, temperature,
and balance of the interaction strengths can alter the fractal dimension and
holistic appearance of fractal domains.

##  IV Diamond Lattice

Geometric frustration naturally occurs on the diamond lattice Ramirez  _et
al._ (1999), a three-dimensional structure of connected tetrapods. When the
legs of the tetrapods are interacting units, such as hydrogen atoms in water
ice or exchange coupled electron orbitals with nonzero spin in magnetic
pyrochlore materials, the symmetry of the interactions on the tetrapods gives
rise to ice physics via the energy equivalence of 6 ”two-in, two-out”
configurations. This inspired the two-dimensional equivalent, the square
lattice, but unless perturbed the diamond lattice does not admit the same
asymmetry between $J_{1}$ and $J_{2}$ interaction energies. Recent interest in
pure magnetic pyrochlore crystals explores the spectra of their magnetic noise
as a function of the applied field and found that fractal structures of moment
reversal mediated by the path of effective monopoles explains features of the
noise Hallén  _et al._ (2022). Though the model presented here is not valid
for single electron spins, the noise profile of artificial spin ice has
recently been explored as a function of field strength and direction Goryca
_et al._ (2022) and nanomagnetic pyrochlore systems have recently been
constructed May  _et al._ (2021) and show promising potential to image their
moments under field reversal.

We perform another set of field reversal experiments on a pyrochlore system
with symmetric nearest-neighbor interactions to explore what structures emerge
from the magnetic moments and the spectra of the magnetization time-series.
Once more, a central, nucleating spin is set to $m_{i}(0)=1$. The scaled
field, $h^{*}=(h-q)/2J$, is ran between 0 and 1.1 and applied in the (100)
direction, opposite the initial magnetization in the (-100) direction. Just
after $h^{*}=0$ the moments adjacent to the nucleating moment can begin to
reverse at all. A unit cell of 16 moments is tiled in a 5x5x5 cubic grid with
periodic boundary conditions for a total of 2000 moments. Higher than
$h^{*}=1$ and all moments collectively reverse, immediately overcoming their
local interactions. As temperature increases, so does the potential to
collectively reverse via the softening of the moments. To characterize the
shape of the reversal cluster, its fractal dimension is calculated. As
observed in Fig. 4a), there is a stark cutoff between regions in the $h^{*}$
and $T$ space where the fractal dimension is 2.63 and 3. A dimension of 3
corresponds to collective reversal: the whole space is filled with flipped
moments. Those with a dimension of 2.63 are similar in structure to shapes
formed by diffusion-limited aggregation in three dimensions,
characteristically having a dimension of around 2.5 Tokuyama and Kawasaki
(1984). Once more, the process here is distinct from diffusion-limited
aggregation in that it is entirely deterministic, but the similarity in
fractal dimension suggests parallels in their physics. The cluster grows by
following a complex series of corridors that are forbidden or permitted based
on the initial state of unchanged, saturated moments and the evolving moments
from the growing cluster. The time-series of the magnetic moment does not
produce harmonics due to less consistent time delay but only drops off with a
$1/\omega^{2}$ scaling below $h^{*}=1$ (Fig 4b). The scaling jumps to
$1/\omega^{4}$ above $h^{*}=1$ as the system reversal is a single event. The
structures of the field reversal clusters are less diverse than those of the
square lattice reversals. As visualized at $h^{*}=0.1$ and $T=0$ in Fig 4c-d),
the fractal cluster follows largely the same path as field and temperature are
increased, until the point where the entire system flips all at once.

##  V Conclusion

The adaptation of nanomagnets for device applications Arava  _et al._ (2019);
Gartside  _et al._ (2022); Saccone  _et al._ (2022) requires a better
understanding of collective out-of-equilibrium dynamics. The dynamical modes
modeled here confirm one-dimensional moment reversal avalanches Bingham  _et
al._ (2021) can arise from the Glauber mean-field dynamics and produce
harmonic, superdiffusive spectra. The avalanches are stable in that their
velocity is constant, an ideal characteristic for information processing
Shriki  _et al._ (2013). The frequency of the harmonics is tunable by the
strength of the external field, making the system sensitive to different
strengths of inputs. When interaction strengths are near equal in the square
lattice and inherently for the pyrochlore lattice, the reversal clusters
become fractals. Fractals like these were seen in previous simulations of
artificial spin ice field reversal Chern  _et al._ (2014), but field reversal
experiments only produced one-dimensional avalanches. Referring to Figure
3a-b, our results demonstrate that the range of fields producing fractal field
reversal clusters shrinks as the ratio of interaction strengths becomes less
degenerate, suggesting the experimental window for fractal behavior was narrow
in previous studies Bingham  _et al._ (2021). Pyrochlore systems lack this
degeneracy, but their field reversal experiments constrained reversals to one
dimension due to higher interaction energy at the surface as a result of
fabrication May  _et al._ (2021). Future experiments with finer field control
and tuned interaction energies may reveal fractal field reversal clusters.
However, the details of single island reversal are coarse-grained in this
model. Real field reversal is facilitated by magnonic modes Zeissler  _et al._
(2016), potentially altering the precise field reversal phase diagrams of
experiments.

An implicit choice in constructing a nanomagnetic device is whether useful
information is encoded in a low-energy state or the dynamic response of the
system. The prevalence of energy-based logic gates for computation has
diminished and dynamic responses for reservoir computing are on the rise
Gartside  _et al._ (2022); Vidamour  _et al._ (2022). Directly driving a
system overcomes some limitations of nanomagnets, allowing them to evolve at
lower temperatures without freezing Morley  _et al._ (2017) and avoiding the
critical slowing down of glassy systems stopping dynamics Souletie and
Tholence (1985). Real, disordered materials will nucleate field reversal
clusters from relatively susceptible magnets and the complexity of the
emerging structures will depend on field duration, direction, and strength.
Initial studies on reservoir computing Jensen and Tufte (2020); Hon  _et al._
(2021) already show specific windows of field and temperature in which the
quality of the reservoir is high. Future analyses of systemic nanomagnet
dynamics will further explore what ranges of parameters, system geometries,
and the modes they facilitate are ideal for computation.

###### Acknowledgements.

We thank Cristiano Nisoli and Will Branford for their insightful discussions.
The work of F.C. and M.S. was carried out under the NNSA of the U.S., DoE at
LANL, Contract No. DE-AC52-06NA25396 (LDRD Grant No. PRD20190195), LA-
UR-23-21161.

##  Appendix A High Temperature Limit

When temperature increases, the constituent spins in nanomagnets become
paramagnetic. Mathematically, $\beta$ and its product with the interaction
strengths and external field become much less than one, allowing for the
approximation

$\dot{\vec{m}}=(\beta Q-\mathbb{1})\vec{m}+\beta\vec{h}(t).$ (11)

Recognizing this as an inhomogeneous system of linear differential equations,
the equations admit a homogeneous and particular solution. The magnetization
over time may be decomposed into a linear combination of eigenvectors of the
matrix $B=\beta Q-\mathbb{1}$, $\mathbf{v}_{k}$, with corresponding
eigenvalues $\lambda_{k}$, making the moments over time
$\vec{m}(t)=\sum_{k}C_{k}(t)\vec{v}_{k}$. The homogenous solution is

$\vec{m}_{h}(t)=\sum_{k}C_{k}(0)\vec{v}_{k}e^{\lambda_{k}t}.$ (12)

For the ring specifically, the eigenvectors and eigenvalues for the $k$th mode
are $v_{k,j}=\exp(i\frac{2\pi k}{N}j)$ and $\lambda_{k}=-1+\beta q+2\beta
J\cos(2\pi k/N)$. All eigenvalues are negative since $\beta q<<1$ and $\beta
J<<1$, meaning that all spatial modes that exist through initial conditions
decay. The decay rate is slightly faster for low wavelength modes when $J$ is
positive (ferromagnetism) and faster for long wavelength modes when $J$ is
negative (antiferromagnetism). Ultimately, these homogeneous solutions vanish
over time, removing the history of initial conditions and even the interaction
structure of the magnets themselves.

The particular solution depends, as the name suggests, on the field protocol
being applied. Considering a periodic field, it may be represented by Fourier
series decomposition: $\vec{h}(t)=\sum_{n}\vec{H}_{n}\exp(i\omega_{n}t)$,
where $\omega_{n}=2\pi n/T$, with period $T$. The particular solution is then

$\vec{m}_{p}=\beta\sum_{n=1}\vec{H}_{n}\frac{1}{\omega_{n}}\exp\bigg{(}i\omega_{n}t-i\pi/2\bigg{)}.$
(13)

The response leads the signal by a phase of $\pi/2$ and diminishes in
magnitude proportional to the frequency, exactly the behavior of an RL
circuit. Because this does not rely on the interaction matrix, we note that
solutions to the mean-field Glauber equation at a high-temperature generally
respond to external fields with uncoupled inductors like oscillations with
exponentially fading hysteresis.

##  Appendix B Low-Temperature Solution on a Tree

We begin with the differential equation

$\dot{\vec{m}}=-\vec{m}+\tanh\big{(}\beta(Q\vec{m}+\vec{h})\big{)}.$ (14)

We consider the limit $\beta\rightarrow\infty$, and we assume $Q_{ii}=0$ and
$Q_{ij}=QA_{ij}$, with $A_{ij}$ the adjacency matrix of the graph. In this
case, we have

$\dot{\vec{m}}=-\vec{m}+\text{sign}\big{(}QA\vec{m}+\vec{h}\big{)}.$ (15)

We now perform the change of variables $\vec{m}(t)=e^{-t}\vec{r}(t)$, and the
differential equation becomes

$\dot{\vec{r}}=\text{sign}\big{(}e^{-t}QA\vec{r}+\vec{h}\big{)}.$ (16)

Since $A_{ii}=0$, the solution of the differential equation above can be
written in the form

$\vec{r}(t)=\vec{r}(t_{0})+\int_{t_{0}}^{t}\sinh\big{(}e^{-s}QA\vec{r}(s)+\vec{h}\big{)}ds$
(17)

and thus

$\vec{m}(t)=e^{-(t-t_{0})}\vec{m}(t_{0})+e^{-t}\int_{t_{0}}^{t}e^{s}\text{sign}\big{(}QA\vec{m}(s)+\vec{h}\big{)}ds$
(18)

We now make the following assumptions. First, we assume that $A_{ij}$ is the
adjacency matrix of a loopless graph, e.g. a tree $\mathcal{T}$. Initially,
the system is in the all up or all down states, e.g. $\vec{m}(t_{0})=\pm 1$.
Since the equation is invariant under the transformation
$\vec{m}\rightarrow-\vec{m}$, $\vec{h}\rightarrow-\vec{h}$, we can choose an
initial state, and obtain the solution for the, say, down state $-1$, and
obtain the evolution for the other solution by flipping the sign of the
solution and of the external field. If all initial spins are in the $-1$
state, a negative external field $h$ leaves the state invariant.

We then assume that an initial spin $i$ is flipped $m_{i}(t_{0})=+1$, and
$h_{i}>0$. In the following, let us define

$F_{i}(t;a,b)\equiv-e^{-(t-t_{0})}m_{i}(t)+e^{-t}\int_{t_{0}}^{t}e^{s}\text{sign}\big{(}Qa(s)+b(s)\big{)}ds.$

Figure 5: Tree $\mathcal{T}$ and its initial perturbation, and the rooted tree
$\mathcal{T}_{i}$ at the perturbation. The end points of the tree are the
children, and from the root to a child, one has to go through a branch.

It is useful then to consider the tree $\mathcal{T}$ rooted at $i$. For a
rooted tree, we set $i=0$ and assign the following coordinates to the spins in
reference to the spin $0$. For each child branch of $0$, we call the
coordinates ${i_{1},i_{2},\cdots}$ where, $i^{1}_{j}$ is the $i^{1}$-th branch
based at the root, ${i_{1},i_{2}}$ represents the $i^{2}$ branch based at the
branch $i^{1}$ and so on and so forth. Then $m_{k;i_{1},i_{2},\cdots}$
represents the $k$-spin on the branch ${i_{1},i_{2},\cdots}$. An example is
shown in Fig. 5.

Since $m_{0}(t_{0})=1$ and $h_{0}>0$, it is easy to see that
$m_{i}(t)=\theta(t-t_{0})$ from eqn. (18), if $h_{0}>Q$. This is the condition
that perturbations are unstable. It is not hard to see that for $h>2Q$, all
spins flip independently because the magnetic field.

We then have the chain of equations

$\displaystyle m_{0}(t)$ $\displaystyle=$ $\displaystyle 1$ $\displaystyle
m_{1;j}(t)$ $\displaystyle=$ $\displaystyle F_{1;j}(t;m_{0}+m_{2;j},h_{1;1})$
$\displaystyle\vdots$ $\displaystyle m_{k_{j};j}(t)$ $\displaystyle=$
$\displaystyle F_{1;j}(t;m_{k_{j}-1;j}+\sum_{r}m_{1;j,r},h_{k_{j};j})$
$\displaystyle\vdots$ $\displaystyle m_{1;j,j^{\prime}}(t)$ $\displaystyle=$
$\displaystyle
F_{1;j,j^{\prime}}(t;m_{k_{j};j}+m_{2;j,j^{\prime}}(t),h_{1;j,j^{\prime}})$
$\displaystyle\vdots$

Since we are on a tree, what this shows is that we can consider any path from
the root to the child, and these equations are decoupled subsets each
topologically equivalent to the flipping of a spin on a chain graph, in which
the spin at the beginning or the end of the chain is flipped at the initial
state. Thus, solving it on the chain graph then provides the solution on each
branch, as in the Lila subgraph in Fig. 5, provided that we properly change
$h_{i}$ at the junction nodes, which we will discuss later. Let us first focus
on a chain graph. This is a chain of equations of the form :

$\displaystyle m_{0}(t)$ $\displaystyle=$ $\displaystyle 1$ $\displaystyle
m_{1}(t)$ $\displaystyle=$ $\displaystyle F_{1}(t,m_{0}(t)+m_{2}(t),h_{1})$
$\displaystyle m_{2}(t)$ $\displaystyle=$ $\displaystyle
F_{2}(t,m_{1}(t)+m_{3}(t),h_{2})$ $\displaystyle m_{3}(t)$ $\displaystyle=$
$\displaystyle F_{3}(t,m_{2}(t)+m_{4}(t),h_{3})$ $\displaystyle\vdots$ (20)

Let us now focus on the second equation. Since $m_{i}\geq-1$, we have the
inequality

$Q(m_{i-1}+m_{i+1})+h\geq Q(m_{i-1}-1)+h_{i},$ (21)

which we will use in a moment, since we assume $Q\geq 0$.

Inside the sign function, we have $Q(1+m_{2}(t))+h_{2}\geq h_{2}$. Since,
initially, $m_{2}(t_{0})=-1$, and we have that $\text{sign}(h_{1})=1$, it is
easy to see that assuming $t_{0}=0$ the solution is given by,

$\displaystyle m_{1}(t)=-\theta(-t)+\theta(t)(1-2e^{-t}).$ (22)

For $m_{2}(t)$, we have

$\displaystyle m_{2}(t)$ $\displaystyle=$ $\displaystyle e^{-t}m_{2}(0)$
$\displaystyle+e^{-t}\int_{0}^{t}e^{s}\text{sign}\big{(}Q(1-2e^{-t}+m_{3}(t))+h_{2}\big{)}ds$

and we have, because of the inequality (21), we have that

$\displaystyle\text{sign}\big{(}Q(1-2e^{-t}+m_{3}(t))+h\big{)}\geq\text{sign}\big{(}Q(-2e^{-t})+h_{2}\big{)}$

and thus, the sign switches when $-2Qe^{-t}+h=0$. We can define
$\tau_{2}=\log(2Q/h)$, and after a little algebra, we find the solution

$\displaystyle
m_{2}(t)=-\theta(\tau_{2}-t)+\theta(t-\tau_{2})\Big{(}1-2e^{-t+\tau_{2}}\Big{)}.$

At this point, we see that for $m_{3}$ we have the same switching condition,
but in which we replace $\tau_{2}$ with $\tau_{3}$, where

$\displaystyle\tau_{3}=\tau_{2}+\log(\frac{2Q}{h_{3}}).$ (23)

The equation above then repeats at $m_{4}$, and we have thus the same set of
solutions, defined by the recursive equations:

$\displaystyle m_{0}(t)$ $\displaystyle=$ $\displaystyle 1$ (24)
$\displaystyle m_{i}(t)=$ $\displaystyle-\theta(\tau_{i}-t)$ (25)
$\displaystyle+\theta(t-\tau_{i})\Big{(}1-2e^{-t+\tau_{i}}\Big{)}$

where

$\displaystyle\tau_{i+1}$ $\displaystyle=$
$\displaystyle\tau_{i}+\log(\frac{2Q}{h_{i+1}})$ (26) $\displaystyle\tau_{1}$
$\displaystyle=$ $\displaystyle 0.$ (27)

or, if $h_{i}=h$, we have simply $\tau_{i}=(i-1)\log(\frac{2Q}{h})$. The speed
of propagation of the signal is then given by

$\displaystyle v=\frac{\Delta x}{\Delta t}=\frac{1}{\log(\frac{2Q}{h})}.$ (28)

The speed of propagation reaches infinity for $h=2Q$, which is the case in
which every single spin is flipped independently from each other, only due to
the effect of the external field. If instead the field is non-homogeneous, the
effective speed is given by

$\displaystyle\bar{v}=\frac{\Delta x}{\Delta
t}=\frac{1}{\log(\frac{2Q}{(\prod_{j=1}^{n}h_{j})^{1/n}})}.$ (29)

The difference between the chain and a tree is the presence of node equations,
such as eqn. (B). These contain terms of the form

$\displaystyle\text{sign}\big{(}Q(m_{k_{\vec{j}};\\{\vec{j}\\}}+\sum_{r}m_{1;\\{\vec{j},r\\}})+h\big{)}.$
(30)

The equation above is an upper bound to

$\displaystyle\text{sign}\big{(}Q(m_{k_{\vec{j}};\\{\vec{j}\\}}-d_{\\{\vec{j}\\}})+1+h\big{)}$
$\displaystyle=\text{sign}\big{(}Q(1-2e^{-t+\tau_{k_{\vec{j}};\\{\vec{j}\\}}}-d_{\\{\vec{j}\\}}+1)+h\big{)}$
(31)

Thus, the key difference between the chain and the tree, is that we need to
replace

$\displaystyle h_{i}\rightarrow h_{i}-(d_{i}-2)Q$ (32)

where $d_{i}$ is the node degree. As a result, for the perturbation to
propagate across nodes, we must have $h_{i}\geq(d_{i}-2)Q$ for all nodes. For
a chain, $d_{i}=2$, and thus reduces to the chain solution above.

## References

  * Zheludev and Kivshar (2012) N. I. Zheludev and Y. S. Kivshar, Nature materials 11, 917 (2012).
  * Kadic  _et al._ (2019) M. Kadic, G. W. Milton, M. van Hecke, and M. Wegener, Nature Reviews Physics 1, 198 (2019).
  * Hess  _et al._ (2012) O. Hess, J. B. Pendry, S. A. Maier, R. F. Oulton, J. M. Hamm, and K. L. Tsakmakidis, Nature materials 11, 573 (2012).
  * Schiffer and Nisoli (2021) P. Schiffer and C. Nisoli, Applied Physics Letters 118, 110501 (2021).
  * Skjærvø  _et al._ (2020) S. H. Skjærvø, C. H. Marrows, R. L. Stamps, and L. J. Heyderman, Nature Reviews Physics 2, 13 (2020).
  * Budrikis  _et al._ (2012) Z. Budrikis, J. Morgan, J. Akerman, A. Stein, P. Politi, S. Langridge, C. Marrows, and R. Stamps, Physical review letters 109, 037203 (2012).
  * Morley  _et al._ (2017) S. Morley, D. A. Venero, J. Porro, S. Riley, A. Stein, P. Steadman, R. Stamps, S. Langridge, and C. Marrows, Physical Review B 95, 104422 (2017).
  * Gartside  _et al._ (2022) J. C. Gartside, K. D. Stenning, A. Vanstone, H. H. Holder, D. M. Arroo, T. Dion, F. Caravelli, H. Kurebayashi, and W. R. Branford, Nature Nanotechnology 17, 460 (2022).
  * Dion  _et al._ (2022) T. Dion, J. Gartside, A. Vanstone, K. Stenning, D. Arroo, H. Kurebayashi, and W. Branford, Physical Review Research 4, 013107 (2022).
  * Vidamour  _et al._ (2022) I. Vidamour, C. Swindells, G. Venkat, P. Fry, A. Welbourne, R. Rowan-Robinson, D. Backes, F. Maccherozzi, S. Dhesi, E. Vasilaki, _et al._ , arXiv preprint arXiv:2206.04446 (2022).
  * Arava  _et al._ (2019) H. Arava, N. Leo, D. Schildknecht, J. Cui, J. Vijayakumar, P. M. Derlet, A. Kleibert, and L. J. Heyderman, Physical Review Applied 11, 054086 (2019).
  * Farhan  _et al._ (2013) A. Farhan, P. M. Derlet, A. Kleibert, A. Balan, R. V. Chopdekar, M. Wyss, J. Perron, A. Scholl, F. Nolting, and L. J. Heyderman, Physical review letters 111, 057204 (2013).
  * Lao  _et al._ (2018) Y. Lao, F. Caravelli, M. Sheikh, J. Sklenar, D. Gardeazabal, J. D. Watts, A. M. Albrecht, A. Scholl, K. Dahmen, C. Nisoli, _et al._ , Nature Physics 14, 723 (2018).
  * Lammert  _et al._ (2010) P. E. Lammert, X. Ke, J. Li, C. Nisoli, D. M. Garand, V. H. Crespi, and P. Schiffer, Nature Physics 6, 786 (2010).
  * Drisko  _et al._ (2015) J. Drisko, S. Daunheimer, and J. Cumings, Physical Review B 91, 224406 (2015).
  * Velo  _et al._ (2020) M. F. Velo, B. M. Cecchi, and K. R. Pirota, Physical Review B 102, 224420 (2020).
  * Saccone  _et al._ (2023) M. Saccone, J. C. Gartside, K. D. Stenning, W. R. Branford, and F. Caravelli, Physics of Fluids 35, 017101 (2023).
  * Bingham  _et al._ (2021) N. Bingham, S. Rooke, J. Park, A. Simon, W. Zhu, X. Zhang, J. Batley, J. Watts, C. Leighton, K. Dahmen, _et al._ , Physical review letters 127, 207203 (2021).
  * Hallén  _et al._ (2022) J. N. Hallén, S. A. Grigera, D. A. Tennant, C. Castelnovo, and R. Moessner, Science 378, 1218 (2022).
  * Zeissler  _et al._ (2016) K. Zeissler, M. Chadha, E. Lovell, L. F. Cohen, and W. R. Branford, Scientific reports 6, 1 (2016).
  * Hügli  _et al._ (2012) R. Hügli, G. Duff, B. O’Conchuir, E. Mengotti, A. F. Rodríguez, F. Nolting, L. Heyderman, and H. Braun, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370, 5767 (2012).
  * Chern  _et al._ (2014) G.-W. Chern, C. Reichhardt, and C. O. Reichhardt, New Journal of Physics 16, 063051 (2014).
  * Shriki  _et al._ (2013) O. Shriki, J. Alstott, F. Carver, T. Holroyd, R. N. Henson, M. L. Smith, R. Coppola, E. Bullmore, and D. Plenz, Journal of Neuroscience 33, 7079 (2013).
  * Mallinson  _et al._ (2019) J. Mallinson, S. Shirai, S. Acharya, S. Bose, E. Galli, and S. Brown, Science advances 5, eaaw8438 (2019).
  * Pike  _et al._ (2020) M. D. Pike, S. K. Bose, J. B. Mallinson, S. K. Acharya, S. Shirai, E. Galli, S. J. Weddell, P. J. Bones, M. D. Arnold, and S. A. Brown, Nano letters 20, 3935 (2020).
  * Contucci and Ghirlanda (2007) P. Contucci and S. Ghirlanda, Quality & quantity 41, 569 (2007).
  * Krapivsky  _et al._ (2010) P. L. Krapivsky, S. Redner, and E. Ben-Naim,  _A kinetic view of statistical physics_ (Cambridge University Press, 2010).
  * Ginzburg (1960) V. L. Ginzburg, Soviet Physics Uspekhi 2, 874 (1960).
  * Dusad  _et al._ (2019) R. Dusad, F. K. Kirschner, J. C. Hoke, B. R. Roberts, A. Eyal, F. Flicker, G. M. Luke, S. J. Blundell, and J. Davis, Nature 571, 234 (2019).
  * Nisoli (2021) C. Nisoli, EPL (Europhysics Letters) 135, 57002 (2021).
  * Möller and Moessner (2006) G. Möller and R. Moessner, Physical Review Letters 96, 237202 (2006).
  * Perrin  _et al._ (2016) Y. Perrin, B. Canals, and N. Rougemaille, Nature 540, 410 (2016).
  * Östman  _et al._ (2018) E. Östman, H. Stopfel, I.-A. Chioar, U. B. Arnalds, A. Stein, V. Kapaklis, and B. Hjörvarsson, Nature Physics 14, 375 (2018).
  * Farhan  _et al._ (2019) A. Farhan, M. Saccone, C. F. Petersen, S. Dhuey, R. V. Chopdekar, Y.-L. Huang, N. Kent, Z. Chen, M. J. Alava, T. Lippert, _et al._ , Science advances 5, eaav6380 (2019).
  * Daccord  _et al._ (1986) G. Daccord, J. Nittmann, and H. E. Stanley, Physical review letters 56, 336 (1986).
  * Mathiesen  _et al._ (2006) J. Mathiesen, I. Procaccia, H. L. Swinney, and M. Thrasher, EPL (Europhysics Letters) 76, 257 (2006).
  * Tokuyama and Kawasaki (1984) M. Tokuyama and K. Kawasaki, Physics Letters A 100, 337 (1984).
  * Ramirez  _et al._ (1999) A. P. Ramirez, A. Hayashi, R. J. Cava, R. Siddharthan, and B. Shastry, Nature 399, 333 (1999).
  * Goryca  _et al._ (2022) M. Goryca, X. Zhang, J. Watts, C. Nisoli, C. Leighton, P. Schiffer, and S. Crooker, Physical Review B 105, 094406 (2022).
  * May  _et al._ (2021) A. May, M. Saccone, A. van den Berg, J. Askey, M. Hunt, and S. Ladak, Nature communications 12, 1 (2021).
  * Saccone  _et al._ (2022) M. Saccone, F. Caravelli, K. Hofhuis, S. Parchenko, Y. A. Birkhölzer, S. Dhuey, A. Kleibert, S. Van Dijken, C. Nisoli, and A. Farhan, Nature Physics 18, 517 (2022).
  * Souletie and Tholence (1985) J. Souletie and J. Tholence, Physical Review B 32, 516 (1985).
  * Jensen and Tufte (2020) J. H. Jensen and G. Tufte, in  _ALIFE 2020: The 2020 Conference on Artificial Life_ (MIT Press, 2020) pp. 376–383.
  * Hon  _et al._ (2021) K. Hon, Y. Kuwabiraki, M. Goto, R. Nakatani, Y. Suzuki, and H. Nomura, Applied Physics Express 14, 033001 (2021).


# Matrix bispectrality and noncommutative algebras: beyond the prolate
spheroidals

F. Alberto Grünbaum
Department of Mathematics
University of California, Berkeley, CA
<EMAIL_ADDRESS>
& Brian D. Vasquez
IMPA, Rio de Janeiro
<EMAIL_ADDRESS>
& Jorge P. Zubelli
Department of Mathematics, Khalifa University
<EMAIL_ADDRESS>

_Keywords_ bispectral algebras $\cdot$ bispectral triple $\cdot$ presentations
of finitely generated algebras $\cdot$ time-and-band limiting $\cdot$ integral
and commuting differential operators

##  1 Preamble

The bispectral problem discussed in this paper has its origin in Duistermaat
and Grünbaum (1986). It is motivated, as mentioned there, by an effort to
understand and extend a remarkable phenomenon in Fourier analysis on the real
line: the operator of time-and-band limiting is an integral operator admiting
a second order differential operator with simple spectrum in its commutator.
This property, which gives a good numerical way to compute the eigenfunctions
of the integral operator, was put to good use in a series of papers by D.
Slepian, H. Landau and H. Pollak at Bells Labs back in the 1960’s, see Slepian
and Pollak (1961); Landau and Pollak (1962, 1961, 1962); Slepian (1964, 1978,
1983, 1976) and is of interest in one other contribution in this issue, see
Connes and Consani (2021); Connes and Marcolli (2008); Connes and Moscovici
(2021). We are thankful to Luc Haine, of Louvain-la-Neuve, Belgium who alerted
one of us (AG) on Nov 2021 about a talk by Alain Connes in the series
Mathematicial Picture Language. This talk delivered on Dec 7th. 2021 can be
seen in Youtube and covers some of the contents of his work with various
collaborators on the zeta function and its relations to the prolate spheroidal
functions. For readers interested in this fascinating connection there is no
better way of learning about this than watching the lecture. See Connes
(2021).

If one tries to extend this property beyond the Fourier case by adding a
potential $V(x)$ to the operator $-\left(\frac{d}{dx}\right)^{2}$ and
replacing an expansion in exponentials by an expansion in term of the
eigenfunctions $\psi(x,z)$ of $-\left(\frac{d}{dx}\right)^{2}+V(x)$ it appears
plausible that a certain property of these eigenfunctions $\psi(x,z)$ could
play a useful role. This property is now known as the bispectral property
which will be formulated below. Its solution in the scalar case was the
purpose of Duistermaat and Grünbaum (1986).

The version of the bispectral problem that we discuss in this paper is a
noncommutative one, obtained by allowing all objects in the original
formulation to be matrix valued. The details are given later.

Now for a bit of history: one of us (AG) gave a couple of talks, one at
Vanderbilt in Nov 2013 and one in IMPA, Rio de Janeiro in March 2014 . The
first talk was at the invitation of Vaughan Jones and started by saying that
the topic was most likely of no interest to him. At some point in the talk
Vaughan said "all of this is about bi-modules and subfactors". The second talk
had both Vaughan and another one of us (JZ) in the audience, there was some
more mention of these topics and then in Grünbaum (2014), written in April
2014, a reference is made to a future joint paper with Vaughan "On the
bimodule structure of the bispectral problem, in preparation". The occasion of
this second talk was a visit that Vaughan (and for part of it with Wendy) did
to Argentina, Uruguay and Brazil. In Buenos Aires he delivered one lecture, in
Rio he delivered a series of talks but in Montevideo besides giving a lecture
he was received in a private audience by President Jose Mujica, described as
"the world humblest head of state" by Wikipedia. Vaughan whose command of the
spanish language was quite high lamented that he had troubles understanding
President Mujica.

The examples in this short paper Grünbaum (2014) were discussed with Vaughan
who showed some level of interest but we never managed to get him fully on
board. Of course we assumed that we always had time to get him onto this
project. The short paper Grünbaum (2014) makes three separate conjectures
going with three different examples. These conjectures were proved to be
correct in the thesis of the third of us (BVC) under the supervision of (JZ).

And now that sadly we have no chance of benefiting from the insights that
Vaughan would have brought in, we present the problem and some results to a
wider audience in the hope that someone may surmise what it was that Vaughan
had in mind. Getting someone involved in this effort would be a nice way of
honoring his memory.

##  2 Commuting integral and differential operators

The bispectral problem, introduced in the next section is motivated by the
following very concrete problem in signal communication: a signal of support
in the interval $\left[-T,T\right]$ is transmitted over a channel that has
bandwidth $\left[-W,W\right]$ i.e., all frequencies in the signal beyond
absolute value of $W$ cannot be sent over. A mathematical formulation is as
follows: an arbitrary signal in $L^{2}(\mathbf{\mathbb{R}})$ is chopped to the
interval $\left[-T,T\right]$ and then its Fourier transform is chopped to the
interval $\left[-W,W\right]$. Denoting for simplicity these two chopping
operations by $T$ and $W$ we are dealing with the operator

$E=W\mathcal{F}T$

where $\mathcal{F}$ stands for the Fourier transform. The spectral analysis of
this operator, i.e., a look at its singular functions and singular values
requires the consideration of the operator $E^{*}E$. It is easy to see that
$E^{*}E$ is an integral operator acting in $L^{2}(-T,T)$ whose kernel is given
by

$\frac{\sin W(t-s)}{t-s}$

and this bounded operator acts on a function in the space $L^{2}(-T,T)$ by

$(Kf)(s)=\int_{-T}^{T}\frac{\sin W(t-s)}{t-s}f(t)dt$

for $f\in L^{2}(-T,T)$ and $s\in(-T,T)$. This $K$ commutes with the operator

$(Lf)(x)=\left(-\frac{d}{dx}\right)(T^{2}-x^{2})\left(\frac{df}{dx}\right)+W^{2}x^{2}f(x)$

defined on $C^{2}$ functions.

One can show that this densely defined operator has a unique selfadjoint
extension in $L^{2}(-T,T)$ with eigenfunctions and eigenvalues that depend on
the parameter $W$.

Its eigenfunctions are known as the prolate spheroidal wave functions, since
this is one of the differential operators resulting in separating varables
when solving for the eigenfunctions of the Laplacian on a prolate spheroid.

What other naturally appearing integral operators allow for commuting
differential operators? Two other examples are the Bessel and Airy kernels, as
in the work of Tracy and Widom, Tracy and Widom (1994a, b) in the context of
Random Matrix Theory. For the Bessel case see also Slepian (1964); Grünbaum et
al. (1982). There are other examples, but the search is nowhere close to
finished. The bispectral property, to be formulated below was put forward in
Duistermaat and Grünbaum (1986) as an important ingredient in the search for
more examples of this commuting property.

##  3 The bispectral problem

The problem was posed and solved in Duistermaat and Grünbaum (1986). It is as
follows:

Find all nontrivial instances where a function $\psi(x,z)$ satisfies

$L\left(x,\frac{d}{dx}\right)\psi(x,z)\equiv\left(-\left(\frac{d}{dx}\right)^{2}+V(x)\right)\psi(x,z)=z\psi(x,z)$

as well as

$B\left(z,\frac{d}{dz}\right)\psi(x,z)\equiv\left(\sum_{i=0}^{M}b_{i}(z)\left(\frac{d}{dz}\right)^{i}\right)\psi(x,z)=\Theta(x)\psi(x,z).$

All the functions $V(x),b_{i}(k),\Theta(x)$ are, in principle, arbitrary
except for smoothness assumptions. Notice that here $M$ is arbitrary (finite).
The operator $L$ could be of higher order, but in Duistermaat and Grünbaum
(1986) attention is restricted to order two.

The complete solution is given as follows:

Theorem  If $M=2$, then $V(x)$ is (except for translation) either
$\frac{c}{x^{2}}$ or $ax$, i.e. we have a Bessel or an Airy case. If $M>2$,
there are two families of solutions

  * a)

$L$ is obtained from $L_{0}=-\left(\frac{d}{dx}\right)^{2}$ by a finite number
of Darboux transformations $(L=AA^{*}\rightarrow{\tilde{L}}=A^{*}A)$. In this
case $V$ is a rational solution of the Korteweg-deVries hierarchy of
equations. Here A is a first order differential operator.

  * b)

$L$ is obtained from $L_{0}=-\left(\frac{d}{dx}\right)^{2}+\frac{1}{4x^{2}}$
after a finite number of (rational) Darboux transformations.

In all cases we have a solution of the ad-conditions, a complicated system of
nonlinear differential equations. These conditions are necesary and
sufficient.

Notice that the solutions organize themselves into nice manifolds.

The simplest example of case $a)$ follows from
$L_{0}=-\left(\frac{d}{dx}\right)^{2}$ by two Darboux transformations, one
gets the operator

$L_{2}=-\left(\frac{d}{dx}\right)^{2}+\frac{6(x^{4}+12t_{3}x)}{(x^{3}-t_{3})^{2}}\
.$

In this case $\Theta(x)=x^{4}-4t_{3}x$ and the differential operator in the
spectral parameter is

$B_{2}\left(z,\frac{d}{dz}\right)=\left(-\left(\frac{d}{dz}\right)^{2}+\frac{6}{z^{2}}\right)^{2}+4it_{3}\left(\frac{d}{dz}\right)$

The potential in the operator $L_{2}=-\left(\frac{d}{dx}\right)^{2}+V(x,t)$
above satisfies the KdV equation.

It was later observed by Magri and Zubelli, see Zubelli and Magri (1991) that
in case $b)$ we are dealing with rational solutions of the Virasoro equations
(i.e. master symmetries of KdV). The bigger picture became more apparent in
the work Zubelli and Valerio Silva (2000) where it is shown that the generic
rational potentials that decay at infinity and remain rational by all the
flows of the master-symmetry KdV hierarchy are bispectral potentials for the
Schrödinger operator.

In case $a)$ the space of common solutions has dimension one and in case $b)$
it has dimension two. One refers to these as the rank one and rank two
situations.

Observe that the “trivial cases” when $M=2$ are self-dual in the sense that
since the eigenfunctions $\psi(x,z)$ are functions either of the product $xz$
or of the sum $x+z$, one gets $B$ by replacing $z$ for $x$ in $L$. The
bispectral involution introduced in Wilson (1993) shows how this can be
adapted in the "higher order cases".

##  4 The non-commutative version of the bispectral problem

A first noncommutative (or matrix) version of the bispectral problem was
considered in J. Zubelli Ph.D. thesis at Berkeley, see also Zubelli (1990,
1992)Zubelli (1992) in the situation where both the physical space and
spectral operators act on the same side of the eigenfunction and the
eigenvalues are both scalar valued. Later on, several other versions were
considered. See Sakhnovich and Zubelli (2001); Kasman (2015); Boyallian and
Liberati (2008); Grünbaum and Iliev (2003); Grünbaum (2014); Geiger et al.
(2017) and references therein. The noncommutative version of the bispectral
problem displayed interesting connections with soliton equations as well.
Indeed, in Zubelli (1992) it was shown that a large class of rational
solutions to the AKNS hierarchy Ablowitz et al. (1974) led to matrix
differential operators that displayed the bispectral property. Among the
important equations in Mathematical Physics that are covered by the AKNS
hierarchy one finds the modified KdV and the nonlinear Schrödinger equation.
The matrix differential operator that appeared in this case was in turn
related to Dirac operators. The connection between bispectrality and another
important topic in Mathematical Physics, namely Huygens’ principle in the
strict sense Berest and Veselov (1994) turned out to appear also in the
context of Matrix Bispectrality. Indeed, in Chalub and Zubelli (2006, 2005,
2004, 2001), it was shown that rational solutions to the AKNS hierarchy led to
Dirac operators which satisfy Huygens’ principle in the stric sense. In other
words, the fundamental solutions of the perturbed Dirac equation in a suitably
high space-time dimension had its support precisely on the surface of the
light cone and not in its interior. Another interesting connection between
Matrix bispectrality to soliton equations of Mathematical Physics was explored
in Sakhnovich and Zubelli (2001).

In this present paper we take the bi-module structure of the problem into
account and let the operators act on different sides as well as allow both
eigenvalues to be matrix valued.

We consider triplets $(L,\psi,B)$ satisfying the equations

$L\psi(x,z)=\psi(x,z)F(z)\hskip 28.45274pt(\psi B)(x,z)=\theta(x)\psi(x,z)$
(1)

with $L=L(x,\frac{d}{dx})$, $B=B(z,\frac{d}{dz})$ linear matrix differential
operators, i.e.,
$L\psi=\sum_{i=0}^{l}a_{i}(x)\cdot\left(\frac{d}{dx}\right)^{i}\psi$, $\psi
B=\sum_{j=0}^{m}\left(\frac{d}{dz}\right)^{j}\psi\cdot b_{j}(z)$. The
functions $a_{i}:U\subset\mathbf{\mathbb{C}}\rightarrow
M_{N}(\mathbf{\mathbb{C}})$, $b_{j}:V\subset\mathbf{\mathbb{C}}\rightarrow
M_{N}(\mathbf{\mathbb{C}}),$ $F:V\subset\mathbf{\mathbb{C}}\rightarrow
M_{N}(\mathbf{\mathbb{C}}),$ $\theta:U\subset\mathbf{\mathbb{C}}\rightarrow
M_{N}(\mathbf{\mathbb{C}})$ and the nontrivial common eigenfunction
$\psi:U\times V\subset\mathbf{\mathbb{C}}^{2}\rightarrow
M_{N}(\mathbf{\mathbb{C}})$ are in principle compatible sized meromorphic
matrix valued functions defined in suitable open subsets
$U,V\subset\mathbf{\mathbb{C}}$.

A triplet $(L,\psi,B)$ satisfying (1) is called a bispectral triplet.

The study of the structure of the algebra of possible $\theta(x)$ going with a
fixed bispectral $\psi(x,z)$ was first raised in Castro and Grünbaum (2006)
and analyzed in Tirao (2011); Grünbaum and Tirao (2007). See also Casper
(2018) and Zurrián (2017).

We consider now the examples and conjectures given in Grünbaum (2014) as well
as their validation and further description given in Vasquez and Zubelli
(2021).

For the benefit of the reader we give a few definitions before giving some
explicit results in the next section.

######  Definition 1.

Let $\mathbb{K}$ be a field, $C$ be a $\mathbb{K}$-algebra, $A$ a subring of
$C$ and $S\subset C$. We define

$A\cdot<S>=span_{\mathbb{K}}\left\\{\prod_{j=1}^{n}s_{j}\mid
s_{1},...,s_{n}\in S\cup A,n\in\mathbf{\mathbb{N}}\right\\}\mbox{ ,}$

where the noncommutative product is understood from left to right, i.e.,
$\prod_{j=1}^{n+1}s_{j}:=(\prod_{j=1}^{n}s_{j})s_{n+1},$ for $n=0,1,2,\cdots.$
For completion, $\prod_{j=1}^{0}s_{j}:=1.$

The set $A\cdot<S>$ is called the subalgebra generated by $S$ over $A$ and we
call an element $f\in A\cdot<S>$ a noncommutative polynomial with coefficients
in $A$ and set of variables $S$.

######  Definition 2.

Let $C$ be a noncommutative ring and $A$ a subring of $C$. We say that an
element $\alpha\in C$ is integral over $A$ if there exists a noncommutative
polynomial $f$ with coefficients in $A$ such that $f(\alpha)=0$. Furthermore,
we say that $\beta\in C$ is integral over $\alpha\in C$ if $\beta$ is integral
over $A\cdot<\alpha>$. Finally, $\alpha$ and $\beta$ are associated integral
if $\alpha$ is integral over $\beta$ and $\beta$ is integral over $\alpha$.

In order to characterize the algebraic structure of bispectrality in the
present noncommutative context, we start with the following definitions.

######  Definition 3.

Let $\mathbb{K}$ be a field, we denote by $\mathbb{K}\langle
x_{\lambda}\mid\lambda\in\Lambda\rangle$ the free algebra generated by the
letters $x_{\lambda}$, $\lambda\in\Lambda$ i.e.,

$\mathbb{K}\langle
x_{\lambda}\mid\lambda\in\Lambda\rangle=\bigoplus_{F\subset\Lambda,F\text{finite}}\bigoplus_{\lambda\in
F}\mathbb{K}\cdot x_{\lambda}.$

######  Definition 4.

Let $A$ be a $\mathbb{K}$-algebra. A presentation for an algebra $A$ is a
triple $(\mathbb{K}\langle x_{\lambda}\mid\lambda\in\Lambda\rangle,f,I)$ such
that $I\subset A$ is an ideal and $f:\mathbb{K}\langle
x_{\lambda}\mid\lambda\in\Lambda\rangle/I\rightarrow A$ is an isomorphism.
Furthermore, we say that $A$ is finitely generated if there exists a
presentation with $\Lambda$ finite and finitely presented if there exists a
presentation with $\Lambda$ finite and the ideal $I$ is generated by finitely
many elements.

##  5  Three examples

Take for $\Psi(x,z)$ the matrix valued function

$\Psi(x,z)=e^{xz}\begin{pmatrix}z-1/x&1/x^{2}\\\ 0&z-1/x\end{pmatrix}$

and consider all instances of matrix-valued polynomials $\theta(x)$ and
differential operators $B$ (with matrix coefficients $b_{i}(z)$) such that

$\Psi B\equiv\sum_{i=1}^{m}\left(\frac{d}{dz}\right)^{i}\Psi\
b_{i}=\theta(x)\Psi(x,z).$

In this case one has

$L\Psi=-z^{2}\Psi$

with

$L=-\left(\frac{d}{dx}\right)^{2}+2\begin{pmatrix}1/x^{2}&-2/x^{3}\\\
0&1/x^{2}\end{pmatrix}.$

In other words for this specific differential operator in the variable $x$ we
are asking for all bispectral "partners" of $L$.

One finds that one such pair $(B,\theta)$ is given by

$B=\left(\frac{d}{dz}\right)^{3}-3\left(\frac{d}{dz}\right)^{2}\cdot\frac{1}{z}+3\left(\frac{d}{dz}\right)\cdot\frac{1}{z^{2}}+3\begin{pmatrix}0&1/z^{2}\\\
0&0\end{pmatrix}$

and $\theta(x)$ the scalar-valued polynomial

$\theta(x)=x^{3}.$

The set of all possible $\theta(x)$ is given by the following subalgebra
$\mathbb{A}$. The complete statement is given by

######  Theorem 1.

Let $\Gamma$ be the sub-algebra of $M_{2}(\mathbf{\mathbb{C}})[x]$ of the form

$\left(\begin{matrix}r_{0}^{11}&r_{0}^{12}\\\
0&r_{0}^{11}\end{matrix}\right)+\left(\begin{matrix}r_{1}^{11}&r_{1}^{12}\\\
0&r_{1}^{11}\end{matrix}\right)x+\left(\begin{matrix}r_{2}^{11}&r_{2}^{12}\\\
r_{1}^{11}&r_{2}^{22}\end{matrix}\right)x^{2}+\left(\begin{matrix}r_{3}^{11}&r_{3}^{12}\\\
r_{2}^{22}+r_{2}^{11}-r_{1}^{12}&r_{3}^{22}\end{matrix}\right)x^{3}+x^{4}p(x),$

where $p\in M_{2}(\mathbf{\mathbb{C}})[x]$ and all the variables
$r_{0}^{11},r_{0}^{12},r_{1}^{11},r_{1}^{12},r_{2}^{11},r_{2}^{22},r_{3}^{11},r_{3}^{12},r_{3}^{22}\in\mathbf{\mathbb{C}}$.
Then $\Gamma=\mathbb{A}$. Moreover, for each $\theta$ we have an explicit
expression for the operator $B$.

Furthermore, we have the presentation
$\mathbb{A}=\mathbf{\mathbb{C}}\cdot\langle\alpha_{0},\alpha_{1}\mid
I=0\rangle$ with the ideal $I$ given by

$I:=\langle\alpha_{0}^{2},\alpha_{1}^{3}+\alpha_{0}\alpha_{1}\alpha_{0}-3\alpha_{1}\alpha_{0}\alpha_{1}+\alpha_{0}\alpha_{1}^{2}+\alpha_{1}^{2}\alpha_{0}\rangle\mbox{
.}$

This is an example of an algebra with an integral element over a nilpotent
one.

For the next example take for $\psi(x,z)$ the matrix-valued function

$\psi(x,z)=\left[\frac{d}{dx}-\begin{pmatrix}1/x&-1/x^{2}&1/x^{3}\\\
0&1/x&-1/x^{2}\\\ 0&0&1/x\end{pmatrix}\right]\
e^{xz}I=e^{xz}\begin{pmatrix}z-1/x&1/x^{2}&-1/x^{3}\\\ 0&z-1/x&1/x^{2}\\\
0&0&z-1/x\end{pmatrix}.$

Here one can see that

$L\psi=-z^{2}\psi$

with

$L=-\left(\frac{d}{dx}\right)^{2}+2\begin{pmatrix}1/x^{2}&-2/x^{3}&3/x^{4}\\\
0&1/x^{2}&-2/x^{3}\\\ 0&0&1/x^{2}\end{pmatrix}.$

The results in this case about the set of all possible $\theta(x)$ are given
below.

######  Theorem 2.

Let $\Gamma$ the sub-algebra of $M_{3}(\mathbf{\mathbb{C}})[x]$ of the form

$\left(\begin{matrix}r_{0}^{11}&r_{0}^{12}&r_{0}^{13}\\\
0&r_{0}^{22}&r_{0}^{23}\\\
0&0&r_{0}^{11}\end{matrix}\right)+\left(\begin{matrix}r_{1}^{11}&r_{1}^{12}&r_{1}^{13}\\\
r_{0}^{22}-r_{0}^{11}&r_{1}^{22}&r_{1}^{23}\\\
0&r_{0}^{22}-r_{0}^{11}&r_{1}^{11}+r_{0}^{23}-r_{0}^{12}\end{matrix}\right)x$
$+\left(\begin{matrix}r_{2}^{11}&r_{2}^{12}&r_{2}^{13}\\\
r_{1}^{22}-r_{1}^{11}-r_{0}^{23}+r_{0}^{12}&r_{2}^{22}&r_{2}^{23}\\\
r_{0}^{22}-r_{0}^{11}&r_{1}^{22}-r_{1}^{11}&r_{2}^{11}+r_{1}^{23}-r_{1}^{12}\end{matrix}\right)x^{2}+\left(\begin{matrix}r_{3}^{11}&r_{3}^{12}&r_{3}^{13}\\\
r_{3}^{21}&r_{3}^{22}&r_{3}^{23}\\\
r_{1}^{22}-2r_{1}^{11}-r_{0}^{23}+r_{0}^{12}&r_{3}^{32}&r_{3}^{33}\end{matrix}\right)x^{3}$
$+\left(\begin{matrix}r_{4}^{11}&r_{4}^{12}&r_{4}^{13}\\\
r_{4}^{21}&r_{4}^{22}&r_{4}^{23}\\\
r_{3}^{32}+r_{3}^{21}-r_{2}^{22}-r_{2}^{11}+r_{1}^{12}&r_{4}^{22}&r_{4}^{33}\end{matrix}\right)x^{4}$
$+\left(\begin{matrix}r_{5}^{11}&r_{5}^{12}&r_{5}^{13}\\\
r_{5}^{21}&r_{5}^{22}&r_{5}^{23}\\\
r_{4}^{32}+r_{4}^{21}-r_{3}^{33}-r_{3}^{22}-r_{3}^{11}+r_{2}^{23}+r_{2}^{12}-r_{1}^{13}&r_{5}^{32}&r_{5}^{33}\end{matrix}\right)x^{5}+x^{6}p(x)\mbox{
, }$

where $p\in M_{3}(\mathbf{\mathbb{C}})[x]$ and all the variables
$r_{0}^{11},r_{0}^{12},...,r_{5}^{33}\in\mathbf{\mathbb{C}}$ are arbitrary.

Then, $\Gamma=\mathbb{A}$ and for each $\theta$ we have an explicit expression
for the operator $B$.

Furthermore, we have the presentation
$\mathbb{A}=\mathbf{\mathbb{C}}\cdot\langle\alpha_{2},\alpha_{3}\mid
I=0\rangle$ with

$I=\langle\alpha_{2}^{3},\alpha_{3}^{2}-\alpha_{3},(\alpha_{3}\alpha_{2})^{2}\alpha_{3}-4\alpha_{3}\alpha_{2}^{2}\alpha_{3}\rangle\mbox{
.}$

This is an example of an algebra with nilpotent and idempotent associated
elements.

As the last example we consider a case when both "eigenvalues" $F$ and
$\theta$ are matrix valued. Let

$\psi(x,z)=\frac{e^{xz}}{(x-2)xz}\left(\begin{matrix}\frac{x^{3}z^{2}-2x^{2}z^{2}-2x^{2}z+3xz+2x-2}{xz}&\frac{1}{x}\\\
\frac{xz-2}{z}&x^{2}z-2xz-x+1\end{matrix}\right)$

and

$L=\left(\begin{matrix}0&0\\\
0&1\end{matrix}\right).\left(\frac{d}{dx}\right)^{2}+\left(\begin{matrix}0&\frac{1}{(x-2)x^{2}}\\\
-\frac{1}{x-2}&0\end{matrix}\right).\left(\frac{d}{dx}\right)+\left(\begin{matrix}-\frac{1}{x^{2}(x-2)^{2}}&\frac{x-1}{x^{3}(x-2)^{2}}\\\
\frac{2x-1}{x(x-2)^{2}}&-\frac{2x^{2}-4x+3}{x^{2}(x-2)^{2}}\end{matrix}\right),$

then $L\psi=\psi F$ with

$F(z)=\left(\begin{matrix}0&0\\\ 0&z^{2}\end{matrix}\right).$

It is easy to check that $\psi B=\theta\psi$ for

$B=\left(\frac{d}{dz}\right)^{3}.\left(\begin{matrix}0&0\\\
1&0\end{matrix}\right)+\left(\frac{d}{dz}\right)^{2}.\left(\begin{matrix}0&0\\\
-\frac{2z+1}{z}&0\end{matrix}\right)+\left(\frac{d}{dz}\right).\left(\begin{matrix}1&0\\\
\frac{2(z-1)}{z^{2}}&1\end{matrix}\right)+\left(\begin{matrix}-z^{-1}&0\\\
6z^{-3}&z^{-1}\end{matrix}\right)$

and

$\theta(x)=\left(\begin{matrix}x&0\\\ x^{2}(x-2)&x\end{matrix}\right).$

In this case we characterize the algebra $\mathbb{A}$ of all polynomial $F$
such that there exist $L=L\left(x,\frac{d}{dx}\right)$ with $L\psi=\psi F$ as
follows

######  Theorem 3.

Let $\Gamma$ be the sub-algebra of $M_{2}(\mathbf{\mathbb{C}})[z]$ of the form

$\left(\begin{matrix}a&0\\\
b-a&b\end{matrix}\right)+\left(\begin{matrix}c&c\\\
a-b-c&-c\end{matrix}\right)z+\left(\begin{matrix}a-b-c&c+a-b\\\
d&e\end{matrix}\right)\frac{z^{2}}{2}+z^{3}p(z),$

where $p\in M_{2}(\mathbf{\mathbb{C}})[z]$ and all the variables $a,b,c,d,e$
are arbitrary. Then $\Gamma=\mathbb{A}$.

Furthermore, we have the presentation
$\mathbb{A}=\mathbf{\mathbb{C}}\cdot\langle\theta_{1},\theta_{3},\theta_{4},\theta_{5}\mid
I=0\rangle$ with

$I=\langle\theta_{1}^{2}-\theta_{1},\theta_{4}^{2},\theta_{4}\theta_{5},\theta_{4}\theta_{1}+\theta_{4}\theta_{3}-2\theta_{4}-\theta_{5}\theta_{4}-\theta_{5}^{2},\theta_{3}^{2}-\theta_{3}+\theta_{5}-3\theta_{3}\theta_{4}\theta_{3}\theta_{5}-\theta_{1}\theta_{4}-\theta_{5}\theta_{1},$
$\theta_{3}\theta_{1}-\theta_{1}-\theta_{4}-\frac{1}{2}\theta_{4}\theta_{1}+\frac{1}{2}\theta_{4}\theta_{3}+\theta_{5}\theta_{1}-\frac{1}{2}\theta_{5}\theta_{4}+\frac{1}{2}\theta_{5}^{2}+\theta_{3}\theta_{4}-\theta_{1}\theta_{5}-\theta_{3}\theta_{5},$
$\theta_{1}\theta_{3}-\theta_{3}+\theta_{4}+\theta_{5}-\frac{3}{2}\theta_{4}\theta_{1}+\frac{3}{2}\theta_{4}\theta_{3}-2\theta_{5}\theta_{1}-\frac{3}{2}\theta_{5}\theta_{4}+\frac{3}{2}\theta_{5}^{2}+3\theta_{3}\theta_{4}+\theta_{3}\theta_{5},$
$\theta_{5}\theta_{3}-\theta_{4}\theta_{1}+\theta_{4}\theta_{3}-\theta_{5}\theta_{1}-\theta_{5}\theta_{4}+\theta_{5}^{2},\theta_{5}\theta_{1}\theta_{5}-\theta_{5}^{2}\theta_{1}-\theta_{5}\theta_{4},\theta_{5}\theta_{4}\theta_{1}-\theta_{5}^{3}+\theta_{5}\theta_{1}\theta_{4}+\theta_{5}^{2}\theta_{1},$
$\theta_{4}\theta_{1}\theta_{5}+\theta_{4}\theta_{3}\theta_{5}-\theta_{3}^{3},\theta_{5}\theta_{3}\theta_{4}+\theta_{5}\theta_{1}\theta_{4}\rangle.$

This is an example of an algebra with two integer elements over one nilpotent
and one idempotent. This is linked to the spin Calogero systems whose relation
with bispectrality can be found in Bergvelt et al. (2009); Kasman (2015).

Theorems 1, 2 and 3 give positive answers to the Conjectures 1, 2 and 3 of
Grünbaum (2014) about three bispectral full rank 1 algebras. Moreover, these
algebras are Noetherian and finitely generated because they are contained in
the $N\times N$ matrix polynomial ring $M_{N}(\mathbb{K}[x])$ (See Vasquez and
Zubelli (2021)).

We close this section by remarking the important role played by Darboux
transformations, which goes back to Duistermaat and Grünbaum (1986) and
Zubelli and Magri (1991). Indeed, the three examples in this section are
instances of rational Darboux transformations from the scalar matrix
exponential functions. All such Darboux transformations were shown to be
bispectral in Theorems 1.1 and 1.2 of Geiger et al. (2017).

##  6 A few extensions of the problems discussed above

Since one of the goals of this paper is to serve as an “invitation” to look at
this problem extended to a wider audience, we give a road map with some
selected references. Indeed, the bispectral problem has many different
incarnations, and in our opinion we are still far from having a unified
theory.

In the scalar case here are some early papers that should be mentioned are
Kasman and Rothstein (1997); Wilson (1993, 1998); Bakalov et al. (1996, 1997,
1998a, 1998b); Grünbaum and Yakimov (2002); Geiger et al. (2017).

One natural issue concerns the numerical aspects involving the prolate
spheroidal functions. In this case, the reader may want to consult the book
Osipov et al. (2013).

Another direction concerns, the purely discrete (actually finite) version of
time-band-limiting. In an effort to better understand the commuting property
in question, one of us looked at the case when the real line is replaced by
the N roots of unity. See Grünbaum (1981) as well as Perlstadt (1984); Perline
(1987). Expression (11) in Grünbaum (1981) contains a small typo: the $r(2)$
in the denominator should be replaced by $r(1)$.

Moving on to the discrete-continuous version of the bispectral problem, we
have that for the scalar case, involving orthogonal polynomials satisfying
differential equations the problem had already been considered by S. Bochner
(and even earlier). A very good introduction to this is given in Haine (2001)
and its references. See also Grünbaum (2001); Grünbaum and Haine (1999);
Grünbaum et al. (1997, 1999).

For the matrix valued case there are two sources of early examples, one
resulting from the theory of matrix valued spherical functions, see Grünbaum
et al. (2001, 2002a, 2002b, 2003, 2004, 2005), and another one see Durán and
Grünbaum (2004, 2007); Durán and de la Iglesia (2008). See also Grünbaum and
Iliev (2003) as well as Boyallian and Liberati (2008); Miranian (2005); Reach
(1988).

Solutions of the bispectral problem can be used to obtain integral operators
which reflect some ordinary differential operator in the sense of (for
instance) Casper et al. (2019). This fact generalizes the commuting property
in the scalar case. It would be interesting to see whether this could be
extended to the matrix case.

In the (noncommutative) matrix case, by considering operators in the physical
and spectral variables acting from opposite directions, we maintain the Ad-
conditions that played a substantial role in Duistermaat and Grünbaum (1986).
In this case, this leads to the embedding of the bispectral algebras of
eigenvalues into the matrix polynomial algebra
$M_{N}(\mathbf{\mathbb{C}}[x])$. See Vasquez and Zubelli (2021).

Another natural direction would be to look for a characterization for algebras
relating to the spin-Calogero system for matrices of arbitrary size of matrix
$N$. The examples $N=1$ and $N=2$ were generalized to arbitrary matrix size
$N$ and was characterized as a subalgebra of $M_{N}(\mathbf{\mathbb{C}}[x])$
using a family of maps $\left\\{P_{k}\right\\}_{k\in\mathbf{\mathbb{N}}}$
satisfying some nice properties such as translation and a product rule similar
to the Leibniz rule. See Vasquez and Zubelli (2021).

Finally, let us go back to the original problem. We recall that both in the
scalar and in the matrix case the motivation behind the bispectral problem was
a desire to understand "what is behind" the remarkable commutativity property
between the operator of time-and-band limiting in the Fourier, Bessel and Airy
cases. It was hard to suspect that this problem would have connections with
many of the recent developments in integrable systems. There has been some
progress in connecting the bispectral problem with the commutativity property
mentioned above, and here again-at least in the scalar case- there are some
connections with Integrable systems, see for instance Grünbaum (2021); Casper
and Yakimov (2020); Casper et al. (2019, 2021a, 2021b). For connections
between the bispectral problem and the commuting property see Castro and
Grünbaum (2017); Grünbaum et al. (2017, 2018).

## Acknowledgments

BVDC was supported by CAPES grants 88882 332418/2019-01 as well as IMPA. JPZ
was supported by CNPq grants 302161 and 47408, as well as by FAPERJ under the
program Cientistas do Nosso Estado grant E-26/202.927/2017, Brazil. BDVC and
JPZ acknowledge the support from the FSU-2020-09 grant from Khalifa
University, UAE.

## References

  * Duistermaat and Grünbaum (1986) J.J. Duistermaat and F.A. Grünbaum.  Differential equations in the spectral parameter.  _Communications in Mathematical Physics_ , 103(2):177–240, 1986.
  * Slepian and Pollak (1961) D. Slepian and H. O. Pollak.  Prolate spheroidal wave functions, Fourier analysis and uncertainty. I.  _Bell Syst. Tech. J._ , 40:43–63, 1961.  ISSN 0005-8580.  doi:10.1002/j.1538-7305.1961.tb03976.x.
  * Landau and Pollak (1962) H. J. Landau and H. O. Pollak.  Prolate spheroidal wave functions, Fourier analysis and uncertainty, III: The dimension of the space of essentially time- and band-limited signals.  _Bell Syst. Tech. J._ , 41:1295–1336, 1962.  ISSN 0005-8580.  doi:10.1002/j.1538-7305.1962.tb03279.x.
  * Landau and Pollak (1961) H. J. Landau and H. O. Pollak.  Prolate spheroidal wave functions, Fourier analysis and uncertainty. II.  _Bell Syst. Tech. J._ , 40:65–84, 1961.  ISSN 0005-8580.  doi:10.1002/j.1538-7305.1961.tb03977.x.
  * Slepian (1964) D. Slepian.  Prolate spheroidal wave functions, Fourier analysis and uncertainty. IV: Extensions to many dimensions; generalized prolate spheroidal functions.  _Bell Syst. Tech. J._ , 43:3009–3057, 1964.  ISSN 0005-8580.  doi:10.1002/j.1538-7305.1964.tb01037.x.
  * Slepian (1978) D. Slepian.  Prolate spheroidal wave functions, Fourier analysis, and uncertainty - V: The discrete case.  _Bell Syst. Tech. J._ , 57:1371–1430, 1978.  ISSN 0005-8580.  doi:10.1002/j.1538-7305.1978.tb02104.x.
  * Slepian (1983) David Slepian.  Some comments on Fourier analysis, uncertainty and modeling.  _SIAM Review_ , 25(3):379–393, 1983.  ISSN 00361445.  URL http://www.jstor.org/stable/2029386.
  * Slepian (1976) David S. Slepian.  On bandwidth.  _Proceedings of the IEEE_ , 64:292–300, 1976.
  * Connes and Consani (2021) Alain Connes and Caterina Consani.  Spectral triples and zeta-cycles.  _arXiv preprint:2106.01715_ , 2021.
  * Connes and Marcolli (2008) Alain Connes and Matilde Marcolli.  _Noncommutative geometry, quantum fields and motives_ , volume 47.  AMS Colloquium Publications, Volume 55, New Delhi: Hindustan Book Agency, 2008.  ISBN 978-81-85931-80-7.
  * Connes and Moscovici (2021) Alain Connes and Henri Moscovici.  Prolate spheroidal operator and zeta.  _arXiv preprint:2112.05500_ , 2021.
  * Connes (2021) Alain Connes.  Prolate spheroidal functions and zeta, 2021.  URL https://www.youtube.com/watch?v=d8LgqGApvlo.
  * Grünbaum (2014) F. Alberto Grünbaum.  Some noncommutative matrix algebras arising in the bispectral problem.  _SIGMA, Symmetry Integrability Geom. Methods Appl._ , 10:paper 078, 9, 2014.  ISSN 1815-0659.  doi:10.3842/SIGMA.2014.078.
  * Tracy and Widom (1994a) Craig A. Tracy and Harold Widom.  Fredholm determinants, differential equations and matrix models.  _Commun. Math. Phys._ , 163(1):33–72, 1994a.  ISSN 0010-3616.  doi:10.1007/BF02101734.
  * Tracy and Widom (1994b) Craig A. Tracy and Harold Widom.  Level-spacing distributions and the Airy kernel.  _Commun. Math. Phys._ , 159(1):151–174, 1994b.  ISSN 0010-3616.  doi:10.1007/BF02100489.
  * Grünbaum et al. (1982) F. A. Grünbaum, L. Longhi, and M. Perlstadt.  Differential operators commuting with finite convolution integral operators: Some non-Abelian examples.  _SIAM J. Appl. Math._ , 42:941–955, 1982.  ISSN 0036-1399.  doi:10.1137/0142067.
  * Zubelli and Magri (1991) Jorge P. Zubelli and Franco Magri.  Differential equations in the spectral parameter, Darboux transformations and a hierarchy of master symmetries for KdV.  _Commun. Math. Phys._ , 141(2):329–351, 1991\.  ISSN 0010-3616.  doi:10.1007/BF02101509.
  * Zubelli and Valerio Silva (2000) Jorge P. Zubelli and D. S. Valerio Silva.  Rational solutions of the master symmetries of the KdV equation.  _Commun. Math. Phys._ , 211(1):85–109, 2000\.  ISSN 0010-3616.  doi:10.1007/s002200050803.
  * Wilson (1993) George Wilson.  Bispectral commutative ordinary differential operators.  _J. Reine Angew. Math._ , 442:177–204, 1993.  ISSN 0075-4102.  doi:10.1515/crll.1993.442.177.
  * Zubelli (1990) J.P. Zubelli.  Differential equations in the spectral parameter for matrix differential operators.  _Physica D: Nonlinear Phenomena_ , 43(2-3):269–287, 1990.
  * Zubelli (1992) Jorge P. Zubelli.  On a zero curvature problem related to the ZS-AKNS operator.  _J. Math. Phys._ , 33(11):3666–3675, 1992.  ISSN 0022-2488.
  * Zubelli (1992) Jorge P. Zubelli.  Rational solutions of nonlinear evolution equations, vertex operators, and bispectrality.  _J. Differ. Equations_ , 97(1):71–98, 1992\.  ISSN 0022-0396.  doi:10.1016/0022-0396(92)90084-Z.
  * Sakhnovich and Zubelli (2001) Alexander Sakhnovich and Jorge P. Zubelli.  Bundle bispectrality for matrix differential equations.  _Integral Equations Oper. Theory_ , 41(4):472–496, 2001.  ISSN 0378-620X.  doi:10.1007/BF01202105.
  * Kasman (2015) Alex Kasman.  Bispectrality of $N$-component KP wave functions: a study in non-commutativity.  _SIGMA, Symmetry Integrability Geom. Methods Appl._ , 11:paper 087, 22, 2015.  ISSN 1815-0659.  doi:10.3842/SIGMA.2015.087.
  * Boyallian and Liberati (2008) Carina Boyallian and Jose I. Liberati.  Matrix-valued bispectral operators and quasideterminants.  _J. Phys. A, Math. Theor._ , 41(36):11, 2008\.  ISSN 1751-8113.  doi:10.1088/1751-8113/41/36/365209.  Id/No 365209.
  * Grünbaum and Iliev (2003) F. Alberto Grünbaum and Plamen Iliev.  A noncommutative version of the bispectral problem.  _J. Comput. Appl. Math._ , 161(1):99–118, 2003\.  ISSN 0377-0427.  doi:10.1016/S0377-0427(03)00548-X.
  * Geiger et al. (2017) Joel Geiger, Emil Horozov, and Milen Yakimov.  Noncommutative bispectral Darboux transformations.  _Trans. Am. Math. Soc._ , 369(8):5889–5919, 2017.  ISSN 0002-9947.  doi:10.1090/tran/6950.
  * Ablowitz et al. (1974) Mark J. Ablowitz, David J. Kaup, Alan C. Newell, and Harvey Segur.  The inverse scattering transform-Fourier analysis for nonlinear problems.  _Stud. Appl. Math._ , 53:249–315, 1974.  ISSN 0022-2526.  doi:10.1002/sapm1974534249.
  * Berest and Veselov (1994) Yu Yu Berest and A P Veselov.  Huygens’ principle and integrability.  _Russian Mathematical Surveys_ , 49(6):5, 1994\.  URL http://stacks.iop.org/0036-0279/49/i=6/a=R02.
  * Chalub and Zubelli (2006) Fabio A. C. C. Chalub and Jorge P. Zubelli.  Huygens’ principle for hyperbolic operators and integrable hierarchies.  _Phys. D_ , 213(2):231–245, 2006.  ISSN 0167-2789.  doi:10.1016/j.physd.2005.11.008.  URL http://dx.doi.org/10.1016/j.physd.2005.11.008.
  * Chalub and Zubelli (2005) Fabio A. C. C. Chalub and Jorge P. Zubelli.  Huygens’ principle, Dirac operators, and rational solutions of the AKNS hierarchy.  _Math. Phys. Anal. Geom._ , 8(3):187–210, 2005\.  ISSN 1385-0172.  doi:10.1007/s11040-004-6683-3.  URL http://dx.doi.org/10.1007/s11040-004-6683-3.
  * Chalub and Zubelli (2004) Fabio A. C. C. Chalub and Jorge P. Zubelli.  Matrix bispectrality and Huygens’ principle for Dirac operators.  In _Partial differential equations and inverse problems_ , volume 362 of _Contemp. Math._ , pages 89–112. Amer. Math. Soc., Providence, RI, 2004.  doi:10.1090/conm/362/06607.  URL http://dx.doi.org/10.1090/conm/362/06607.
  * Chalub and Zubelli (2001) Fabio A. C. C. Chalub and Jorge P. Zubelli.  On Huygens’ principle for Dirac operators and nonlinear evolution equations.  _J. Nonlinear Math. Phys._ , 8(suppl.):62–68, 2001.  ISSN 1402-9251.  Nonlinear evolution equations and dynamical systems (Kolimbary, 1999).
  * Castro and Grünbaum (2006) Mirta M. Castro and F. Alberto Grünbaum.  The algebra of differential operators associated to a family of matrix-valued orthogonal polynomials: five instructive examples.  _Int. Math. Res. Not._ , 2006(7):33, 2006.  ISSN 1073-7928.  doi:10.1155/IMRN/2006/47602.  Id/No 47602.
  * Tirao (2011) Juan Tirao.  The algebra of differential operators associated to a weight matrix: a first example.  In _Groups, algebras and applications. XVIII Latin American algebra colloquium, São Pedro, Brazil, August 3–8, 2009. Proceedings._ , pages 291–324. Providence, RI: American Mathematical Society (AMS), 2011.  ISBN 978-0-8218-5239-2.
  * Grünbaum and Tirao (2007) F. Alberto Grünbaum and Juan Tirao.  The algebra of differential operators associated to a weight matrix.  _Integral Equations Oper. Theory_ , 58(4):449–475, 2007.  ISSN 0378-620X.  doi:10.1007/s00020-007-1517-x.
  * Casper (2018) W. Riley Casper.  Elementary examples of solutions to Bochner’s problem for matrix differential operators.  _J. Approx. Theory_ , 229:36–71, 2018.  ISSN 0021-9045.  doi:10.1016/j.jat.2018.02.004.
  * Zurrián (2017) Ignacio Zurrián.  The algebra of differential operators for a matrix weight: an ultraspherical example.  _Int. Math. Res. Not._ , 2017(8):2402–2430, 2017.  ISSN 1073-7928.  doi:10.1093/imrn/rnw104.
  * Vasquez and Zubelli (2021) Brian D. Vasquez and Jorge P. Zubelli.  Matrix bispectrality of full rank one algebras.  _arXiv preprint:2108.06997_ , 2021.
  * Bergvelt et al. (2009) Maarten Bergvelt, Michael Gekhtman, and Alex Kasman.  Spin Calogero particles and bispectral solutions of the matrix KP hierarchy.  _Math. Phys. Anal. Geom._ , 12(2):181–200, 2009\.  ISSN 1385-0172.  doi:10.1007/s11040-009-9058-y.
  * Zubelli and Magri (1991) Jorge P. Zubelli and Franco Magri.  Differential equations in the spectral parameter, Darboux transformations and a hierarchy of master symmetries for KdV.  _Comm. Math. Phys._ , 141(2):329–351, 1991.  ISSN 0010-3616.
  * Kasman and Rothstein (1997) Alex Kasman and Mitchell Rothstein.  Bispectral Darboux transformations: The generalized Airy case.  _Physica D_ , 102(3-4):159–176, 1997.  ISSN 0167-2789.  doi:10.1016/S0167-2789(96)00208-4.
  * Wilson (1998) George Wilson.  Collisions of Calogero-Moser particles and an adelic Grassmannian (with an appendix by I. G. Macdonald).  _Invent. Math._ , 133(1):1–41, 1998.  ISSN 0020-9910.  doi:10.1007/s002220050237.
  * Bakalov et al. (1996) B. Bakalov, E. Horozov, and M. Yakimov.  General methods for constructing bispectral operators.  _Phys. Lett., A_ , 222(1-2):59–66, 1996.  ISSN 0375-9601.  doi:10.1016/0375-9601(96)00624-X.
  * Bakalov et al. (1997) B. Bakalov, E. Horozov, and M. Yakimov.  Bispectral algebras of commuting ordinary differential operators.  _Commun. Math. Phys._ , 190(2):331–373, 1997\.  ISSN 0010-3616.  doi:10.1007/s002200050244.
  * Bakalov et al. (1998a) B. Bakalov, E. Horozov, and M. Yakimov.  Highest weight modules over the $W_{1+\infty}$ algebra and the bispectral problem.  _Duke Math. J._ , 93(1):41–72, 1998a.  ISSN 0012-7094.  doi:10.1215/S0012-7094-98-09302-4.
  * Bakalov et al. (1998b) Bojko Bakalov, Emil Horozov, and Milen Yakimov.  Automorphisms of the Weyl algebra and bispectral operators.  In _The bispectral problem. Proceedings of the CRM workshop, Montréal, Canada, March 1997_ , pages 3–10. Providence, RI: American Mathematical Society, 1998b.  ISBN 0-8218-0949-0.
  * Grünbaum and Yakimov (2002) F. Alberto Grünbaum and Milen Yakimov.  Discrete bispectral Darboux transformations from Jacobi operators.  _Pac. J. Math._ , 204(2):395–431, 2002.  ISSN 1945-5844.  doi:10.2140/pjm.2002.204.395.
  * Osipov et al. (2013) Andrei Osipov, Vladimir Rokhlin, and Hong Xiao.  _Prolate spheroidal wave functions of order zero. Mathematical tools for bandlimited approximation_ , volume 187.  New York, NY: Springer, 2013.  ISBN 978-1-4614-8258-1; 978-1-4614-8259-8.  doi:10.1007/978-1-4614-8259-8.
  * Grünbaum (1981) F. Alberto Grünbaum.  Eigenvectors of a Toeplitz matrix: Discrete version of the prolate spheroidal wave functions.  _SIAM J. Algebraic Discrete Methods_ , 2:136–141, 1981\.  ISSN 0196-5212.  doi:10.1137/0602017.
  * Perlstadt (1984) Marci Perlstadt.  A property of orthogonal polynomial families with polynomial duals.  _SIAM J. Math. Anal._ , 15:1043–1054, 1984.  ISSN 0036-1410.  doi:10.1137/0515081.
  * Perline (1987) Ronald Keith Perline.  Discrete time-band limiting operators and commuting tridiagonal matrices.  _SIAM J. Algebraic Discrete Methods_ , 8:192–195, 1987\.  ISSN 0196-5212.  doi:10.1137/0608016.
  * Haine (2001) Luc Haine.  The Bochner-Krall problem: Some new perspectives.  In _Special functions 2000: current perspective and future directions. Proceedings of the NATO Advanced Study Institute, Tempe, AZ, USA, May 29–June 9, 2000_ , pages 141–178. Dordrecht: Kluwer Academic Publishers, 2001\.  ISBN 0-7923-7119-4.
  * Grünbaum (2001) F. Alberto Grünbaum.  The bispectral problem: An overview.  In _Special functions 2000: current perspective and future directions. Proceedings of the NATO Advanced Study Institute, Tempe, AZ, USA, May 29–June 9, 2000_ , pages 129–140. Dordrecht: Kluwer Academic Publishers, 2001\.  ISBN 0-7923-7119-4.
  * Grünbaum and Haine (1999) F. Alberto Grünbaum and Luc Haine.  Associated polynomials, spectral matrices and the bispectral problem.  _Methods Appl. Anal._ , 6(2):209–224, 1999\.  ISSN 1073-2772.  doi:10.4310/MAA.1999.v6.n2.a6.
  * Grünbaum et al. (1997) F. Alberto Grünbaum, Luc Haine, and Emil Horozov.  On the Krall-Hermite and the Krall-Bessel polynomials.  _Int. Math. Res. Not._ , 1997(19):953–966, 1997\.  ISSN 1073-7928.  doi:10.1155/S1073792897000615.
  * Grünbaum et al. (1999) F. Alberto Grünbaum, Luc Haine, and Emil Horozov.  Some functions that generalize the Krall-Laguerre polynomials.  _J. Comput. Appl. Math._ , 106(2):271–297, 1999\.  ISSN 0377-0427.
  * Grünbaum et al. (2001) F. A. Grünbaum, I. Pacharoni, and J. Tirao.  A matrix-valued solution to Bochner’s problem.  _J. Phys. A, Math. Gen._ , 34(48):10647–10656, 2001.  ISSN 0305-4470.  doi:10.1088/0305-4470/34/48/328.
  * Grünbaum et al. (2002a) F. A. Grünbaum, I. Pacharoni, and J. Tirao.  Matrix valued spherical functions associated to the complex projective plane.  _J. Funct. Anal._ , 188(2):350–441, 2002a.  ISSN 0022-1236.  doi:10.1006/jfan.2001.3840.
  * Grünbaum et al. (2002b) F. A. Grünbaum, I. Pacharoni, and J. Tirao.  Matrix valued spherical functions associated to the three dimensional hyperbolic space.  _Int. J. Math._ , 13(7):727–784, 2002b.  ISSN 0129-167X.  doi:10.1142/S0129167X02001460.
  * Grünbaum et al. (2003) F. A. Grünbaum, I. Pacharoni, and J. Tirao.  Matrix valued orthogonal polynomials of the Jacobi type.  _Indag. Math., New Ser._ , 14(3-4):353–366, 2003.  ISSN 0019-3577.  doi:10.1016/S0019-3577(03)90051-6.
  * Grünbaum et al. (2004) F. Alberto Grünbaum, Inés Pacharoni, and Juan Tirao.  An invitation to matrix-valued spherical functions: Linearization of products in the case of complex projective space $P_{2}(\mathbb{C})$.  In _Modern signal processing._ , pages 147–160. Cambridge: Cambridge University Press, 2004.  ISBN 0-521-82706-X.
  * Grünbaum et al. (2005) F. Alberto Grünbaum, Inés Pacharoni, and Juan Alfredo Tirao.  Matrix valued orthogonal polynomials of Jacobi type: the role of group representation theory.  _Annales de l’Institut Fourier_ , 55(6):2051–2068, 2005.  doi:10.5802/aif.2151.  URL http://www.numdam.org/articles/10.5802/aif.2151/.
  * Durán and Grünbaum (2004) Antonio J. Durán and F. Alberto Grünbaum.  Orthogonal matrix polynomials satisfying second-order differential equations.  _Int. Math. Res. Not._ , 2004(10):461–484, 2004\.  ISSN 1073-7928.  doi:10.1155/S1073792804132583.
  * Durán and Grünbaum (2007) Antonio J. Durán and F. Alberto Grünbaum.  Matrix orthogonal polynomials satisfying second-order differential equations: coping without help from group representation theory.  _J. Approx. Theory_ , 148(1):35–48, 2007.  ISSN 0021-9045.  doi:10.1016/j.jat.2007.02.003.
  * Durán and de la Iglesia (2008) Antonio J. Durán and Manuel D. de la Iglesia.  Second-order differential operators having several families of orthogonal matrix polynomials as eigenfunctions.  _Int. Math. Res. Not._ , 2008:24, 2008.  ISSN 1073-7928.  doi:10.1093/imrn/rnn084.  Id/No rnn084.
  * Miranian (2005) L. Miranian.  Matrix-valued orthogonal polynomials on the real line: some extensions of the classical theory.  _J. Phys. A, Math. Gen._ , 38(25):5731–5749, 2005.  ISSN 0305-4470.  doi:10.1088/0305-4470/38/25/009.
  * Reach (1988) Michael Reach.  Generating difference equations with the Darboux transformation.  _Communications in Mathematical Physics_ , 119(3):385–402, September 1988.  doi:10.1007/BF01218079.
  * Casper et al. (2019) W. Riley Casper, F. Alberto Grünbaum, Milen Yakimov, and Ignacio Zurrián.  Reflective prolate-spheroidal operators and the KP/KdV equations.  _Proc. Natl. Acad. Sci. USA_ , 116(37):18310–18315, 2019.  ISSN 0027-8424.  doi:10.1073/pnas.1906098116.
  * Grünbaum (2021) F. Alberto Grünbaum.  Commuting integral and differential operators and the master symmetries of the Korteweg-de Vries equation.  _Inverse Probl._ , 37(8):21, 2021.  ISSN 0266-5611.  doi:10.1088/1361-6420/ac11c6.  Id/No 085010.
  * Casper and Yakimov (2020) W. Riley Casper and Milen T. Yakimov.  Integral operators, bispectrality and growth of Fourier algebras.  _J. Reine Angew. Math._ , 766:151–194, 2020.  ISSN 0075-4102.  doi:10.1515/crelle-2019-0031.
  * Casper et al. (2021a) W. Riley Casper, F. Alberto Grünbaum, Milen Yakimov, and Ignacio Zurrián.  Reflective prolate-spheroidal operators and the adelic grassmannian.  _to appear in Comm. Pure and Applied Mathematics_ , 2021a.
  * Casper et al. (2021b) W. R. Casper, F. A. Grünbaum, M. Yakimov, and I. Zurrián.  Algebras of differential operators for kernels of Airy type.  _in preparation_ , 2021b.
  * Castro and Grünbaum (2017) M. Castro and F. A. Grünbaum.  Time-and-band limiting for matrix orthogonal polynomials of Jacobi type.  _Random Matrices Theory Appl._ , 6(4):12, 2017\.  ISSN 2010-3263.  doi:10.1142/S2010326317400019.  Id/No 1740001.
  * Grünbaum et al. (2017) F. A. Grünbaum, I. Pacharoni, and I. Zurrián.  Time and band limiting for matrix valued functions: an integral and a commuting differential operator.  _Inverse Probl._ , 33(2):14, 2017.  ISSN 0266-5611.  doi:10.1088/1361-6420/aa53b8.  Id/No 025005.
  * Grünbaum et al. (2018) F. Alberto Grünbaum, Luc Vinet, and Alexei Zhedanov.  Algebraic Heun operator and band-time limiting.  _Commun. Math. Phys._ , 364(3):1041–1068, 2018\.  ISSN 0010-3616.  doi:10.1007/s00220-018-3190-0.