# A Kitaev–type spin liquid on a quasicrystal

M. A. Keskiner<EMAIL_ADDRESS>Department of Physics, Bilkent
University, Ankara, 06800, TURKEY  O. Erten<EMAIL_ADDRESS>Department of
Physics, Arizona State University, Tempe, AZ 85287, USA  M. Ö. Oktel
<EMAIL_ADDRESS>Department of Physics, Bilkent University, Ankara,
06800, TURKEY

(August 28, 2024)

###### Abstract

We develop an exactly solvable model with Kitaev-type interactions and study
its phase diagram on the dual lattice of the quasicrystalline Ammann-Beenker
lattice. Our construction is based on the $\Gamma$-matrix generalization of
the Kitaev model and utilizes the cut-and-project correspondence between the
four-dimensional simple cubic lattice and the Ammann-Beenker lattice to
designate four types of bonds. We obtain a rich phase diagram with gapped
(chiral and abelian) and gapless spin liquid phases via Monte Carlo
simulations and variational analysis. We show that the ground state can be
further tuned by the inclusion of an onsite term that selects 21 different
vison configurations while maintaining the integrability of the model. Our
results highlight the rich physics at the intersection of quasicrystals and
quantum magnetism.

##  I Introduction

Quantum spin liquids (QSLs) are the disordered phases of magnetic systems and
exhibit exceptional properties such as fractionalized excitations and long-
range entanglement due to their underlying topological order[1, 2, 3, 4, 5, 6,
7, 8]. The Kitaev model on the honeycomb lattice is a foundational model as it
is the first exactly solvable spin model that showcases a spin liquid ground
state with both gapless and gapped phases, featuring abelian and non-abelian
anyonic excitations[9]. Recently, there has been a tour de force effort in
discovering new materials with strong Kitaev interactions such as iridates and
$\alpha-$RuCl3[10, 11, 12]. Kitaev interactions may be strong in other van der
Waals (vdW) materials like CrI3 [13].

Most models of quantum magnetism, particularly models exhibiting a QSL ground
state, have been explored in the context of periodic systems. While QSLs lack
long-range magnetic order, they are in general defined on models with perfect
translational symmetry. For instance, recent works have explored the
consequences of translational symmetry breaking for the QSL state, either by
considering defects in a periodic lattice[14, 15, 16, 17] or by defining
randomly generated amorphous lattices[18, 19]. Solids do not have to be either
periodic structures or random glasses. A third possibility is
quasiperiodicity.

Quasicrystals are unique classes of materials that exhibit a regular atomic
structure with non-repeating patterns, forming a contrast between the
disordered arrangement of glasses and the fully periodic arrangement of
crystals. They can be characterized by the presence of long-range order in
terms of their translational and orientational symmetries, but the patterns do
not repeat at regular intervals[20, 21, 22]. Due to this unique arrangement of
the sites, they exhibit rare features such as strictly-localized states [23,
24], or electronic states which are neither localized nor extended [25, 26].
Quasicrystals can also exhibit different symmetries, including five- and
eight-fold rotations, forbidden in conventional crystal structures. These
forbidden symmetries lead to unprecedented topological states that are not
allowed in crystalline materials[27, 28, 29].

Figure 1: (a) Ammann-Beenker lattice (ABL) and the coloring scheme consisting
of four type of bonds obtained via projecting from the hyperspace. (b) The
dual lattice of ABL (dABL) and the corresponding coloring scheme. The pink and
blue dots show the position of the lattice sites in the ABL and dABL
respectively.

The understanding of the elementary excitation spectrum of quasicrystals is
quite limited compared to periodic solids. Experimentally, most quasicrystals
are found to be poor conductors [30]. However, non-Fermi liquid behavior[31]
and superconductivity [32] have also been observed. Recent experiments with
high-quality samples have demonstrated that ferromagnetic long-range order is
possible in a quasicrystal [33]. The experimental realization of
quasicrystalline systems goes beyond the synthesis of alloys. Vapor deposition
on metallic surfaces[34], as well as atomic force microscopy assembly of
molecules [35] have created large two-dimensional quasicrystals. Synthetic
quasicrystals for light [36], and ultracold atoms [37] promise to probe
bosonic models in quasicrystals. In a more recent development, quasicrystals
can also be formed in moiré superlattices of vdW materials[38]. Experimental
possibilities for testing the excitation spectrum of quasicrystals are rapidly
increasing.

Classical models of magnetism, such as the Ising model, have been studied in
some detail for quasicrystals [39, 40]. The non-trivial structural properties
of quasicrystals force a reconsideration of the basic models, such as dimer
coverings on the simplest two-dimensional lattices [41, 42]. However, quantum
magnetism research in quasicrystals is at its early stages. Numerical work on
small systems [43] and two-dimensional lattices [44, 45, 46] indicate elements
of frustration and self-similar magnetic ground states. Recently, Ref. 47
showed that exact dimer wave functions can be constructed in certain
quasicrystals for generalized Heisenberg Hamiltonians.

Here, we investigate the interplay between quasiperiodicity and spin liquid
order by constructing an exactly solvable model with Kitaev-type interactions
on a quasicrystal. The Kitaev model and its generalizations are integrable due
to the presence of conserved quantities at each mesh of the lattice. The
conservation of such mesh loops requires two structural properties. First, the
lattice must have the same coordination number $z$ on each site. Second, the
link lattice itself must be $z-$partite, i.e., the links of the lattice can be
labeled by one of the $z$ colors so that all links starting from any site have
different colors. The two most common quasicrystal models, the Penrose [48]
(PL) and Ammann-Beenker [49, 50] (ABL) lattices, have varying coordination
numbers for their vertices. However, both of these models have quadrilateral
tiles and their dual lattices have a constant coordination number $z=4$. The
tight-binding models on the dual lattices are also referred to as the center
tight-binding models [51]. Although the condition of uniform coordination is
easily satisfied by considering the dual lattices, the second condition for
the z-partite coloring of the lattice is non-trivial. For the dual of the
Ammann-Beenker lattice (dABL), we find a z-partite coloring rule by lifting
the ABL back into the four-dimensional space. We show that the four-
dimensional cubic lattice with nearest neighbor links can be colored in a way
so that each square has sides of 4 different colors. This rule allows us to
define the Kitaev-type Hamiltonian on the dABL. It is worth remarking that the
same method cannot be extended to a five-dimensional simple cubic lattice, and
we have not been able to define a similar model on the PL.

Our coloring scheme for the ABL obtained by projection and corresponding dABL
are shown in Fig. 1(a), and Fig. 1(b). In this model, an exact solution of the
spin model can be achieved via a Majorana fermion representation of the
$\Gamma$ matrices. dABL consists of plaquettes with both even- and odd-
numbered sides. The odd-numbered plaquettes break time reversal symmetry[9]
and allow for chiral spin liquid phases that are classified by a Chern number,
$\nu$. We obtain the ground state phase diagram via Monte Carlo simulations
and variational analysis. We find that the ground state follows a flux
configuration that is an extension of the Lieb’s theorem [52], even though the
theorem was not proved for non-bipartite lattices. Depending on the coupling
constants ($J_{\mu},\leavevmode\nobreak\ \mu=1,2,3,4$), we find that the
ground state can be gapped with $\nu=2$ or $\nu=0$. These phases are separated
by a gapless phase, as shown in Fig. 4. Upon inclusion of an onsite term that
commutes with the flux operators, the ground state vison configuration changes
as shown in Fig. 6. We find that up to 21 different flux configurations can be
stabilized as a function of the onsite field strength.

The paper is organized as follows. Section II briefly overviews the Ammann-
Beenker lattice and its periodic approximants. Section III focuses on
constructing our model with Kitaev-type interactions, including our coloring
scheme. In section IV, we present our results and discuss the phase diagram.
We conclude with a discussion and a summary of our results.

##  II Ammann-Beenker Lattice

We use one of the most well-known quasicrystals, the Ammann-Beenker Lattice
(ABL), to construct our model. All the points forming the ABL can be written
as

$\vec{R}_{ABL}=k_{0}\hat{e}_{0}+k_{1}\hat{e}_{1}+k_{2}\hat{e}_{2}+k_{3}\hat{e}_{3},$
(1)

where $k_{j}$ are integers and the star-vectors are

$\hat{e}_{0}=\hat{x},\,\hat{e}_{1}=\frac{1}{\sqrt{2}}\left(\hat{x}+\hat{y}\right),\,\hat{e}_{2}=\hat{y},\,\hat{e}_{3}=\frac{1}{\sqrt{2}}\left(-\hat{x}+\hat{y}\right).$
(2)

All the bonds in the ABL are parallel to one of the above star vectors, making
the meshes of ABL either squares or $\pi/4$ rhombi.

If the integers $k_{j}$ in Eq. 1 are allowed to vary independently, the set of
lattice points will uniformly fill the plane and create points that are
infinitesimally close to each other. The definition of the quasicrystal
lattice can be seen as a way of constraining the set
$k_{0},k_{1},k_{2},k_{3}$. Although alternative methods exist for defining
this constraint, we use the cut-and-project method which relates the ABL to
the four-dimensional cubic lattice. Each quadruple $k_{0},k_{1},k_{2},k_{3}$,
defines a unit tessaract. ABL is formed by using $k_{j}$ which has an
intersection with a specific two-plane. The two-plane is chosen so that the
projection of unit vectors in four dimensions onto this plane defines the star
vectors given above [50].

The first condition necessary for an exactly solvable Kitaev model is to
obtain a lattice with a uniform coordination number. ABL has vertices with
coordination numbers varying from 3 to 8, hence does not satisfy this
condition. However, as each mesh of ABL is a quadrilateral, its dual lattice
(dABL) is uniformly coordinated. dABL is obtained by placing vertices to the
center of each mesh and forming a link between two sites if their meshes share
an edge. So each dABL site is expressed as
$\vec{R}_{dABL}=\vec{R}_{ABL}\pm\frac{\hat{e}_{m}}{2}\pm\frac{\hat{e}_{n}}{2}$,
where the edges of the mesh including this vertex are parallel to
$\hat{e}_{m}$ and $\hat{e}_{n}$. Tight binding models have been studied on
dABL and other dual lattices, where such models are referred to as center
models [51].

The second condition we need to define a Kitaev-type model is to make sure
that the uniformly coordinated dABL is 4-partite. Each of the links in the
dABL has to be assigned an index from 1 to 4, with the condition that all 4
links meeting at every lattice site has a different index. The same condition
can be expressed in terms of the direct ABL as coloring all the links so that
each mesh has 4 edges of different colors. By considering the four-dimensional
simple cubic lattice, we establish a coloring rule that satisfies this
condition. As the cut-and-project construction shows every mesh in the ABL
corresponds to a mesh in a two-plane of the four dimensional cubic lattice.
The details of this coloring scheme for the four-dimensional lattice which
satisfies the coloring condition is presented in the Appendix A. The resulting
coloring for the ABL and dABL is shown in Fig. 1. Our coloring rule cannot be
trivially generalized to higher dimensions, hence to other quasicrystals such
as the PL.

While our model is defined for the infinite quasicrystal, numerical
calculations require a finite region. We use two kinds of regions for our
calculations. First, we use a finite region with open boundary conditions. We
choose different local configurations for the finite region to make sure that
our results are not specific to a particular patch of the quasicrystal. We use
a periodic approximant for the ABL [53] as the second method. Periodic
approximants are obtained by choosing a cut plane with a rational slope in the
cut-and-project construction. The approximants are indexed by an integer $s$
which increases the number of sites within the unit cell and creates the
infinite ABL in the limit $s\rightarrow\infty$. For $s=1,2,3$ the periodic
unit cell has $N=7,41,239$ sites respectively. The unit cell for the ABL
approximants is a square, and we double the unit cell in both directions for
compatibility with our coloring rules. In Fig. 2 we display the enlarged unit
cell for the $s=2$ approximant with the dABL coloring.

Figure 2: Four unit cells of the approximant, represented by s=2 for dABL. The
black dotted lines show the original ABL and the purple dashed lines show a
single unit cell of the s=2 approximant. Figure 3: Ground state flux
configuration for $J_{5}=0$ for a s=3 approximant. The $0$-flux, $\pi$-flux,
$\frac{\pi}{2}$-flux, and $-\frac{\pi}{2}$-flux are represented by yellow,
blue, red, and green colors, respectively. The ground state breaks time-
reversal symmetry due to odd-numbered plaquettes. Figure 4: Ground state phase
diagram for $J_{1}+J_{2}+J_{3}+J_{4}=1$ and $J_{5}=0$. (a), (b), (c), (d)
panels have $J_{4}$=0, 0.25, 0.5, and 0.75 respectively. The gap value is
normalized by the bandwidth which varies depending on the parameters. The
phase diagram contains two gapped phases with $\nu=2$ and $\nu=0$ which are
separated by a gapless region.

##  III Microscopic Model

An essential requirement for solving the Kitaev model exactly is the
anticommutation relations of the Pauli matrices,
$\\{\sigma_{i},\sigma_{j}\\}=2\delta_{ij}$. Since there are only three Pauli
matrices, this approach is applicable only to lattices with a coordination
number of $z=3$, such as the honeycomb, hyperhoneycomb, and hyperoctagon
lattices. However, an extension of Kitaev’s method is possible using $\Gamma$
matrices that follow the Clifford algebra
$\\{\Gamma_{i},\Gamma_{j}\\}=2\delta_{ij}$ [54, 55]. For instance, in a four-
dimensional representation of the Clifford algebra, there exist five
$\Gamma^{\mu}$ operators, along with ten
$\Gamma^{\mu\nu}=\frac{i}{2}[\Gamma^{\mu},\Gamma^{\nu}]$ operators and an
identity matrix, which form the basis of the local Hilbert space.
Consequently, Kitaev’s construction can be expanded to lattices with a
coordination number of up to $z=5$ [54, 55]. dABL has coordination number
$z=4$, represented by four different colors for each type of bond. Therefore,
we use this representation and consider the Hamiltonian,

$\displaystyle H_{K}=-\sum_{\langle
ij\rangle_{\mu}}J_{\mu}(\Gamma_{i}^{\mu}\Gamma_{j}^{\mu}+\Gamma_{i}^{\mu
5}\Gamma_{j}^{\mu 5})+J_{5}\sum_{j}\Gamma_{j}^{5}$ (3)

where $\mu$ is the type if the bond, $\mu=\\{1,2,3,4\\}$. Since the
$\Gamma^{5}$ operator is not used as a bond operator, it can be included as an
onsite term. Similar models have been considered for other $z=4$ lattices such
as square lattice[56, 57, 55]. Note that four-dimensional $\Gamma$ matrices
can be represented by two sets of Pauli matrices [58, 59, 60] with Kugel-
Khomskii interactions[61] or $J=3/2$ operators[55]. For the latter,
$\Gamma^{\mu}$ operators correspond to quadrupolar operators. For each
plaquette, we define $W_{p}=\prod_{(j,k)\in
p}\Gamma_{j}^{\mu}\Gamma_{k}^{\mu}$ where the product is taken in the counter-
clockwise direction. These plaquette operators commute with the Hamiltonian,
$[H,W_{p}]=0$ and each other $[W_{p},W_{p^{\prime}}]=0$. Therefore this leads
to infinitely many conserved quantities for the model. An exact solution can
be obtained by a Majorana representation of the $\Gamma$ matrices [54],

$\Gamma_{j}^{\mu}=ib_{j}^{\mu}c_{j},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\Gamma_{j}^{\mu\nu}=ib_{j}^{\mu}b_{j}^{\nu}$ (4)

where we introduce a total of 6 Majorana fermions per site. Relabeling
$b_{j}^{5}$ as $b_{j}^{5}\rightarrow d_{j}$, the Hamiltonian can be
reexpressed as

$\displaystyle H=\sum_{\langle
ij\rangle_{\mu}}J_{\mu}iu_{ij}^{\mu}(c_{i}c_{j}+d_{i}d_{j})+J_{5}\sum_{i}id_{i}c_{i}$
(5)

where $u_{ij}^{\mu}=ib_{i}^{\mu}b_{j}^{\mu}$. $u_{ij}^{\mu}$ also commute with
the Hamiltonian, $[H,u_{ij}^{\mu}]=0$ and thus it is also a constant of
motion. This representation is redundant, and the physical states need to be
eigenstates of $D_{i}=ib^{1}_{i}b^{2}_{i}b^{3}_{i}b^{4}_{i}c_{i}d_{i}$ with
eigenvalues $+1$. These constraints can be implemented by a projection
operator $P=\prod_{i}(1+D_{i})/2$. A $\mathbb{Z}_{2}$ gauge transformation at
site $i$ involves flipping the signs of the Majorana fermions and bond
operators,
$\\{c_{i},d_{i}\\}\rightarrow\\{-c_{i},-d_{i}\\};\leavevmode\nobreak\
u^{\mu}_{ij}\rightarrow-u^{\mu}_{ij}$. The plaquette operators can be
expressed in terms of the bond operators as $W_{p}=(-i)^{n}\prod_{(j,k)\in
p}u_{jk}^{\mu}$ where n is the number of the links on the boundary of the
plaquette p and the product is taken in the counter-clockwise direction. The
eigenvalues of $W_{p}$ for even (odd) plaquettes can take values $\pm 1$ ($\pm
i$). Therefore the solution of the model involves two flavors of free Majorana
fermions hopping in the background of static $\mathbb{Z}_{2}$ fluxes. Since
there is no Lieb’s theorem for quasicrystals, we obtain the ground state phase
diagram via Monte Carlo simulations for approximants on varying sizes, s=1, 2,
and 3. We use the Metropolis algorithm for the flux degrees of freedom, and
for a given vison configuration, we perform exact diagonalization to obtain
the corresponding energy. We perform 5000 Monte Carlo steps for a given
temperature and reduce the temperature down to $10^{-3}J$. Motivated by the
Monte Carlo results, we also construct various variational states and compare
their energies. For instance, in the case of finite $J_{5}$ calculations, we
construct 32 vison configurations for which 21 of them are stabilized as a
function of $J_{5}/J$.

##  IV Results and Discussion

We first discuss the phase diagram in the absence of the onsite anisotropy
term ($J_{5}=0$). Even though Lieb’s theorem[52] does not apply to
quasicrystals, we observe that the ground state follows a simple rule that is
an extension of Lieb’s theorem, and the fluxes through the plaquettes,
$W_{p}$, are given by

$\displaystyle\phi_{p}=-(\pm i)^{n}$ (6)

which recovers the well-known result that honeycomb (square, octagon)
plaquettes exhibit 0 ($\pi$) fluxes whereas the triangle, heptagon (pentagon)
plaquettes have $\pm\pi/2$ ($\mp\pi/2$) fluxes in the ground state as shown in
Fig. 3. We will refer this vison configuration ‘L’ for the remainder of the
article. Albeit there is no proof for the extension of Lieb’s theorem for
graphs that has plaquettes with varying number of edge sites, n, similar
behaviors have also been observed in amorphous and polycrystalline Kitaev
models[18, 19] and Kitaev models in closed-geometries (i. e.
tetrahedrons)[62]. Eq. 6 implies that the ground state flux configuration
spontaneously breaks time-reversal symmetry due to the odd-numbered plaquettes
(n=3, 5, 7). This phenomenon was originally pointed out by Kitaev[9] and was
first discovered in a model by Yao and Kivelson[55]. A time reversal operation
on a ground state flips the sign of the fluxes on the odd-numbered plaquettes
and generates a new ground state.

Figure 5: (a) Edge states without the onsite term for a point in the phase
space characterized by the Chern number 2. (b) Edges states for the same point
with onsite term splits the edge states but preserves $\nu=2$. Figure 6: Phase
diagram as a function of the onsite anisotropy, $J_{5}/J$ for
$J_{1}=J_{2}=J_{3}=J_{4}=J$. The y-axis represents the gap, $\Delta$. We
obtain 21 different vison configurations shown in Table I.

We find that the L vison configuration remains the ground state even with
anisotropic couplings, $J_{\mu}$, as long as $J_{5}=0$. However, similar to
Kitaev’s original model, the Majorana band structure depends on the values of
$J_{\mu}$’s. To investigate that, we study the phase diagram for
$J_{1}+J_{2}+J_{3}+J_{4}=1$, which forms a tetrahedron phase space. For the
sake of a simpler and clearer presentation, we take four cuts on the
tetrahedron as shown in Fig. 4. We find two gapped phases: a chiral spin
liquid with $\nu=2$ and an abelian spin liquid with $\nu=0$. The chiral edge
modes for the $\nu=2$ phase originate from the two flavors of Majorana
fermions both contributing a single chiral edge mode. For $J_{5}=0$, the two
Majorana fermions are completely decoupled and have an identical spectrum. We
identify the Chern number by counting the chiral edge modes on a slab
geometry. An example of the chiral edge modes are shown in Fig. 5. The gapped
phases are connected by a gap-closing transition and in certain cases a finite
region of the gapless phase (see Fig 4(a). We note that our phase diagram
shows similarities to the amorphous Kitaev model[18, 19] which has $\nu=0$ and
$\nu=1$ phases.

Flux sector | $n=3$ | $n=4$ | $n=5$ | $n=6$ | $n=7$ | $n=8$
---|---|---|---|---|---|---
$L$ | $\pi/2$ | $\pi$ | $-\pi/2$ | 0 | $+\pi/2$ | $\pi$
$R_{1}$ | $\pi/2$ | $\pi$ | $-\pi/2$ | 0 | $+\pi/2$ | 0
$R_{2}$ | $\pi/2$ | $\pi$ | $-\pi/2$ | 0 | $-\pi/2$ | 0
$R_{3}$ | $\pi/2$ | $\pi$ | $-\pi/2$ | $\pi$ | $-\pi/2$ | 0
$R_{4}$ | $\pi/2$ | $\pi$ | $+\pi/2$ | $\pi$ | $-\pi/2$ | 0
$R_{5}$ | $\pi/2$ | $\pi$ | $+\pi/2$ | $\pi$ | $-\pi/2$ | $\pi$
$R_{6}$ | $\pi/2$ | $\pi$ | $+\pi/2$ | $\pi$ | $+\pi/2$ | $\pi$
$R_{7}$ | $\pi/2$ | $0$ | $-\pi/2$ | 0 | $-\pi/2$ | 0
$R_{8}$ | $\pi/2$ | $0$ | $-\pi/2$ | 0 | $-\pi/2$ | $\pi$
$R_{9}$ | $\pi/2$ | $0$ | $-\pi/2$ | $\pi$ | $-\pi/2$ | $\pi$
$R_{10}$ | $\pi/2$ | $0$ | $+\pi/2$ | 0 | $-\pi/2$ | $\pi$
$R_{11}$ | $\pi/2$ | $0$ | $+\pi/2$ | $\pi$ | $-\pi/2$ | 0
$R_{12}$ | $\pi/2$ | $0$ | $+\pi/2$ | $\pi$ | $+\pi/2$ | $\pi$
$R_{13}$ | $\pi/2$ | $0$ | $+\pi/2$ | $\pi$ | $+\pi/2$ | 0
$R_{14}$ | $\pi/2$ | $0$ | $+\pi/2$ | 0 | $+\pi/2$ | 0
$R_{15}$ | $\pi/2$ | $0$ | $+\pi/2$ | 0 | $-\pi/2$ | 0
$R_{16}$ | $\pi/2$ | $0$ | $+\pi/2$ | 0 | $-\pi/2$ | $\pi$
$R_{17}$ | $\pi/2$ | $0$ | $-\pi/2$ | 0 | $-\pi/2$ | $\pi$
$R_{18}$ | $\pi/2$ | $0$ | $-\pi/2$ | 0 | $+\pi/2$ | $\pi$
$R_{19}$ | $\pi/2$ | $0$ | $-\pi/2$ | 0 | $+\pi/2$ | 0
$R_{20}$ | $\pi/2$ | $0$ | $-\pi/2$ | $\pi$ | $+\pi/2$ | 0
Table 1: Flux Configurations that are stabilized as a function of $J_{5}/J$ in
Fig. 6.

Next, we discuss the effects of the $J_{5}$ term which is an onsite term that
couples to $\Gamma^{5}$. Since $\Gamma^{\mu}$ operators can be expressed as
$J=3/2$ quadruple operators, $J_{5}$ term does not break time-reversal
symmetry but it is akin to an anisotropy term. Previous studies[57, 63, 64]
show that similar onsite terms can change the ground state vison
configuration. In terms of Majorana fermion description, $J_{5}$ term couples
the two flavors, $c$ and $d$ (see Eq. 5). We consider uniform bond coupling
constants, $J_{1}=J_{2}=J_{3}=J_{4}=J$ and study the $J_{5}/J$ phase diagram.
Our Monte Carlo simulations suggest that new vison configurations can be
stabilized. However, these configurations all obey a simple rule that
plaquettes with the same $n$ have the same flux configuration. Since dABL has
$n=\\{3,4,5,6,7,8\\}$ plaquettes, there are $2^{6}$ vison configuration
possibilities. As each configuration has a time-reversal partner and there are
a total of 32 distinct possible vison configurations that follow this rule.
Since the energy differences between the vison configurations can be quite
small, obtaining the true ground state from Monte Carlo can be challenging due
to phase separation. Therefore, we construct 32 variational configurations and
compare their energies. In Fig. 6, we show that 21 of these states can be
realized as a function of $J_{5}/J_{1}$. The vison configurations of these
phases are shown in Table I. Note that their time-reversal partners are also
ground states. For $J_{5}/J_{1}>3.5$ all vison configurations become
degenerate. This phase is unstable to confinement and has been observed in
similar models with onsite anisotropy[57]. For $J_{1}=J_{2}=J_{3}=J_{4}$ and
$J_{5}=0$, the ground state has $\nu=\pm 2$. We find that $J_{5}$ term
hybridizes the edge states but does not change $\nu$. Remarkably, all 21
states in Table I have the same Chern number $\nu=\pm 2$ even though they have
different vison configurations.

##  V Conclusion

We developed an exactly solvable spin liquid model on a quasicrystal lattice.
The structural properties required for Kitaev type integrability are uniform
coordination number $z$ throughout the lattice and the bond lattice being
$z$-partite. The dual lattice (center model) of the commonly used ABL is
uniformly coordinated with 4 bonds at each site. We showed that a coloring
rule applied to the four dimensional simple cubic lattice generates the dABL
with the required 4-partite property. With these structural properties we used
the $\Gamma$ matrix generalization of the Kitaev interactions to define our
Hamiltonian. A Majorana fermion representation of the interactions reduce the
Hamiltonian to a tight binding model on the quasicrystal that is coupled to a
static $Z_{2}$ gauge field. We used both finite size quasicrystals with open
boundary conditions and periodic approximants for our numerical calculations.

Figure 7: A Tesseract colored with four different colors so that no square
mesh has colors of the same side.

We obtained the ground state phase diagram via the Monte Carlo simulations and
variational analysis and showed that it follows a flux configuration that is
an extension of the Lieb’s theorem, albeit the Lieb’s theorem does not
directly apply here. As a function of coupling strengths, we find that the
ground state can be gapped with $\nu=\pm 2$ or $\nu=0$. These phases are
separated by a gapless phase. Subsequently, including an onsite term that
preserves the integribility of the model, we found that 21 different vison
configurations can be stabilized for isotropic exchange constants. Notably,
all of these phases also share the same Chern number, $\nu=2$. Interesting
future directions include extending our formalism to the emerging field of
moiré quasicrystals.

##  VI Acknowledgements

We thank Yuan-Ming Lu and Johannes Knolle for fruitful discussions. OE
acknowledges support from NSF Award No. DMR 2234352. MAK and MOO are supported
by the TUBITAK 1001 program Grant no 122F346.

##  Appendix A Coloring Rules

The cut-and-project construction defines the relation between the four-
dimensional simple cubic lattice and the ABL. Both the square and rhombus
meshes of the ABL correspond to squares in the two-planes of the four-
dimensional lattice. Therefore, we derive the coloring rules by considering a
four-dimensional hypercube shown in Fig. 7. We first color the intersecting
edges at any corner of the tesseract with four different colors. This coloring
scheme is then extended to the full simple cubic lattice by repeating
reflected copies of the tesseract.

A link parallel to $\hat{e}_{n}$ connects two sites that have the same three
indices $k_{m\neq n}$. We obtain the coloring rules in Table 2.

Table 2: Coloring rules for the ABL so that each mesh has four edges with
different colors.

$\hat{e}_{0}$ Bonds
---
Color | $k_{1}$ | $k_{2}$ | $k_{3}$
1 | even | odd | odd
| odd | even | even
2 | odd | odd | odd
| even | even | even
3 | even | odd | even
| odd | even | odd
4 | odd | odd | even
| even | even | odd

$\hat{e}_{1}$ Bonds
---
Color | $k_{0}$ | $k_{2}$ | $k_{3}$
1 | even | even | odd
| odd | odd | even
2 | odd | even | odd
| even | odd | even
3 | odd | even | even
| even | odd | odd
4 | even | even | even
| odd | odd | odd

$\hat{e}_{2}$ Bonds
---
Color | $k_{0}$ | $k_{1}$ | $k_{3}$
1 | odd | odd | odd
| even | even | even
2 | even | even | odd
| odd | odd | even
3 | odd | even | odd
| even | odd | even
4 | even | odd | odd
| odd | even | even

$\hat{e}_{3}$ Bonds
---
Color | $k_{0}$ | $k_{1}$ | $k_{2}$
1 | even | odd | odd
| odd | even | even
2 | odd | even | odd
| even | odd | even
3 | odd | odd | odd
| even | even | even
4 | odd | odd | even
| even | even | odd

In all our figures, 1 $\equiv$ black, 2 $\equiv$ yellow, 3 $\equiv$ green, 4
$\equiv$ red. The color of a link in the dual lattice is the color of the link
that intersects it in the ABL.

## References

  * Zhou  _et al._ [2017] Y. Zhou, K. Kanoda, and T.-K. Ng, Quantum spin liquid states, Rev. Mod. Phys. 89, 025003 (2017).
  * Clark and Abdeldaim [2021] L. Clark and A. H. Abdeldaim, Quantum spin liquids from a materials perspective, Annual Review of Materials Research 51, 495 (2021), https://doi.org/10.1146/annurev-matsci-080819-011453 .
  * Wen  _et al._ [2019] J. Wen, S.-L. Yu, S. Li, W. Yu, and J.-X. Li, Experimental identification of quantum spin liquids, npj Quantum Materials 4, 12 (2019).
  * Wen [2002] X.-G. Wen, Quantum orders and symmetric spin liquids, Phys. Rev. B 65, 165113 (2002).
  * Balents [2010] L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 (2010).
  * Wen [2017] X.-G. Wen, Colloquium: Zoo of quantum-topological phases of matter, Rev. Mod. Phys. 89, 041004 (2017).
  * Knolle and Moessner [2019] J. Knolle and R. Moessner, A field guide to spin liquids, Annual Review of Condensed Matter Physics 10, 451 (2019).
  * Broholm  _et al._ [2020] C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil, Quantum spin liquids, Science 367, eaay0668 (2020).
  * Kitaev [2006] A. Kitaev, Anyons in an exactly solved model and beyond, Annals of Physics 321, 2 (2006), january Special Issue.
  * Hwan Chun  _et al._ [2015] S. Hwan Chun, J.-W. Kim, J. Kim, H. Zheng, C. C. Stoumpos, C. D. Malliakas, J. F. Mitchell, K. Mehlawat, Y. Singh, Y. Choi, T. Gog, A. Al-Zein, M. M. Sala, M. Krisch, J. Chaloupka, G. Jackeli, G. Khaliullin, and B. J. Kim, Direct evidence for dominant bond-directional interactions in a honeycomb lattice iridate na2iro3, Nature Physics 11, 462 (2015).
  * Takagi  _et al._ [2019] H. Takagi, T. Takayama, G. Jackeli, G. Khaliullin, and S. E. Nagler, Concept and realization of kitaev quantum spin liquids, Nature Reviews Physics 1, 264 (2019).
  * Kitagawa  _et al._ [2018] K. Kitagawa, T. Takayama, Y. Matsumoto, A. Kato, R. Takano, Y. Kishimoto, S. Bette, R. Dinnebier, G. Jackeli, and H. Takagi, A spin–orbital-entangled quantum liquid on a honeycomb lattice, Nature 554, 341 (2018).
  * Lee  _et al._ [2020] I. Lee, F. G. Utermohlen, D. Weber, K. Hwang, C. Zhang, J. van Tol, J. E. Goldberger, N. Trivedi, and P. C. Hammel, Fundamental spin interactions underlying the magnetic anisotropy in the kitaev ferromagnet ${\mathrm{cri}}_{3}$, Phys. Rev. Lett. 124, 017201 (2020).
  * Kao and Perkins [2021] W.-H. Kao and N. B. Perkins, Disorder upon disorder: Localization effects in the kitaev spin liquid, Annals of Physics 435, 168506 (2021), special issue on Philip W. Anderson.
  * Dantas and Andrade [2022] V. Dantas and E. C. Andrade, Disorder, low-energy excitations, and topology in the kitaev spin liquid, Phys. Rev. Lett. 129, 037204 (2022).
  * Kao and Perkins [2022] W.-H. Kao and N. B. Perkins, Spin- and flux-gap renormalization in the random kitaev spin ladder, Phys. Rev. B 106, L100402 (2022).
  * Singhania  _et al._ [2023] A. Singhania, J. van den Brink, and S. Nishimoto, Disorder effects in the kitaev-heisenberg model, Phys. Rev. Res. 5, 023009 (2023).
  * Cassella  _et al._ [2022] G. Cassella, P. D’Ornellas, T. Hodson, W. M. H. Natori, and J. Knolle, An exact chiral amorphous spin liquid (2022), arXiv:2208.08246 [cond-mat.str-el] .
  * Grushin and Repellin [2023] A. G. Grushin and C. Repellin, Amorphous and polycrystalline routes toward a chiral spin liquid, Phys. Rev. Lett. 130, 186702 (2023).
  * Shechtman  _et al._ [1984] D. Shechtman, I. Blech, D. Gratias, and J. W. Cahn, Metallic phase with long-range orientational order and no translational symmetry, Phys. Rev. Lett. 53, 1951 (1984).
  * Levine and Steinhardt [1986] D. Levine and P. J. Steinhardt, Quasicrystals. i. definition and structure, Phys. Rev. B 34, 596 (1986).
  * Levine and Steinhardt [1984] D. Levine and P. J. Steinhardt, Quasicrystals: A new class of ordered structures, Phys. Rev. Lett. 53, 2477 (1984).
  * Arai  _et al._ [1988] M. Arai, T. Tokihiro, T. Fujiwara, and M. Kohmoto, Strictly localized states on a two-dimensional penrose lattice, Phys. Rev. B 38, 1621 (1988).
  * Mirzhalilov and Oktel [2020] M. Mirzhalilov and M. O. Oktel, Perpendicular space accounting of localized states in a quasicrystal, Phys. Rev. B 102, 064213 (2020).
  * Sutherland [1986] B. Sutherland, Self-similar ground-state wave function for electrons on a two-dimensional penrose lattice, Phys. Rev. B 34, 3904 (1986).
  * Kalugin and Katz [2014] P. Kalugin and A. Katz, Electrons in deterministic quasicrystalline potentials and hidden conserved quantities, Journal of Physics A: Mathematical and Theoretical 47, 315206 (2014).
  * Varjas  _et al._ [2019] D. Varjas, A. Lau, K. Pöyhönen, A. R. Akhmerov, D. I. Pikulin, and I. C. Fulga, Topological phases without crystalline counterparts, Phys. Rev. Lett. 123, 196401 (2019).
  * Fan and Huang [2021] J. Fan and H. Huang, Topological states in quasicrystals, Frontiers of Physics 17, 13203 (2021).
  * Chen  _et al._ [2020] R. Chen, C.-Z. Chen, J.-H. Gao, B. Zhou, and D.-H. Xu, Higher-order topological insulators in quasicrystals, Phys. Rev. Lett. 124, 036803 (2020).
  * Poon [1992] S. Poon, Electronic properties of quasicrystals an experimental review, Advances in Physics 41, 303 (1992), https://doi.org/10.1080/00018739200101513 .
  * Deguchi  _et al._ [2012] K. Deguchi, S. Matsukawa, N. K. Sato, T. Hattori, K. Ishida, H. Takakura, and T. Ishimasa, Quantum critical state in a magnetic quasicrystal, Nature Materials 11, 1013 (2012).
  * Kamiya  _et al._ [2018] K. Kamiya, T. Takeuchi, N. Kabeya, N. Wada, T. Ishimasa, A. Ochiai, K. Deguchi, K. Imura, and N. K. Sato, Discovery of superconductivity in quasicrystal, Nature Communications 9, 154 (2018).
  * Tamura  _et al._ [2021] R. Tamura, A. Ishikawa, S. Suzuki, T. Kotajima, Y. Tanaka, T. Seki, N. Shibata, T. Yamada, T. Fujii, C.-W. Wang, M. Avdeev, K. Nawa, D. Okuyama, and T. J. Sato, Experimental observation of long-range magnetic order in icosahedral quasicrystals, Journal of the American Chemical Society 143, 19938 (2021), pMID: 34786934, https://doi.org/10.1021/jacs.1c09954 .
  * Förster  _et al._ [2013] S. Förster, K. Meinel, R. Hammer, M. Trautmann, and W. Widdra, Quasicrystalline structure formation in a classical crystalline thin-film system, Nature 502, 215 (2013).
  * Collins  _et al._ [2017] L. C. Collins, T. G. Witte, R. Silverman, D. B. Green, and K. K. Gomes, Imaging quasiperiodic electronic states in a synthetic penrose tiling, Nature Communications 8, 15961 (2017).
  * Vardeny  _et al._ [2013] Z. V. Vardeny, A. Nahata, and A. Agrawal, Optics of photonic quasicrystals, Nature Photonics 7, 177 (2013).
  * Viebahn  _et al._ [2019] K. Viebahn, M. Sbroscia, E. Carter, J.-C. Yu, and U. Schneider, Matter-wave diffraction from a quasicrystalline optical lattice, Phys. Rev. Lett. 122, 110404 (2019).
  * Uri  _et al._ [2023] A. Uri, S. C. de la Barrera, M. T. Randeria, D. Rodan-Legrain, T. Devakul, P. J. D. Crowley, N. Paul, K. Watanabe, T. Taniguchi, R. Lifshitz, L. Fu, R. C. Ashoori, and P. Jarillo-Herrero, Superconductivity and strong interactions in a tunable moiré quasiperiodic crystal, arXiv:2302.00686 (2023).
  * Komura and Okabe [2016] Y. Komura and Y. Okabe, High-precision monte carlo simulation of the ising models on the penrose lattice and the dual penrose lattice, Journal of the Physical Society of Japan 85, 044004 (2016), https://doi.org/10.7566/JPSJ.85.044004 .
  * Thiem and Chalker [2015] S. Thiem and J. T. Chalker, Long-range magnetic order in models for rare-earth quasicrystals, Phys. Rev. B 92, 224409 (2015).
  * Flicker  _et al._ [2020] F. Flicker, S. H. Simon, and S. A. Parameswaran, Classical dimers on penrose tilings, Phys. Rev. X 10, 011005 (2020).
  * Lloyd  _et al._ [2022] J. Lloyd, S. Biswas, S. H. Simon, S. A. Parameswaran, and F. Flicker, Statistical mechanics of dimers on quasiperiodic ammann-beenker tilings, Phys. Rev. B 106, 094202 (2022).
  * Konstantinidis [2005] N. P. Konstantinidis, Antiferromagnetic heisenberg model on clusters with icosahedral symmetry, Phys. Rev. B 72, 064453 (2005).
  * Wessel  _et al._ [2003] S. Wessel, A. Jagannathan, and S. Haas, Quantum antiferromagnetism in quasicrystals, Phys. Rev. Lett. 90, 177205 (2003).
  * Jagannathan  _et al._ [2007] A. Jagannathan, A. Szallas, S. Wessel, and M. Duneau, Penrose quantum antiferromagnet, Phys. Rev. B 75, 212407 (2007).
  * Koga and Tsunetsugu [2017] A. Koga and H. Tsunetsugu, Antiferromagnetic order in the hubbard model on the penrose lattice, Phys. Rev. B 96, 214402 (2017).
  * Ghosh [2023] P. Ghosh, Exact quantum ground state of a two-dimensional quasicrystalline antiferromagnet (2023), arXiv:2301.11331 [cond-mat.str-el] .
  * Penrose [1974] R. Penrose, The role of aesthetics in pure and applied mathematical research, Bull.Inst.Math.Appl. 10, 266 (1974).
  * Ammann  _et al._ [1992] R. Ammann, B. Grünbaum, and G. C. Shephard, Aperiodic tiles, Discrete & Computational Geometry 8, 1 (1992).
  * Beenker [1982] F. Beenker,  _Algebraic theory of non-periodic tilings of the plane by two simple building blocks : a square and a rhombus_ , EUT report. WSK, Dept. of Mathematics and Computing Science (Eindhoven University of Technology, 1982).
  * Fujiwara  _et al._ [1988] T. Fujiwara, M. Arai, T. Tokihiro, and M. Kohmoto, Localized states and self-similar states of electrons on a two-dimensional penrose lattice, Phys. Rev. B 37, 2797 (1988).
  * Lieb [1994] E. H. Lieb, Flux phase of the half-filled band, Phys. Rev. Lett. 73, 2158 (1994).
  * Duneau [1989] M. Duneau, Approximants of quasiperiodic structures generated by the inflation mapping, Journal of Physics A: Mathematical and General 22, 4549 (1989).
  * Wu  _et al._ [2009] C. Wu, D. Arovas, and H.-H. Hung, $\Gamma$-matrix generalization of the kitaev model, Phys. Rev. B 79, 134427 (2009).
  * Yao and Kivelson [2007] H. Yao and S. A. Kivelson, Exact chiral spin liquid with non-abelian anyons, Phys. Rev. Lett. 99, 247203 (2007).
  * Nakai  _et al._ [2012] R. Nakai, S. Ryu, and A. Furusaki, Time-reversal symmetric kitaev model and topological superconductor in two dimensions, Phys. Rev. B 85, 155119 (2012).
  * Chulliparambil  _et al._ [2021] S. Chulliparambil, L. Janssen, M. Vojta, H.-H. Tu, and U. F. P. Seifert, Flux crystals, majorana metals, and flat bands in exactly solvable spin-orbital liquids, Phys. Rev. B 103, 075144 (2021).
  * Yao and Lee [2011] H. Yao and D.-H. Lee, Fermionic magnons, non-abelian spinons, and the spin quantum hall effect from an exactly solvable spin-$1/2$ kitaev model with su(2) symmetry, Phys. Rev. Lett. 107, 087205 (2011).
  * Vijayvargia  _et al._ [2023] A. Vijayvargia, E. M. Nica, R. Moessner, Y.-M. Lu, and O. Erten, Magnetic fragmentation and fractionalized goldstone modes in a bilayer quantum spin liquid, Phys. Rev. Res. 5, L022062 (2023).
  * Nica  _et al._ [2023] E. Nica, M. Akram, A. Vijayvargia, R. Moessner, and O. Erten, Kitaev spin-orbital bilayers and their moiré superlattices, npj Quantum Mater. 8, 9 (2023).
  * Kugel' and Khomskiĭ [1982] K. I. Kugel' and D. I. Khomskiĭ, The jahn-teller effect and magnetism: transition metal compounds, Soviet Physics Uspekhi 25, 231 (1982).
  * Mellado  _et al._ [2015] P. Mellado, O. Petrova, and O. Tchernyshyov, Projective symmetry of partons in the kitaev honeycomb model, Phys. Rev. B 91, 041103 (2015).
  * Chua  _et al._ [2011] V. Chua, H. Yao, and G. A. Fiete, Exact chiral spin liquid with stable spin fermi surface on the kagome lattice, Phys. Rev. B 83, 180412 (2011).
  * Akram  _et al._ [2023] M. Akram, E. M. Nica, Y.-M. Lu, and O. Erten, Vison crystals, chiral and crystalline phases in the yao-lee model, arXiv:2304.09891 (2023).


# Geometric Integration on Symmetric Spaces

Hans Munthe-Kaas  Lie Størmer Center, UiT The Arctic University of Norway
and Department of Mathematics, University of Bergen<EMAIL_ADDRESS>
(Date: September 3, 2024)

###### Abstract.

We consider geometric numerical integration algorithms for differential
equations evolving on symmetric spaces. The integrators are constructed from
canonical operations on the symmetric space, its Lie triple system (LTS), and
the exponential from the LTS to the symmetric space. Examples of symmetric
spaces are $n$-spheres and Grassmann manifolds, the space of positive definite
symmetric matrices, Lie groups with a symmetric product, and elliptic and
hyperbolic spaces with constant sectional curvatures. We illustrate the
abstract algorithm with concrete examples. In particular for the $n$-sphere
and the $n$-dimensional hyperbolic space the resulting algorithms are very
simple and cost only ${\mathcal{O}}(n)$ operations per step.

##  1\. Introduction

_Symmetric spaces_ are fundamental geometric objects in mathematics
originating from the seminal work on non-Euclidean geometries by Gauss,
Bolyai, Lobachevsky, Riemann, Klein and Lie in the 19’th century. Examples
include:

  * •

The $n$-sphere $S^{n}$, the surface of the unit ball in $\mathbb{R}^{n+1}$.

  * •

Hyperbolic spaces $H^{n}$ with constant negative sectional curvature.

  * •

Symmetric positive definite matrices with $A\cdot B=AB^{-1}A$.

  * •

Lie groups $G$ with a symmetric product $g\cdot h=gh^{-1}g$ and subsets of $G$
closed under $\cdot$.

  * •

Homogeneous spaces ${\mathcal{M}}=G/H$ where the Lie algebra of $G$ splits as
${\mathfrak{g}}={\mathfrak{h}}\oplus{\mathfrak{m}}$ such that ${\mathfrak{h}}$
is the Lie algebra of $H$ and ${\mathfrak{m}}$ is a Lie triple system.

  * •

The real line $\mathbb{R}$ with $x\cdot y=2x-y$.

  * •

Grassman manifolds: the $p$-dimensional subspaces of $\mathbb{R}^{n}$.

Symmetric spaces have also appeared in various numerical algorithms, such as
e.g. matrix computations and splines [16, 3, 5].

In this paper we are concerned with numerical algorithms for integrating
dynamical systems evolving on symmetric spaces. Over the last two decades
extensive theories have been developed for numerical integration of Lie groups
and homogeneous spaces, so-called _Lie group integrators_ [9]. A symmetric
space can always be constructed as a homogeneous space, so Lie group
integrators _can_ be used also for symmetric spaces. However, this approach
neglects the finer geometric structures of the symmetric space. In particular,
the natural metric of a symmetric space is preserved by a symmetric product,
but generally not by a (left- or right) Lie group action. Furthermore, the
tangent space of the symmetric space is a triple system closed under double
Lie brackets, but not under the single bracket of the Lie algebra. The
exponential map from the triple system to the symmetric space is often much
cheaper to compute than the corresponding exponential map from the Lie algebra
to the Lie group.

So, there are good reasons to consider _intrinsic_ integration algorithms on
symmetric spaces. In this paper we propose integration schemes entirely built
from operations which can be canonically obtained from the definition of a
symmetric space $({\mathcal{M}},\cdot)$ as a manifold with a symmetric
product. Algorithms are written out explicitly for several important examples.
For the $n$-sphere and hyperbolic $n$-space, the computational cost of the
algorithms are just ${\mathcal{O}}(n)$ per step, using Rodrigues type
formulae. The concrete algorithms for these cases can easily be implemented
without comprehending the theoretical framework of this paper.

##  2\. Canonical integration on symmetric spaces

###  2.1. An outline of the algorithm

We discuss the basic idea of canonical symmetric space integration in a
general geometric language, assuming some familiarity with symmetric spaces.
Detailed mathematical definitions are given in Section 2.2 and concrete
examples in Section 3.

The basic idea of the algorithm is very similar to the basic formulation of
the RKMK algorithm on a Lie group $G$, as presented in [15]. In Lie group case
we seek the solution of a differential equation $y^{\prime}(t)=f(y(t))y(t)$
for $y(t)\in G$, where $f\colon G\rightarrow{\mathfrak{g}}$ and
${\mathfrak{g}}=T_{e}G$ is the Lie algebra. The equation is pulled back from
$G$ to the linear space ${\mathfrak{g}}$ by the ansatz $y(t)=\exp(\theta(t))$
for $\theta(t)\in{\mathfrak{g}}$. A crucial point is that $\theta(t)$ is the
solution of the ’dexpinv’ equation
$\theta^{\prime}(t)=\operatorname{dexp}^{-1}_{\theta}f(\exp(\theta))$, where
$\operatorname{dexp}$ denotes a trivialised tangent of the exponential mapping
and $\operatorname{dexp}^{-1}_{\theta}$ can be computed by commutators on
${\mathfrak{g}}$. The RKMK Lie group integrator is obtained by solving the
dexpinv equation by a Runge–Kutta method, moving back and forth between
${\mathfrak{g}}$ and $G$ using $\exp$.

Unlike a Lie group, which has the identity $e$ as a special point, a symmetric
space has no special point. We can choose any point $o\in{\mathcal{M}}$ as a
_base point_ for various constructions, as long as we ensure that the
construction is geometrically independent of the chosen point. It should,
however, be remarked that from the perspective matrix representations and
efficient numerical linear algebra, certain base points might be preferred. We
call $({\mathcal{M}},\cdot,o)$ a _pointed symmetric space_. The tangent space
${\mathfrak{m}}:=T_{o}{\mathcal{M}}$ has the algebraic structure of a _Lie
Triple System_ , $({\mathfrak{m}},[\\_,\\_,\\_])$, where
$[\\_,\\_,\\_]\colon{\mathfrak{m}}\times{\mathfrak{m}}\times{\mathfrak{m}}\rightarrow{\mathfrak{m}}$
is the _triple bracket_ , defined below.

On ${\mathcal{M}}$ there exists a canonical connection $\nabla$ which is
torsion free $T=0$ and has constant (parallel) curvature $\nabla R=0$. The
algebra describing vector fields with the canonical connection on a symmetric
space, called a Lie Admissible Triple algebra, is studied in [18]. The
geodesic of this connection, starting at a tangent $V\in{\mathfrak{m}}$ is
denoted $\operatorname{Exp}\colon{\mathfrak{m}}\rightarrow{\mathcal{M}}$. The
connection defines _parallel transport_ of tangent vectors along the geodesic.
For a given $V\in{\mathfrak{m}}$, parallel transport yields a linear
isomorphism $\Gamma_{V}\colon{\mathfrak{m}}\rightarrow
T_{\operatorname{Exp}(V)}{\mathcal{M}}$, which transports a tangent vector
along the geodesic of $V$ from $t=0$ to $t=1$.

Crucial to our approach is the following characterisation of the tangent of
the geodesic exponential $\operatorname{Exp}$ trivialised by parallel
transport with $\Gamma$:

######  Proposition 2.1.

For $V,W\in{\mathfrak{m}}$ we have

(1)
$\left.\frac{d}{dt}\right|_{t=0}\operatorname{Exp}(V+tW)=\Gamma_{V}\operatorname{dExp}_{V}W,$

where $\operatorname{dExp}_{V}\colon{\mathfrak{m}}\rightarrow{\mathfrak{m}}$
is given as

(2)
$\operatorname{dExp}_{V}=\left.\frac{\sinh(\sqrt{x})}{\sqrt{x}}\right|_{x={\operatorname{ad}}^{2}_{V}}=\left.\sum_{n=0}^{\infty}\frac{x^{n}}{(2n+1)!}\right|_{x={\operatorname{ad}}^{2}_{V}},$

for ${\operatorname{ad}}^{2}_{V}W:=[W,V,V]$.

Proof: See [7], Th. 4.1 p.215. A proof is also provided here, in Section 3.1.

Consider the differential equation for $y(t)\in{\mathcal{M}}$ given a vector
field $F\colon{\mathcal{M}}\rightarrow T{\mathcal{M}}$:

(3) $y^{\prime}(t)=F(y(t)),\quad y(0)=o.$

Assuming $y(t)\in U$ is within a normal neighbourhood $o\in
U\subset{\mathcal{M}}$, where $\operatorname{Exp}$ is well-defined, we write

(4) $y(t)=\operatorname{Exp}(\theta(t)),\quad\theta(t)\in{\mathfrak{m}},$

and trivialise $F$ as

(5)
$F(\operatorname{Exp}(\theta))=\Gamma_{\theta}f(\operatorname{Exp}(\theta)),\quad\mbox{for
$f\colon U\rightarrow{\mathfrak{m}}$.}$

######  Proposition 2.2.

The curve $\theta(t)\in{\mathfrak{m}}$ satisfies the differential equation

(6)
$\theta^{\prime}(t)={\operatorname{dExp}}^{-1}_{\theta}f(\operatorname{Exp}(\theta)),\quad\theta(0)=0,$

for

(7)
${\operatorname{dExp}}^{-1}_{\theta}=\frac{\sqrt{x}}{\sinh(\sqrt{x})}=1-\sum_{n=1}^{\infty}(2^{2n}-2)\frac{B_{2n}}{(2n)!}x^{n}=1-\frac{x}{6}+\frac{7x^{2}}{360}-\frac{31x^{3}}{15120}+{\mathcal{O}}(x^{4}),$

where $x={\operatorname{ad}}^{2}_{\theta}=[\\_,\theta,\theta]$ and $B_{2n}$
are the Bernoulli numbers.

We obtain the first version of our algorithm by applying a Runge–Kutta (RK)
method to $\theta(t)$:

Result: Evolve $y^{\prime}(t)=F(y(t))$ from $y_{0}=o$ to $y_{1}\approx
y(t=h)$.

Choose time step $h$ and $\\{a_{i,j}\\}_{i,j=1}^{r}$, $\\{b_{j}\\}_{j=1}^{r}$
coefficients of an RK method.

for  _$i=1,\ldots,r$_ do

$\theta_{i}=\sum_{j=1}^{r}a_{i,j}K_{j}$

$K_{i}=h\
{\operatorname{dExp}}^{-1}_{\theta_{i}}\Gamma_{\theta_{i}}^{-1}F(\operatorname{Exp}(\theta_{i}))$

end for

$y_{1}=\operatorname{Exp}(\sum_{j=1}^{r}b_{j}K_{j})$

Algorithm 1 Canonical integration step on symmetric space
$({\mathcal{M}},\cdot,o)$, with $y_{0}=o$.

_Remark._ The perhaps most important case is _explicit RK methods_ , where the
coefficient matrix $\\{a_{i,j}\\}$ is strictly lower triangular. In this case
$\theta_{1}=0$ and we find the others as
$\theta_{i}=\sum_{j=1}^{i-1}a_{i,j}K_{j}$. If the method is _implicit_ , a
system of equations must be solved to find $K_{i}$ and $\theta_{i}$.

_Remark._ The ${\operatorname{dExp}}^{-1}$ series for the symmetric space
exponential starts with the double bracket, while for the Lie group
exponential the series starts with a single bracket. In [14] we showed that a
naïve Lie group method, with no correction for ${\operatorname{dExp}}^{-1}$
can generally only achieve order 2. The symmetric space algorithm implemented
without a ${\operatorname{dExp}}^{-1}$ correction can achieve order 3. With
higher order truncations of the ${\operatorname{dExp}}^{-1}$ expansion, this
method achieves the same order as the underlying RK method.

This basic version of the method chooses the base point $o=y_{0}$. If we want
to choose a different base point $o\neq y_{0}$, we can pull back the equation
using an automorphism of ${\mathcal{M}}$ sending $o$ to $y_{0}$. The
automorphisms $\operatorname{Aut}({\mathcal{M}})$ are all mappings
$\tau\colon{\mathcal{M}}\rightarrow{\mathcal{M}}$ such that $\tau(x\cdot
y)=\tau(x)\cdot\tau(y)$. Let $T\tau\colon T{\mathcal{M}}\rightarrow
T{\mathcal{M}}$ denote the tangent map.

Result: Evolve $y^{\prime}(t)=f(y(t))$ from $y_{0}$ to $y_{1}\approx y(t=h)$.

Choose any $\tau_{y_{0}}\in\operatorname{Aut}({\mathcal{M}})$ such that
$\tau_{y_{0}}(o)=y_{0}$.

Compute $\tilde{y}_{1}$ by applying Algorithm 1 on the pull-back equation
$\tilde{y}^{\prime}(t)=T\tau_{y_{o}}^{-1}F(\tau(\tilde{y}(t))),\quad\tilde{y}(0)=o.$

We return to a more detailed description of these algorithms in Section 2.2,
after a basic introduction to symmetric spaces.

###  2.2. Defining the canonical operations

We start with a pointed symmetric space $({\mathcal{M}},\cdot,o)$ and define
all the operations of the algorithm from this structure. Concrete examples are
found in Section 3. For a comprehensive treatment of the background theory, we
refer to [12, 7].

######  Definition 2.1 (Symmetric space).

A symmetric space is a manifold $({\mathcal{M}},\cdot)$ with a product
${\cdot\colon{\mathcal{M}}\times{\mathcal{M}}\rightarrow{\mathcal{M}}}$ such
that

$\displaystyle x\cdot x$ $\displaystyle=x$ $\displaystyle x\cdot(x\cdot y)$
$\displaystyle=y$ $\displaystyle x\cdot(y\cdot z)$ $\displaystyle=(x\cdot
y)\cdot(x\cdot z)$

and every $x$ has a neighbourhood $U$ such that $x\cdot y=y$ implies $y=x$ for
all $y\in U$.

As an example, any Lie group $L$ with the symmetric product $x\cdot
y:=xy^{-1}x$ is a symmetric space denoted $(L^{+},\cdot)$. An other example is
the sphere $S^{n}=\\{x\in\mathbb{R}^{n+1}\colon x^{T}x=1\\}$ with the product
$x\cdot y=y-2xx^{T}y$.

For each point $x\in{\mathcal{M}}$ we have the symmetry
$\sigma_{x}\colon{\mathcal{M}}\rightarrow{\mathcal{M}}$ defined as
$\sigma_{x}(y):=x\cdot y$. The symmetry
$\sigma_{x}\in\operatorname{Aut}({\mathcal{M}})$ is an involutive automorphism
of the symmetric space $({\mathcal{M}},\cdot)$, having $x$ as an isolated fix
point111To get a feeling for this structure, consider a Riemannian metric
space with its geodesics. Moving from $y$ to $x$ along the joining geodesic
twice the distance from $y$ to $x$, we arrive at $\sigma_{x}(y)$. So
$\sigma_{x}$ is a reflection of the geodesics through the point $x$. A
Riemannian symmetric space has for each $x\in{\mathcal{M}}$ such a symmetry
$\sigma_{x}$ which is a metric isometry and has $x$ as an isolated fix-point..
Note: ’involutive’: $\sigma_{x}(\sigma_{x}(y))=y$, ’automorphism’:
$\sigma_{x}(y\cdot z)=\sigma_{x}(y)\cdot\sigma_{x}(z)$, isolated ’fix-point
$x$’: $\sigma_{x}(x)=x$, are exactly the contents of Definition 2.1

The product of two reflections is called a _displacement_. The _group of
displacements_ is generated by all $\sigma_{x}\sigma_{y}$ for
$x,y\in{\mathcal{M}}$ and denoted $G=G({\mathcal{M}})$. This is a normal
subgroup of $\operatorname{Aut}({\mathcal{M}})$.

Let $o\in{\mathcal{M}}$ be a chosen point called the _base point_ , and let
${\mathfrak{m}}:=T_{o}{\mathcal{M}}$. We call $({\mathcal{M}},\cdot,o)$ a
_pointed symmetric space_. The map $Q\colon{\mathcal{M}}\rightarrow
G({\mathcal{M}})$ defined as $x\mapsto Q_{x}:=\sigma_{x}\sigma_{o}$, called
the _quadratic representation_. $Q_{\mathcal{M}}$ generates all of $G$ and it
is a homomorphism of $({\mathcal{M}},\cdot,o)$ onto $(G^{+},\cdot,e)$ with the
symmetric product $g\cdot h:=gh^{-1}g$, satisfying

$Q_{x\cdot y}=Q_{x}Q_{y}^{-1}Q_{x}=Q_{x}\cdot Q_{y},\quad Q_{o}=e\in G.$

The quadratic representation defines an action of ${\mathcal{M}}$ on itself,
$(x,y)\mapsto
Q_{x}y\colon{\mathcal{M}}\times{\mathcal{M}}\rightarrow{\mathcal{M}}$. By
differentiating $Q$ at $o$ we find the infinitesimal generator of the action
as the vector field $\xi_{V}\in{\mathcal{X}}{\mathcal{M}}$, defined for
$V\in{\mathfrak{m}}$ at $y\in{\mathcal{M}}$ as

(8)
$\xi_{V}(y)=\left.\frac{d}{dt}\right|_{t=0}Q_{\gamma(t)}(y),\quad\gamma(0)=o,\quad\gamma^{\prime}(0)=\frac{1}{2}V.$

The factor $\frac{1}{2}$ is due to the quadratic nature of $Q$. With this
scaling we have $\xi_{V}(o)=V$, so it defines an extension of a tangent at $o$
to a vector field on all of ${\mathcal{M}}$. Also recall that $Q_{y}o$ sends
$o$ beyond $y$ to the point on the opposite side of $y$ on the geodesic, so we
have to travel with half the speed to arrive at $y$ for $t=1$.

######  Definition 2.2 (Geodesic $\operatorname{Exp}$).

The geodesic exponential at the base point is a mapping
$\operatorname{Exp}\colon{\mathfrak{m}}\rightarrow{\mathcal{M}}$ defined as
the $t=1$ solution of the differential equation

(9) $y^{\prime}(t)=\xi_{V}(y(t)),\quad y(0)=o.$

Let $(\mathbb{R},\cdot)$ be the real line with symmetric product $s\cdot
t=2s-t$. It can be shown that ${\operatorname{Exp}(s\cdot
tV)=\operatorname{Exp}(sV)\cdot\operatorname{Exp}(tV)}$. We summarise the most
important properties of the symmetric space exponential:

######  Proposition 2.3.

$\operatorname{Exp}\colon(\mathbb{R},\cdot,0)\rightarrow({\mathcal{M}},\cdot,o)$
is the unique pointed symmetric space homomorphism with derivative
$V\in{\mathfrak{m}}$ at $0$. It relates to the quadratic representation via

$\operatorname{Exp}(V)=Q_{\operatorname{Exp}(\frac{1}{2}V)}o.$

Note that $Q_{\operatorname{Exp}(\frac{1}{2}tV)}o$ is a geodesic through $o$,
while $Q_{\operatorname{Exp}(\frac{1}{2}tV)}p$ is generally _not_ a geodesic
through $p$ for $p\neq o$. The maps
$Q_{\operatorname{Exp}(\frac{1}{2}tV)}\colon{\mathcal{M}}\rightarrow{\mathcal{M}}$
are $t$-parametrised families of automorphisms (isometries for Riemannian
symmetric spaces) of ${\mathcal{M}}$. In the case of the sphere $S^{2}$, the
family $Q_{\operatorname{Exp}(\frac{1}{2}tV)}\colon S^{2}\rightarrow S^{2}$
rotates the sphere around an axis $V$ perpendicular to the north pole.
Starting at the north pole, this yields a great circle (geodesic). Starting at
other points it yields (possibly) smaller circles orthogonal to the rotation
axis $V$.

For a fixed $y\in{\mathcal{M}}$, the tangent of the automorphism $p\mapsto
Q_{y}p$ is denoted $TQ_{y}\colon T_{p}{\mathcal{M}}\rightarrow
T_{Q_{y}p}{\mathcal{M}}$

(10)
$TQ_{y}(W)=\left.\frac{d}{dt}\right|_{t=0}TQ_{y}(\gamma(t)),\quad\gamma(0)=p,\gamma^{\prime}(0)=W.$

Restricted to $W\in{\mathfrak{m}}$, this defines parallel transport along the
geodesic:

######  Definition 2.3 (Parallel transport of $W$ along the geodesic of $V$).

For $V,W\in{\mathfrak{m}}$ define $\Gamma_{V}\colon{\mathfrak{m}}\rightarrow
T_{\operatorname{Exp}(V)}{\mathcal{M}}$ as

$\Gamma_{V}W=TQ_{\operatorname{Exp}(\frac{1}{2}V)}(W).$

In the introduction we mentioned that symmetric spaces are equipped with a
canonical connection
$\nabla\colon{\mathcal{X}}{\mathcal{M}}\times{\mathcal{X}}{\mathcal{M}}\rightarrow{\mathcal{X}}{\mathcal{M}}$
for which the name _’geodesic curves’_ acquires geometric meaning. For the
Levi–Civita connection on a Riemannian symmetric space, the geodesics define
(locally) shortest paths between points. The canonical connection is torsion
free and has constant curvature. We do not need $\nabla$ in this paper, so we
will be brief on its construction. For two vector fields
$F,G\in{\mathcal{X}}{\mathcal{M}}$ define the connection at the base point $o$
as the time derivative of the parallel transport,

$\nabla_{F}G(o)=\left.\frac{d}{dt}\right|_{t=0}\Gamma_{tF(o)}^{-1}G(\operatorname{Exp}(tF(o))).$

The definition can be extended to an arbitrary point $y\in{\mathcal{M}}$ by
moving around with $\tau\in\operatorname{Aut}({\mathcal{M}})$ such that
$\tau(o)=y$. We refer to [12, 10] for details.

We define the last operation we need to integrate on ${\mathcal{M}}$, the _Lie
triple system_ on ${\mathfrak{m}}$.

######  Definition 2.4 (Triple bracket on ${\mathfrak{m}}$).

The tri-linear bracket
$[\\_,\\_,\\_]\colon{\mathfrak{m}}\times{\mathfrak{m}}\times{\mathfrak{m}}\rightarrow{\mathfrak{m}}$
is defined as

(11) $[V,W,Z]=[[\xi_{V},\xi_{W}]_{J},\xi_{z}]_{J}(o),$

where $[\\_,\\_]_{J}$ is the Jacobi bracket of vector fields.

The algebra $({\mathfrak{m}},[\\_,\\_,\\_])$ has the following structure:

######  Definition 2.5 (Lie triple system).

A _Lie triple system_ (Lts) $({\mathcal{A}},[\cdot,\cdot,\cdot])$ is a vector
space ${\mathcal{A}}$ with a tri-linear bracket
$[\cdot,\cdot,\cdot]\colon{\mathcal{A}}\times{\mathcal{A}}\times{\mathcal{A}}\rightarrow{\mathcal{A}}$
satisfying for all $x,y,z,t,w\in{\mathcal{A}}$:

(12) $\displaystyle[x,x,z]=0$ (13) $\displaystyle[x,y,z]+[y,z,x]+[z,x,y]=0$
(14) $\displaystyle[x,y,[z,t,w]]=[[x,y,z],t,w]+[z,[x,y,t],w]+[z,t,[x,y,w]].$

We can combine Algorithms 1 and 2 in a single algorithm by using the quadratic
representation to pull back $y_{0}$ to ${\mathfrak{m}}$ and
$\operatorname{Exp}$ to pull back from a neigbourhood of $o\in{\mathcal{M}}$
to ${\mathfrak{m}}$. This yields:

######  Proposition 2.4.

Given $y^{\prime}(t)=F(y(t))$, $y(0)=y_{0}$. Let $s\in{\mathcal{M}}$ be such
that $Q_{s}o=y_{0}$ and let $\theta(t)\in{\mathfrak{m}}$ such that
$y(t)=Q_{s}(\operatorname{Exp}(\theta(t)))$. Then $\theta(t)$ satisfies

$\theta^{\prime}(t)={\operatorname{dExp}}^{-1}_{\theta}\Gamma_{\theta}^{-1}TQ_{s}^{-1}F(y(t)),\theta(0)=0.$

Using a RK method on this equation yields the following algorithm:

Result: Evolve $y^{\prime}(t)=F(y(t))$ from $y(0)=y_{0}$ to $y_{n}\approx
y(t\\!=\\!T)$.

Choose time step $h=T/n$ and $\\{a_{i,j}\\}_{i,j=1}^{r}$,
$\\{b_{j}\\}_{j=1}^{r}$ coefficients of an RK method.

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

solve $Q_{s_{\ell}}o=y_{\ell}$ for $s_{\ell}$

for  _$i=1,\ldots,r$_ do

$\theta_{i}=\sum_{j=1}^{r}a_{i,j}K_{j}$

$K_{i}=h\
{\operatorname{dExp}}^{-1}_{\theta_{i}}\Gamma_{\theta_{i}}^{-1}TQ_{s_{\ell}}^{-1}F\big{(}Q_{s_{\ell}}(\operatorname{Exp}(\theta_{i}))\big{)}$

end for

$y_{\ell+1}=Q_{s_{\ell}}\big{(}\operatorname{Exp}(\sum_{j=1}^{r}b_{j}K_{j})\big{)}$

end for

Algorithm 3 CSSI (Canonical Symmetric Space Integrator)
$({\mathcal{M}},\cdot,o)$.

##  3\. Examples

###  3.1. Lie groups

Let $L$ be a Lie group with Lie algebra
$({\mathfrak{g}},[\\_,\\_]_{{\mathfrak{g}}})$. For $V\in{\mathfrak{g}}$ and
$g\in G$ we write products in matrix style, e.g. $Vg$ instead of the more
elaborate $TR_{g}(V)$.

Let $\exp\colon{\mathfrak{g}}\rightarrow L$ denote the Lie group exponential.
Any Lie group $L$ defines a pointed symmetric space denoted
$L^{+}:=(L,\cdot,o)$, where $o=e$ is the identity of the Lie group and $x\cdot
y:=xy^{-1}x$. We obtain the quadratic representation
$Q_{x}=\sigma_{x}\sigma_{e}$ and

$Q_{x}y=xyx.$

For any $V\in{\mathfrak{g}}$ and $g\in L$, we have
$\xi_{V}(g)=\left.\frac{d}{dt}\right|_{t=0}Q_{\exp(\frac{1}{2}tV)}g=\frac{1}{2}\big{(}Vg+gV\big{)}$.
The solution of $y^{\prime}(t)=\xi_{V}(y(t))$, $y(0)=g$ is
$y(t)=\exp(\frac{1}{2}tV)g\exp(\frac{1}{2}tV)$, in particular for $g=e$ we
have $y(t)=\exp(tV)$, so we find

$\operatorname{Exp}(V)=\exp(V).$

Let $V^{r}(g)=Vg$ and $V^{\ell}(g)=gV$. We have $[V^{r},W^{\ell}]_{J}=0$,
$[V^{r},W^{r}]_{J}=-[V,W]_{\mathfrak{g}}^{r}$ and
$[V^{\ell},W^{\ell}]_{J}=[V,W]_{\mathfrak{g}}^{\ell}$. Since
$\xi_{V}=\frac{1}{2}(V^{r}+V^{\ell})$ we find the triple bracket

$[V,W,Z]=[[\xi_{V},\xi_{W}]_{J},\xi_{Z}]_{J}(o)=\frac{1}{4}[[V,W]_{\mathfrak{g}},Z]_{\mathfrak{g}}.$

The parallel transport is

$\Gamma_{V}W=\left.\frac{d}{dt}\right|_{t=0}Q_{\operatorname{Exp}(\frac{1}{2}V)}\operatorname{Exp}(tW)=\exp(\frac{1}{2}V)W\exp(\frac{1}{2}V).$

We compute $\operatorname{dExp}$ and $\operatorname{dExp}^{-1}$ explicitly.
From (1) we have for $V,W\in{\mathfrak{g}}$

$\operatorname{dExp}_{V}W=\Gamma_{V}^{-1}\left.\frac{d}{dt}\right|_{t=0}\operatorname{Exp}(V+tW)=\exp(-\frac{1}{2}V)\left.\frac{d}{dt}\right|_{t=0}\exp(V+tW)\exp(-\frac{1}{2}V).$

In the following computation, let the linear operators
$x,y:{\mathfrak{g}}\rightarrow{\mathfrak{g}}$ be given as
$xW=[V,W]_{\mathfrak{g}}$ and
$yW=[W,V,V]=\frac{1}{4}[V,[V,W]_{\mathfrak{g}}]_{\mathfrak{g}}=\frac{1}{4}x^{2}W$.
It is well known [9] that
$\left.\frac{d}{dt}\right|_{t=0}\exp(V+tW)=\operatorname{dexp}_{V}W\exp(V)$
where $\operatorname{dexp}_{V}=(e^{x}-1)/x$. This yields, using
$\operatorname{Ad}_{\exp(V)}=\exp({\operatorname{ad}}_{V})$:

$\displaystyle\operatorname{dExp}_{V}W$
$\displaystyle=\exp(-\frac{1}{2}V)\frac{e^{x}-1}{x}W\exp(V)\exp(-\frac{1}{2}V)=\operatorname{Ad}_{\exp(-\frac{1}{2}V)}\frac{e^{x}-1}{x}W=e^{-\frac{x}{2}}\frac{e^{x}-1}{x}W$
$\displaystyle=\frac{\sinh(\frac{x}{2})}{x/2}W=\frac{\sinh(\sqrt{y})}{\sqrt{y}}W.$

This establishes Proposition 2.1. It follows immediately that
${\operatorname{dExp}}^{-1}_{V}=\frac{\sqrt{y}}{\sinh(\sqrt{y})}$.

This results in the following matrix version of Algorithm 3. As a specific
application consider $y^{\prime}(t)=F(y(t))$, where $y(t)$ is a symmetric
positive definite (SPD) matrix and $F(y)$ is symmetric.

Result: Evolve $y^{\prime}(t)=F(y(t))$ from $y(0)=y_{0}$ to $y_{n}\approx
y(t\\!=\\!T)$.

Choose time step $h=T/n$ and $\\{a_{i,j}\\}_{i,j=1}^{r}$,
$\\{b_{j}\\}_{j=1}^{r}$ coefficients of an RK method.

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

solve $Q_{s_{\ell}}o=s_{\ell}s_{\ell}=y_{\ell}$ for $s_{\ell}$ (matrix square
root)

for  _$i=1,\ldots,r$_ do

$\theta_{i}=\sum_{j=1}^{r}a_{i,j}\tilde{K}_{j}$

$K_{i}=h\
\Gamma_{\theta_{i}}^{-1}TQ_{s_{\ell}}^{-1}F\big{(}Q_{s_{\ell}}(\operatorname{Exp}(\theta_{i}))\big{)}=h\
\exp(-\frac{\theta_{i}}{2})s_{\ell}^{-1}F\big{(}s_{\ell}\exp(\theta_{i})s_{\ell}\big{)}s_{\ell}^{-1}\exp(-\frac{\theta_{i}}{2})$

$\tilde{K}_{i}={\operatorname{dExp}}^{-1}_{\theta_{i}}K_{i}=K_{i}-\frac{x}{6}K_{i}+\frac{7x^{2}}{360}K_{i}-\frac{31x^{3}}{15120}K_{i}+\cdots$

end for

$\theta_{\ell}=\sum_{j=1}^{r}b_{j}\tilde{K}_{j}$

$y_{\ell+1}=Q_{s_{\ell}}\big{(}\operatorname{Exp}(\theta_{\ell})\big{)}=s_{\ell}\exp(\theta_{\ell})s_{\ell}$

end for

$({\mathcal{M}},\cdot,o)$, $o=I$, $A\cdot B=AB^{-1}A$.

Algorithm 4 CSGI (Canonical Symmetric Group Integrator)

$({\mathcal{M}},\cdot,o)$, $o=I$, $A\cdot B=AB^{-1}A$.

Here $xK_{i}=\frac{1}{4}[[K_{i},\theta_{i}],\theta_{i}]$ (matrix triple
commutator), and the ${\operatorname{dExp}}^{-1}$ series is truncated to the
order of the RK method. The matrix square root is in general not uniquely
defined. We return to the choice of square root in the end of Section 3.1.1.

####  3.1.1. Automorphisms on Lie groups and symmetric decompositions

For general matrix Lie groups, it is not clear to us if the symmetric space
integration algorithm above has any advantages over the conventional Lie group
integrators. However, there are many important cases where $L^{+}$ contains a
symmetric sub algebra. In these cases the conventional Lie group integrator
destroys this subspace structure, while it is respected by the symmetric
integration algorithm.

Given a Lie group $L$ with Lie algebra ${\mathfrak{g}}$ and the corresponding
symmetric space $L^{+}$. An involutive automorphism on $L$ is $S\colon
L\rightarrow L$ such that $S(xy)=S(x)S(y)$ and $S^{2}=\operatorname{Id}$.
There are two important subsets of $L$

$\displaystyle L^{S}$ $\displaystyle=\\{x\in L\colon S(x)=x\\}$ $\displaystyle
L_{S}$ $\displaystyle=\\{x\in L\colon S(x)=x^{-1}\\}.$

The invariant elements $L^{S}$ is a subgroup of $L$ and the alternating
elements $L_{S}$ is a symmetric subspace of $L^{+}$. Furthermore, $L_{S}$ is
isomorphic to the homogeneous space $L/L^{S}$ and _every_ connected symmetric
space ${\mathcal{M}}$ is isomorphic to a symmetric space of this kind. So this
is really more than just an example, it is a generic case covering all
connected symmetric spaces.

On the algebra level the derivative of $S$ at the identity yields an
involutive Lie algebra automorphism
$dS\colon{\mathfrak{g}}\rightarrow{\mathfrak{g}}$, and the algebra splits as
${\mathfrak{g}}={\mathfrak{h}}\oplus{\mathfrak{m}}$ where ${\mathfrak{h}}$ is
the $+1$ eigenspace of $dS$ and $m$ the $-1$ eigenspace. Here ${\mathfrak{h}}$
is a Lie algebra and ${\mathfrak{m}}$ a Lie triple system. The split parts
induces a $Z_{2}$ grading on ${\mathfrak{g}}$, where the sub spaces satisfy
$[{\mathfrak{h}},{\mathfrak{h}}]\subset{\mathfrak{h}}$,
$[{\mathfrak{h}},{\mathfrak{m}}]\subset{\mathfrak{m}}$ and
$[{\mathfrak{m}},{\mathfrak{m}}]\subset{\mathfrak{h}}$.

There are many structured problems in matrix computations which fit into this
format. A famous example is the group of all invertible matrices
$L=\operatorname{Gl}(n)$ with the automorphism $S(A)=A^{-T}$. Here $L^{S}$ is
the group of orthogonal matrices and $L_{S}$ the symmetric space of symmetric
non-singular matrices. An important sub algebra is the symmetric space of
symmetric positive definite (SPD) matrices. On the algebra level,
$dS(V)=-V^{T}$ splitting the space of matrices as
${\mathfrak{g}l}(n)={\mathfrak{h}}\oplus{\mathfrak{m}}$ where ${\mathfrak{h}}$
are the skew symmetric and ${\mathfrak{m}}$ the symmetric matrices. Other
examples are discussed in [17, 16].

We return to the matrix square root in Algorithm 4. The final step is

$y_{\ell+1}=s_{\ell}\exp(\frac{\theta_{\ell}}{2})\exp(\frac{\theta_{\ell}}{2})s_{\ell}=S(A_{\ell+1})^{-1}A_{\ell+1},$

where $A_{\ell+1}=\exp(\frac{\theta_{\ell}}{2})s_{\ell}$. The _generalised
polar decomposition_ [17] is $A_{\ell+1}=QP$, where $S(Q)=Q$ and
$S(P)=P^{-1}$. For the SPD example, this is the classical polar decomposition
where $Q$ is orthogonal and $P$ SPD. From this we find

$y_{\ell+1}=S(A_{\ell+1})^{-1}A_{\ell+1}=PP,$

and hence $s_{\ell+1}=P$, the polar part of $A_{\ell+1}$. This polar part
matrix square root is unique in many cases, such as for SPD matrices.

We finally note that in a dynamic situation, where $s_{\ell+1}$ is not too far
from $s_{\ell}$, iterative techniques for the polar decomposition could be
considered [8]. There are many such practical issues we will not pursue here,
but we will address these in forthcoming papers.

###  3.2. Spheres

Let $(S^{n},\cdot,o)$ be the pointed symmetric space
$S^{n}=\\{x\in\mathbb{R}^{n+1}\colon x^{T}x=1\\}$ with

$x\cdot y:=2xx^{T}y-y$

The base point could be the north pole $o=(0,0,\ldots,0,1)$ (or any other
point). The automorphisms are the orthogonal group
$\operatorname{Aut}(S^{n})=O(n+1)$ generated by the symmetries

(15) $\sigma_{x}=2{xx^{T}}-I\in O(n+1).$

The group of displacements $G(S^{n})=\operatorname{SO}(n+1)$ is the group of
pure rotations, and it is generated by the quadratic representation

$Q_{x}=\sigma_{x}\sigma_{o}=(2xx^{T}-I)(2oo^{T}-I),$

where $Q_{x}$ is the rotation in the $o-x$ plane through twice the angle from
$o$ to $x$. Let ${\mathfrak{m}}$ be the horizontal vectors at the north pole
(the equatorial plane)

${\mathfrak{m}}=T_{o}S^{n}=\\{v\in\mathbb{R}^{n+1}\colon
v^{T}o=0\\}\simeq\mathbb{R}^{n},$

where we identify with $\mathbb{R}^{n}$ by deleting last $0$ in the vector.
The infinitesimal generator is

$\xi_{V}(p)=\left.\frac{d}{dt}\right|_{t=0}Q_{o+tv}p=\left(vo^{T}-ov^{T}\right)p=\hat{v}p,$

where the hat map
$\hat{v}:=\left(vo^{T}-ov^{T}\right)\colon{\mathfrak{m}}\rightarrow{\mathfrak{s}o}(n+1)$
is the unique map such that $\hat{v}o=v$. From this we compute the triple
bracket and the double adjoint on ${\mathfrak{m}}$:

$\displaystyle[u,v,w]$
$\displaystyle=[[\hat{u},\hat{v}],\hat{w}]o=\big{(}vw^{T}-v^{T}wI\big{)}u$
$\displaystyle{\operatorname{ad}}^{2}_{v}w$
$\displaystyle=[w,v,v]=\big{(}vv^{T}-v^{T}vI\big{)}w.$

For $x={\operatorname{ad}}^{2}v$ we find $x^{2}=-v^{T}vx=(i\varphi)^{2}x$ for
$\varphi=||v||$. Since
$\frac{\sqrt{x}}{\sinh(x)}=1+\sum_{j=1}^{\infty}a_{j}x^{j}$, we find

${\operatorname{dExp}}^{-1}_{v}=\frac{\sqrt{x}}{\sinh(x)}=I+\left(\frac{i\varphi}{\sinh(i\varphi)}-1\right)\frac{x}{-\varphi^{2}}=I+\left(\frac{\varphi}{\sin(\varphi)}-1\right)\pi_{v}^{\perp},$

where $\pi_{v}^{\perp}=\frac{x}{-\varphi^{2}}=I-\frac{vv^{T}}{v^{T}v}$ is
projection on the orthogonal complement of $v$. This formula is not hard to
understand geometrically. In ${\operatorname{dExp}}^{-1}_{v}w$ the normal
component of $w$ is increased by the factor $\frac{\varphi}{\sin\varphi}$.
This is exactly the ratio of the circumference of a circle in the plane
${\mathfrak{m}}$ with radius $\varphi$ and the circumference of the constant
latitude circle on the sphere at angle $\varphi$ from the north pole.

The exponential $\operatorname{Exp}(v)$ can be found by solving
${y^{\prime}(t)=\xi_{v}(y(t))}$, $y(0)=o$ at $t=1$, giving
$\operatorname{Exp}(v)=\exp(\hat{v})o$. We can compute the matrix exponential
$\exp(\hat{v})$ via Rodrigues formula [9]. From
$\hat{v}^{3}=-\varphi^{2}\hat{v}$ for $\varphi=(v^{T}v)^{\frac{1}{2}}$ we find

$\exp(\hat{v})=I+\frac{\sin(\varphi)}{\varphi}\hat{v}+\frac{1}{2}\frac{\sin^{2}(\varphi/2)}{(\varphi/2)^{2}}\hat{v}^{2}.$

This gives

$\operatorname{Exp}(v)=\sin(\varphi)\frac{v}{||v||}+\cos(\varphi)o,$

a rotation of $o$ in the $v$-$o$ plane through the angle $\varphi$.

We obtain the following spherical integrator, where we at each step choose a
new base point $o=y_{\ell}$. The algorithm should never take steps
$\theta_{i}\geq\pi$, because of the singularity in
${\operatorname{dExp}}^{-1}$. If this happens, the step size $h$ must be
reduced.

Result: Evolve $y^{\prime}(t)=F(y(t))$ from $y(0)=y_{0}$ to $y_{n}\approx
y(t\\!=\\!T)$,

where $||y(t)||=1$ and $F(y)^{T}y=0$.

Initialisation: Choose time step $h=T/n$ and $\\{a_{i,j}\\}_{i,j=1}^{r}$,
$\\{b_{j}\\}_{j=1}^{r}$ coefficients of RK.

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

for  _$i=1,\ldots,r$_ do

$\theta_{i}=\sum_{j=1}^{r}a_{i,j}\tilde{K}_{j}$

$\varphi=||\theta_{i}||$

if  _$\varphi==0$_ then

$K_{i}=hF(y_{\ell})$

end if

else

$E_{i}=\operatorname{Exp}(\theta_{i})=\frac{\sin\varphi}{\varphi}\theta_{i}+\cos(\varphi)y_{\ell}$

$s=\operatorname{Exp}(\theta_{i}/2)=(E_{i}+y_{\ell})/||E_{i}+y_{\ell}||$

$v=h\ F\big{(}E_{i}\big{)}$

$K_{i}=\Gamma_{\theta_{i}}^{-1}v=v-2ss^{T}v$

$\tilde{K}_{i}={\operatorname{dExp}}^{-1}_{\theta_{i}}K_{i}=K_{i}+\left(\frac{\varphi}{\sin(\varphi)}-1\right)\left(K_{i}-\frac{\theta_{i}\theta_{i}^{T}}{\varphi^{2}}K_{i}\right)$

end if

end for

$\theta=\sum_{j=1}^{r}b_{j}\tilde{K}_{j}$

$\varphi=||\theta||$

$y_{\ell+1}=\operatorname{Exp}(\theta)=\frac{\sin\varphi}{\varphi}\theta+\cos(\varphi)y_{\ell}$

end for

Algorithm 5 CSI (Canonical Spherical Integrator)

The algorithm is not time symmetric for a time symmetric dynamical system,
even if the underlying RK method is ’self adjoint’ (time symmetric). The
reason for this is the choice of $o=y_{\ell}$ at each step, which is not
symmetric with respect to time reversal. An alternative is to choose the base
point as the geodesic midpoint
$o=(y_{\ell}+y_{\ell+1})/||y_{\ell}+y_{\ell+1}||$. This yields a time
symmetric integrator for time symmetric problems and self adjoint RK methods,
see [19] for such methods in the Lie group case. But, of course, such methods
are necessarily implicit. See remark on relations to the Spherical Midpoint
Method [13] in Section 4.

###  3.3. Hyperbolic spaces

$H^{n}$ is the unique simply connected $n$-dimensional Riemannian manifold of
constant sectional curvature $-1$. We will present this via the “hyperboloid
model”, an isometric embedding of $H^{n}$ in $n+1$ dimensional Minkowski space
[6].

Let $\left\\{\left.{\bf{z}}=(\underline{z};{z}_{0})\in\mathbb{R}^{n+1}\
\right|\ \underline{z}\in\mathbb{R}^{n},{z}_{0}\in\mathbb{R}\right\\}$ be
standard coordinates on Minkowski space, where $(\underline{z};{z}_{0})$
denotes a column vector with ${z}_{0}$ in last position. Define the indefinite
Minkowski inner product

(16) $\langle{\bf y},{\bf
z}\rangle={\bf{y}}^{T}J{\bf{z}}={y}_{0}{z}_{0}-\underline{y}^{T}\underline{z}$

where $J=\operatorname{diag}(-1,\ldots,-1,+1)$. Special relativity is defined
on spacetime $\mathbb{R}^{4}$, where $\underline{z}\in\mathbb{R}^{3}$ is space
and ${z}_{0}\in\mathbb{R}$ is time.

The subset at unit distance from the origin splits in two connected
hyperboloids

$H^{\pm}=\left\\{\left.{\bf{z}}=(\underline{z};{z}_{0})\in\mathbb{R}^{n+1}\
\right|\ \langle{\bf z},{\bf
z}\rangle=1\right\\}=\left\\{\left.(\underline{z};{z}_{0})\ \right|\
{z}_{0}=\pm\sqrt{\underline{z}^{T}\underline{z}+1}\right\\}.$

We identify $H^{n}:=H^{+}$. Let the _Lorentz group_ $O(n,1)$ be the matrix
group preserving the Minkowski inner product:

$\displaystyle O(n,1):=$
$\displaystyle\left\\{\left.A\in\operatorname{GL}(n+1,\mathbb{R})\ \right|\
\langle A{\bf{y}},A{\bf{z}}\rangle=\langle{\bf{y}},{\bf{z}}\rangle\mbox{ for
all ${\bf{y}},{\bf{z}}\in\mathbb{R}^{n+1}$}\right\\}$ $\displaystyle=$
$\displaystyle\left\\{\left.A\in\operatorname{GL}(n+1,\mathbb{R})\ \right|\
JAJ=A^{-T}\right\\}.$

Polar decompositions in Lie groups are derived from involutive automorphisms
on Lie groups [16]. The classical matrix polar decomposition $A=US$, where
$U^{T}U=I$ is orthogonal222We usually write orthogonal matrices ‘$Q$’, but
here we need to distinguish from the quadratic representation. and $S$ is
symmetric positive definite (SPD), is obtained from the involutive
automorphism on $\operatorname{GL}(n+1,\mathbb{R})$ given as
$\alpha(A)=A^{-T}$ by requiring $\alpha(U)=U$ and $\alpha(S)=S^{-1}$. An other
involutive automorphism important for the Lorentz group is $\beta(A)=JAJ$. The
two automorphisms commute, so $\beta\alpha(A)=JA^{-T}J$ is also involutive.
Note that $O(n,1)=\left\\{\left.A\in\operatorname{GL}(n+1,\mathbb{R})\
\right|\ \beta\alpha(A)=A\right\\}$.

######  Proposition 3.1.

The polar decomposition of $A\in O(n,1)$ is given as $A=US_{\underline{s}}$
where

(17) $U=\left(\begin{array}[]{cc}\widetilde{U}&0\\\ 0&u\end{array}\right)$

for an orthogonal $\widetilde{U}\in O(n)$, $u=\pm 1$, and

(18)
$S_{\underline{s}}=\left(\begin{array}[]{cc}\widetilde{S}&\underline{s}\\\
\underline{s}^{T}&{s}_{0}\end{array}\right)$

for some $\underline{s}\in\mathbb{R}^{n}$,
$\widetilde{S}=\sqrt{I_{n}+\underline{s}\underline{s}^{T}}$ (SPD square root)
and ${s}_{0}=\sqrt{1+\underline{s}^{T}\underline{s}}$, thus
$(\underline{s};{s}_{0})\in H^{n}$.

###### Proof.

An alternative and detailed proof is given in [6]. Here we sketch a structural
argument, omitting some minor details. We seek $A=US$ where $\alpha(U)=U$ and
$\alpha(S)=S^{-1}$. We have

$US=A=\beta\alpha(A)=\beta\alpha(U)\beta\alpha(S)=\beta(U)\beta(S^{-1}),$

which yields $\beta(U)=U$ and $\beta(S)=S^{-1}$. From $\beta(U)=U$ follows

$U=\left(\begin{array}[]{cc}\widetilde{U}&0\\\ 0&u\end{array}\right),$

and $\alpha(U)=U$ implies $\widetilde{U}$ is orthogonal and $u^{2}=1$. Writing

$S=\left(\begin{array}[]{cc}\widetilde{S}&\underline{s}\\\
\underline{s}^{T}&{s}_{0}\end{array}\right)$

we find from $\beta(S)=S^{-1}$ that

$S^{-1}=\left(\begin{array}[]{cc}\widetilde{S}&-\underline{s}\\\
-\underline{s}^{T}&{s}_{0}\end{array}\right).$

Multiplying $SS^{-1}=I_{n+1}$ yields
$\widetilde{S}^{2}-\underline{s}\underline{s}^{T}=I_{n}$,
${s}_{0}^{2}-\underline{s}^{T}\underline{s}=1$ and
$\widetilde{S}\underline{s}={s}_{0}\underline{s}$. ∎

Let $O_{0}(n,1)$ denote the matrices in $O(n,1)$ where $u=+1$. This is the
subgroup mapping $H^{+}\mapsto H^{+}$ and $H^{-}\mapsto H^{-}$. The matrices
where $u=-1$ swap the two components of $H^{\pm}$, but they do not form a
subgroup. Note that

$\det(S^{2})=\det(SJS^{-1}J)=\det(SS^{-1})=\det(I)=1,$

and since $S$ is SPD, we have $\det(S)=1$. Hence
$\det(A)=\det(U)=u\det(\widetilde{U})=\pm 1$. The special Lorentz group

$\operatorname{SO}(n,1):=\left\\{\left.A\in O(n,1)\ \right|\
\det(A)=1\right\\}$

has two connected components for $u=\pm 1$. We denote
$\operatorname{SO}_{0}(n,1)$ the connected component containing the identity.
$\operatorname{SO}(n,1)$ consists of those $A$ where $u\det(\widetilde{U})=1$
and $\operatorname{SO}_{0}(n,1)$ those where $u=\det(\widetilde{U})=1$.

We will present $H^{n}$ as a homogeneous space. Let $o=({\underline{0}};1)\in
H^{n}$ be the base point. Since $So=(\underline{s};{s}_{0})$, it follows that
both $O_{0}(n,1)$ and $\operatorname{SO}_{0}(n,1)$ act transitively on
$H^{n}$. The stabiliser subgroup at $o$ consists of $A=US$ such that $Ao=o$.
Multiplied out this implies $S=I_{n+1}$ and

$U=\left(\begin{array}[]{cc}\widetilde{U}&\underline{0}\\\
\underline{0}&1\end{array}\right),\quad\widetilde{U}\in O(n).$

Thus

$H^{n}=O_{0}(n,1)/O(n)\simeq\operatorname{SO}_{0}(n,1)/\operatorname{SO}(n).$

We seek the symmetric product on $H^{n}$. Note that $\sigma_{o}:=J\in
O_{0}(n,1)$ is the isometric reflection around $o$. Reflection around a
general point ${\bf{s}}=(\underline{s};{s}_{0})\in H^{n}$ is obtained by
moving the point to $o$, reflect with $J$ and move back,

$\sigma_{{\bf{s}}}=S_{\underline{s}}\sigma_{o}S_{\underline{s}}^{-1}=\left(\begin{array}[]{cc}C&\underline{s}\\\
\underline{s}^{T}&c\end{array}\right)\left(\begin{array}[]{cc}-I&{\bf 0}\\\
{\bf 0}^{T}&1\end{array}\right)\left(\begin{array}[]{cc}C&-\underline{s}\\\
-\underline{s}^{T}&c\end{array}\right)=\left(\begin{array}[]{cc}-I-2\underline{s}\underline{s}^{T}&2{s}_{0}\underline{s}\\\
-2{s}_{0}\underline{s}^{T}&1+2\underline{s}^{T}\underline{s}\end{array}\right).$

From this we find the more ‘obvious’ expression for the isometric reflection

(19) $\sigma_{{\bf{s}}}=2{\bf{s}}\langle{\bf{s}},\\_\rangle-I.$

The quadratic representation is

(20)
$Q_{{\bf{s}}}=\sigma_{{\bf{s}}}\sigma_{o}=S_{\underline{s}}JS_{\underline{s}}^{-1}J=S_{\underline{s}}^{2}=S_{2{s}_{0}\underline{s}}=\left(\begin{array}[]{cc}I_{n}+2\underline{s}\underline{s}^{T}&2{s}_{0}\underline{s}\\\
2{s}_{0}\underline{s}^{T}&1+2\underline{s}^{T}\underline{s}\end{array}\right).$

The computation of geodesic exponential $\operatorname{Exp}$, triple bracket
and $\operatorname{dExp}^{-1}$ is very similar to the spherical case, Section
3.2. We compute the geodesic exponential from (8)-(9). From the embedding
$H^{n}\subset\mathbb{R}^{n+1}$ we see that
${\mathfrak{m}}=T_{o}H^{n}=\mathbb{R}^{n}\simeq\\{{\bf{v}}=(v;0)\in\mathbb{R}^{n+1}\\}$.
Let $\gamma(t)=o+\frac{1}{2}t{\bf{v}}+\mathcal{O}(t^{2})$. At ${\bf{y}}\in
H^{n}$ we have

(21) $\xi_{v}({\bf{y}})=\left.\frac{\partial}{\partial
t}\right|_{t=0}Q_{\gamma{t}}{\bf{y}}=\left.\frac{\partial}{\partial
t}\right|_{t=0}\left(\begin{array}[]{cc}I_{n}&tv\\\
tv^{T}&1\end{array}\right){\bf{y}}=\left(\begin{array}[]{cc}{\bf 0}&v\\\
v^{T}&0\end{array}\right){\bf{y}}=\hat{{\bf{v}}}{\bf{y}},$

where the hat-map in this case is $\hat{{\bf{v}}}={\bf{v}}\langle
o,\\_\rangle-o\langle{\bf{v}},\\_\rangle$.

_Remark_ : The structure of this hat map can be explained from the involutive
automorphisms $\alpha,\beta\colon GL(n+1,\mathbb{R})\rightarrow
GL(n+1,\mathbb{R})$. Differentiating these we get involutive automorphisms
$d\alpha,d\beta\colon{\mathfrak{g}l}(n+1,\mathbb{R})\rightarrow{\mathfrak{g}l}(n+1,\mathbb{R})$

$d\alpha(A)=-A^{T},\qquad d\beta(A)=JAJ.$

Since $\alpha(S)=S^{-1}$ and $\beta(S)=S^{-1}$, we have $S=\exp(V)$ for
$d\alpha(V)=-V$ and $d\beta(V)=-V$. This implies that $V$ is of the form
$V=\hat{{\bf{v}}}$. Let $d\alpha^{-}$ and $d\beta^{-}$ denote the -1
eigenspaces of these involutive automorphisms on
${\mathfrak{g}l}(n+1,\mathbb{R})$, then
$\hat{{\mathfrak{m}}}:={\mathfrak{g}l}(n+1,\mathbb{R})\cap d\alpha^{-}\cap
d\beta^{-}$ is a Lie triple system (LTS) of symmetric matrices and $\hat{\
}\colon{\mathfrak{m}}\rightarrow\hat{{\mathfrak{m}}}$ is the unique LTS
isomorphism such that ${\bf{v}}=\hat{{\bf{v}}}o$. A matrix of the form of $S$
is called a _Lorentz boost_ in special relativity, and is the exponential of a
matrix in $\hat{{\mathfrak{m}}}$.

From the hat map isomorphism we compute the triple bracket and the double
adjoint on ${\mathfrak{m}}$:

$\displaystyle[u,v,w]$
$\displaystyle=[[\hat{u},\hat{v}],\hat{w}]o=\big{(}v^{T}wI-vw^{T}\big{)}u$
$\displaystyle{\operatorname{ad}}^{2}_{v}w$
$\displaystyle=[w,v,v]=\big{(}v^{T}vI-vv^{T}\big{)}w.$

Note that this is the negative of the triple bracket in the spherical case,
which is not strange since the triple bracket is the Riemannian curvature
tensor.

For $x={\operatorname{ad}}^{2}v$ we find $x^{2}=\varphi^{2}x$ for
$\varphi=\sqrt{-\langle{\bf{v}},{\bf{v}}\rangle}$. As in the spherical case,
we find

(22)
${\operatorname{dExp}}^{-1}_{{\bf{v}}}=\frac{\sqrt{x}}{\sinh(x)}=I+\left(\frac{\varphi}{\sinh(\varphi)}-1\right)\frac{x}{\varphi^{2}}=I+\left(\frac{\varphi}{\sinh(\varphi)}-1\right)\pi_{{\bf{v}}}^{\perp},$

where

(23)
$\pi_{\bf{v}}^{\perp}=I-\frac{\langle{\bf{v}},\\_\rangle}{\langle{\bf{v}},{\bf{v}}\rangle}$

is the Minkowski projection on the orthogonal complement of ${\bf{v}}$.

From (21) we find the the geodesic exponential
$\operatorname{Exp}\colon{\mathfrak{m}}\rightarrow H^{n}$

(24) $\operatorname{Exp}({\bf{v}})=\exp\left(\hat{{\bf{v}}}\right)o.$

Similar to the spherical case, $\hat{{\bf{v}}}^{3}=\varphi^{2}\hat{{\bf{v}}}$
yields a Rodrigues type formula

(25)
$\exp(\hat{{\bf{v}}})=I+\frac{\sinh(\varphi)}{\varphi}\hat{{\bf{v}}}+\frac{1}{2}\frac{\sinh^{2}(\varphi/2)}{(\varphi/2)^{2}}\hat{{\bf{v}}}^{2}.$

This gives the geodesic exponential

(26)
$\operatorname{Exp}({\bf{v}})=\frac{\sinh(\varphi)}{\varphi}{\bf{v}}+\cosh(\varphi)o,$

for $\langle{\bf{v}},o\rangle=0$. This is a hyperbolic rotation of the
${\bf{v}}$-$o$ plane through the hyperbolic angle $\varphi$.

Similar to CSI Algorithm 5 we will express our hyperbolic integrator using at
each step the base point $o={\bf{y}}_{n}$, where each step is performed as in
Algorithm 1. At an arbitrary point $o={\bf{y}}_{n}\in H^{n}$ the tangent space
is

${\mathfrak{m}}=T_{o}H^{n}=\left\\{\left.{\bf{w}}\in\mathbb{R}^{n+1}\ \right|\
\langle{\bf{w}},o\rangle=0\right\\}.$

The geodesic exponential is still given by (26), with inverse differential
(22)-(23). Finally, for $\theta\in{\mathfrak{m}}$, we need
$\Gamma_{\theta}^{-1}$, parallel transport from
$T_{\operatorname{Exp}(\theta)}H^{n}$ to ${\mathfrak{m}}$ along the joining
geodesic. Let ${\bf{s}}=\operatorname{Exp}(\theta/2)$ be the midpoint on the
geodesic. From Definition 2.3 we find

$\Gamma_{\theta}^{-1}=TQ_{{\bf{s}}}^{-1}=T\left(\sigma_{{\bf{s}}}\sigma_{o}\right)^{-1}=T\left(\sigma_{o}\sigma_{{\bf{s}}}\right)=-T\sigma_{{\bf{s}}}$

since $T\sigma_{o}=-I$. From (19) this yields

(27) $\Gamma_{\theta}^{-1}=I-2{\bf{s}}\langle{\bf{s}},\\_\rangle.$

Result: Evolve ${\bf{y}}^{\prime}(t)=F({\bf{y}}(t))$, ${\bf{y}}(t)\in H^{n}$,
from ${\bf{y}}(0)={\bf{y}}_{0}$ to ${\bf{y}}_{n}\approx{\bf{y}}(t\\!=\\!T)$,

where $\langle{\bf{y}}(t),{\bf{y}}(t)\rangle=1$ and $\langle
F({\bf{y}}),{\bf{y}}\rangle=0$.

Initialisation: Choose time step $h=T/n$ and $\\{a_{i,j}\\}_{i,j=1}^{r}$,
$\\{b_{j}\\}_{j=1}^{r}$ coefficients of RK.

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

for  _$i=1,\ldots,r$_ do

$\theta_{i}=\sum_{j=1}^{r}a_{i,j}\tilde{K}_{j}\in{\mathfrak{m}}=T_{{\bf{y}}_{\ell}}H^{n}$

$\varphi=\sqrt{-\langle\theta_{i},\theta_{i}\rangle}$

if  _$\varphi==0$_ then

$K_{i}=hF({\bf{y}}_{\ell})$

end if

else

${\bf{u}}=\operatorname{Exp}(\theta_{i})=\frac{\sinh(\varphi)}{\varphi}\theta_{i}+\cosh(\varphi){\bf{y}}_{\ell}$

${\bf{s}}=\operatorname{Exp}(\theta_{i}/2)=\frac{\sinh(\varphi/2)}{\varphi}\theta_{i}+\cosh(\varphi/2){\bf{y}}_{\ell}$

${\bf{v}}=h\ F\big{(}{\bf{u}}\big{)}$

$K_{i}=\Gamma_{\theta_{i}}^{-1}{\bf{v}}={\bf{v}}-2{\bf{s}}\langle{\bf{s}},{\bf{v}}\rangle$

$\tilde{K}_{i}={\operatorname{dExp}}^{-1}_{\theta_{i}}K_{i}=K_{i}+\left(\frac{\varphi}{\sinh(\varphi)}-1\right)\left(K_{i}-\frac{\langle\theta_{i},K_{i}\rangle}{\varphi^{2}}\theta_{i}\right)$

end if

end for

$\theta=\sum_{j=1}^{r}b_{j}\tilde{K}_{j}$

$\varphi=\sqrt{-\langle\theta,\theta\rangle}$

${\bf{y}}_{\ell+1}=\operatorname{Exp}(\theta)=\frac{\sinh\varphi}{\varphi}\theta+\cosh(\varphi){\bf{y}}_{\ell}$

end for

Algorithm 6 CHI (Canonical Hyperbolic Integrator)

####  3.3.1. Other models of $H^{n}$

There are other well known models of $H^{n}$, most famous are the Poincaré
_half space_ and the _disc models_. These are obtained from the hyperboloid
model by stereographic projections. The half space model is the stereographic
projection from infinity along the edge of the light cone in the $x_{n}$ \-
$x_{0}$ plane, where the hyperbolic space is realised as the upper half space
of $\mathbb{R}^{n}$ given by $x_{n}>0$, or upper half plane of $\mathbb{C}$ in
the case $n=2$. The disc model is the stereographic projection of $H^{n}$ from
the point $(\underline{0};-1)$, realising the hyperbolic space as the unit
ball in $\mathbb{R}^{n}$, or unit disc $\mathbb{C}$ in the case $n=2$. These
models are beautiful, in particular for $n=2$ where the hyperbolic space
becomes a homogeneous space under the action of the projective linear
fractional transformations $\operatorname{PSL}(2,\mathbb{R})$ on $\mathbb{C}$.
However, these models are mathematically equivalent to the model we have
presented above, so we omit the details. Spaces of constant negative curvature
are important in the geometric theory of dynamical systems. E.g. Anosov flows
on compact surfaces $H^{2}/\Gamma$, where $\Gamma$ is a discrete subgroup of
$\operatorname{PSL}(2,\mathbb{R})$ (Fuchsian group) are generic examples of
Axiom A type dynamical systems with chaotic solutions. Variants of the CHI
algorithm should be interesting for numerical studies of Anosov flows. In
numerical analysis, Riemannian spaces of negative curvature have recently
appeared in geometric generalisations of classical stability theory to
Riemannian manifolds [1]. The present algorithms should be of interest in the
development of these ideas.

##  4\. Final remarks

We have in this paper developed the general concept of canonical geometric
integration on symmetric spaces. The important issue of numerical qualities of
these algorithms and their applications to practical computational problems
are left to sequel papers. Here we briefly remark on interesting questions to
be pursued, and possible applications.

Differential equations evolving on spheres are ubiquitous, with important
examples from robotics, computational mechanics, rigid body dynamics and flows
on planetary surfaces. Infinite dimensional spheres are natural domains for
partial differential equations preserving energy or $L_{2}$ norms, such as the
Schrödinger and Korteweg–de Vries equations. It is interesting that the CSI
algorithm on the $n$-sphere is not significantly more expensive than classical
RK methods on the embedding of the sphere in $\mathbb{R}^{n+1}$.

Hyperbolic geometry is the foundation of special relativity. Our CHI algorithm
is formulated in terms of Lorentz transformations and should be interesting
for geometric integration of Maxwell equations and dynamical systems in
special relativity, as it preserves the Lorentz invariance of the equations.

An other interesting class of equations evolving on symmetric spaces are Lie-
Poisson systems. Hamiltonian mechanics are naturally formulated on cotangent
bundles of Lie groups. By symmetry reduction, these can often be reduced to
Lie-Poisson systems evolving on the dual of Lie algebras. The natural
coadjoint action foliates the dual Lie algebras in symplectic leafs, which in
many cases are symmetric spaces. For example spheres are symplectic leaves for
${\mathfrak{s}o}(n)^{*}$ and hyperbolic spaces for ${\mathfrak{s}l}(n)^{*}$.
Questions of symplectic algorithms for such problems are interesting, hard and
largely open. We remark that the CSI integrator based on the implicit midpoint
RK and the geodesic midpoint as base point for each step is _not_ the same as
the celebrated symplectic spherical midpoint method [13]. However the
difference between these two methods is closely related to the difference
between using the exponential and the Cayley transform as coordinate mappings.
For use of Cayley maps, we refer to the pioneering work on symplectic
integration on spheres by Lewis and Simo [11]. Symplectic integration on
hyperbolic symplectic leaves is to our knowledge largely unchartered
territory.

Differential equations evolving on the space of symmetric positive definite
matrices occur in different contexts. One example is in the inverse eigenvalue
problem for SPD Toeplitz matrices, formulated as isospectral flows [4]. A
different example is the tracing of nerve fibres in diffusion tensor imaging
of the brain, where the voxels are symmetric positive definite diffusion
tensors.

We have in this paper not written out the example of Grassman manifolds, which
is yet another example of symmetric spaces. Differential equations evolving on
Grassman manifolds has been considered by several authors, see [2] and
references within.

## References

  * [1] Martin Arnold, Elena Celledoni, Ergys Çokaj, Brynjulf Owren, and Denise Tumiotto.  B-stability of numerical integrators on Riemannian manifolds.  arXiv preprint arXiv:2308.08261, 2023.
  * [2] Margaret Beck, Anastasia Doikou, Simon JA Malham, and Ioannis Stylianidis.  Grassmannian flows and applications to nonlinear partial differential equations.  In Computation and Combinatorics in Dynamics, Stochastics and Control: The Abel Symposium, Rosendal, Norway, August 2016, pages 71–98. Springer, 2018.
  * [3] Geir Bogfjellmo, Klas Modin, and Olivier Verdier.  A numerical algorithm for C^2-splines on symmetric spaces.  SIAM Journal on Numerical Analysis, 56(4):2623–2647, 2018.
  * [4] Mari Calvo, Arieh Iserles, and Antonella Zanna.  Numerical solution of isospectral flows.  Mathematics of computation, 66(220):1461–1486, 1997.
  * [5] M Camarinha, F Silva Leite, and P Crouch.  Splines of class Ck on non-Euclidean spaces.  IMA Journal of Mathematical Control and Information, 12(4):399–410, 1995.
  * [6] Jean Quaintance Gallier and Jocelyn Quaintance.  Differential geometry and Lie groups, volume 12.  Springer, 2020.
  * [7] Sigurdur Helgason.  Differential geometry and symmetric spaces, volume 341.  American Mathematical Soc., 2001.
  * [8] Nicholas J Higham and Robert S Schreiber.  Fast polar decomposition of an arbitrary matrix.  SIAM Journal on Scientific and Statistical Computing, 11(4):648–655, 1990.
  * [9] Arieh Iserles, Hans Z Munthe-Kaas, Syvert P Nørsett, and Antonella Zanna.  Lie-group methods.  Acta numerica, 9:215–365, 2000.
  * [10] S Kobayashi and K Nomizu.  Foundations of differential geometry I, Interscience.  New York, 1963.
  * [11] Debra Lewis and Juan Carlos Simo.  Conserving algorithms for the dynamics of Hamiltonian systems on Lie groups.  Journal of Nonlinear Science, 4:253–299, 1994.
  * [12] Ottmar Loos.  Symmetric spaces: General theory, volume 1.  WA Benjamin, 1969.
  * [13] Robert McLachlan, Klas Modin, and Olivier Verdier.  A minimal-variable symplectic integrator on spheres.  Mathematics of Computation, 86(307):2325–2344, 2017.
  * [14] Hans Munthe-Kaas.  Lie–Butcher theory for Runge-Kutta methods.  BIT Numerical Mathematics, 35:572–587, 1995.
  * [15] Hans Munthe-Kaas.  High order Runge-Kutta methods on manifolds.  Applied Numerical Mathematics, 29(1):115–127, 1999.
  * [16] Hans Z Munthe-Kaas, Gilles Reinout W Quispel, and Antonella Zanna.  Symmetric spaces and Lie triple systems in numerical analysis of differential equations.  BIT Numerical Mathematics, 54:257–282, 2014.
  * [17] Hans Z Munthe-Kaas, GRW Quispel, and Antonella Zanna.  Generalized polar decompositions on Lie groups with involutive automorphisms.  Foundations of Computational Mathematics, 1:297–324, 2001.
  * [18] Hans Z Munthe-Kaas and Jonatan Stava.  Lie Admissible Triple Algebras: The Connection algebra of symmetric spaces.  arXiv preprint arXiv:2306.15582, 2023.
  * [19] Antonella Zanna, Kenth Engø, and Hans Z Munthe-Kaas.  Adjoint and selfadjoint Lie-group methods.  BIT Numerical Mathematics, 41:395–421, 2001.