A Matrix Valued Kuramoto Model

TL;DR

This paper introduces a new matrix-valued Kuramoto model where symmetric matrices evolve to align their eigenframes, analyzing stability and dynamics, including the instability of twist states.

Contribution

It proposes a novel matrix-valued Kuramoto model and provides stability analysis for phase-locked and twist states, extending synchronization theory to non-commuting matrices.

Findings

Phase-locked state is stable when matrices align.

Twist states are dynamically unstable.

Eigenframe alignment leads to commuting matrices.

Abstract

Beginning with the work of Lohe [14,15] there have been a number of papers [3,5,8,9,11] that have generalized the Kuramoto model for phase-locking to a non-commuting situation. Here we propose and analyze another such model. We consider a collection of symmetric matrix-valued variables that evolve in such a way as to try to align their eigenvector frames. The phase-locked state is one where the eigenframes all align, and thus the matrices all commute. We analyze the stability of the phase-locked state and show that it is stable. We also analyze a dynamic analog of the twist states arising in the standard Kuramoto model, and show that these twist states are dynamically unstable.

Equations

$$\frac{d{\theta }_i}{dt}={\omega }_i+\gamma \sum _i\sin ({\theta }_j-{\theta }_i),$$

$$i\dot{{\mathbf{U}}_j}{\mathbf{U}}_j^{\ast }={\mathbf{H}}_j-\frac{iK}{2N}\sum _{k=1}^N{c}_{kj}({\mathbf{U}}_j{\mathbf{U}}_k^{\ast }-{\mathbf{U}}_k{\mathbf{U}}_j^{\ast }),$$

$$\frac{d}{dt}{\mathbf{X}}_i{\dot{\mathbf{X}}}_i^{-1}={\mathbf{\Omega }}_i+\frac{1}{2}\sum _{j=1}^n{\gamma }_{ij}\left(f({\mathbf{X}}_j{\mathbf{X}}_i^{-1})-f({\mathbf{X}}_i{\mathbf{X}}_j^{-1})\right),$$

$$i\hslash \frac{\partial }{\partial t}|{\psi }_i⟩={\mathbf{H}}_i|{\psi }_i⟩+\frac{i\hslash \kappa }{2N}\sum _j{a}_{ij}\left(|{\psi }_j⟩⟨{\psi }_i|{\psi }_i⟩-|{\psi }_i⟩⟨{\psi }_j|{\psi }_i⟩\right).$$

$$u_t=\mathbf{J}{\delta }_{u}H-{\delta }_{u}E-\sum _i{\mu }_i{\delta }_u{C}_i.$$

$$E_d=\frac{\kappa }{4N}\sum _{i,j}⟨{\psi }_i-{\psi }_j|{\psi }_i-{\psi }_j⟩$$

$${\tilde{E}}_d=\frac{\kappa }{4N}\sum _{i,j}⟨{\psi }_i-{\psi }_j|{\psi }_i-{\psi }_j⟩-\sum _{i=1}^N{\mu }_i⟨{\psi }_i|{\psi }_i⟩.$$

$$\frac{\partial }{\partial t}|{\psi }_i⟩=J\frac{\partial }{\partial ⟨{\psi }_i|}H-\frac{\partial }{\partial ⟨{\psi }_i|}{\tilde{E}}_d$$

$$i\hslash \frac{\partial }{\partial t}|{\psi }_i⟩={\mathbf{H}}_i|{\psi }_i⟩+\frac{i\hslash \kappa }{2N}\sum _j{a}_{ij}|{\psi }_j-{\psi }_i⟩-{\mu }_i|{\psi }_i⟩.$$

$$i\hslash \frac{\partial }{\partial t}|{\psi }_i⟩={\mathbf{H}}_i|{\psi }_i⟩+\frac{i\hslash \kappa }{2N}\sum _j{a}_{ij}\left(|{\psi }_j⟩-|{\psi }_i⟩\frac{⟨{\psi }_j|{\psi }_i⟩}{⟨{\psi }_i|{\psi }_i⟩}\right)$$

$$\frac{d{\mathbf{x}}_i}{dt}={\Omega }_i{\mathbf{x}}_i+\frac{\kappa }{N}\sum _j{a}_{ij}({\mathbf{x}}_j-{\mathbf{x}}_i)-{\mu }_i{\mathbf{x}}_i.$$

$$E=\frac{1}{4}\sum _{i,j}{a}_{ij}\mathrm{Tr}\left([{\mathbf{M}}_i,{\mathbf{M}}_j{]}^2\right).$$

$$\frac{d{\mathbf{M}}_i}{dt}=-\frac{\partial E}{\partial {\mathbf{M}}_i}.$$

$$\frac{d{\mathbf{M}}_i}{dt}=-\sum _j{a}_{ij}[{\mathbf{M}}_j,[{\mathbf{M}}_i,{\mathbf{M}}_j]].$$

$$\frac{d{\mathbf{M}}_i}{dt}=[{\mathbf{\Omega }}_i,{\mathbf{M}}_i]-\sum _j{a}_{ij}[{\mathbf{M}}_j,[{\mathbf{M}}_i,{\mathbf{M}}_j]]-{\mu }_i{\mathbf{M}}_i.$$

$$\frac{d{\tilde{\mathbf{M}}}_i}{dt}=[{\mathbf{\Omega }}_i-\tilde{\mathbf{\Omega }},{\tilde{\mathbf{M}}}_i]-\sum _j{a}_{ij}[{\tilde{\mathbf{M}}}_j,[{\tilde{\mathbf{M}}}_i,{\tilde{\mathbf{M}}}_j]]-{\mu }_i{\tilde{\mathbf{M}}}_i.$$

$${\mathbf{M}}_i=\mathbf{R}({\theta }_i){\mathbf{\Lambda }}_i\mathbf{R}(-{\theta }_i)=\left[\begin{array}{cc}\cos ({\theta }_i)& \sin ({\theta }_i)\\ -\sin ({\theta }_i)& \cos ({\theta }_i)\end{array}\right]\left[\begin{array}{cc}{\lambda }_i^1& 0\\ 0& {\lambda }_i^2\end{array}\right]\left[\begin{array}{cc}\cos ({\theta }_i)& -\sin ({\theta }_i)\\ \sin ({\theta }_i)& \cos ({\theta }_i)\end{array}\right].$$

$${\mathbf{M}}_i=\frac{1}{2}({\lambda }_i^1+{\lambda }_i^2)\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]+\frac{1}{2}({\lambda }_i^2-{\lambda }_i^1)\left[\begin{array}{cc}-\cos (2{\theta }_i)& \sin (2{\theta }_i)\\ \sin (2{\theta }_i)& \cos (2{\theta }_i)\end{array}\right].$$

$$\frac{d{\mathbf{M}}_i}{dt}=\frac{1}{2}(\dot{{\lambda }_i^1}+\dot{{\lambda }_i^2})\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]+\frac{1}{2}(\dot{{\lambda }_i^2}-\dot{{\lambda }_i^1})\left[\begin{array}{cc}-\cos (2{\theta }_i)& \sin (2{\theta }_i)\\ \sin (2{\theta }_i)& \cos (2{\theta }_i)\end{array}\right]+({\lambda }_i^2-{\lambda }_i^1)\dot{{\theta }_i}\left[\begin{array}{cc}\sin (2{\theta }_i)& \cos (2{\theta }_i)\\ \cos (2{\theta }_i)& -\sin (2{\theta }_i)\end{array}\right].$$

$$[{\mathbf{M}}_j,[{\mathbf{M}}_i,{\mathbf{M}}_j]]=\frac{1}{2}({\lambda }_j^1-{\lambda }_j^2{)}^2({\lambda }_i^1-{\lambda }_i^2)\sin 2({\theta }_i-{\theta }_j)\left[\begin{array}{cc}\sin (2{\theta }_j)& \cos (2{\theta }_j)\\ \cos (2{\theta }_j)& -\sin (2{\theta }_j)\end{array}\right].$$

$$\dot{{\lambda }_i^1}+\dot{{\lambda }_i^2}=0$$

$$\dot{{\lambda }_i^2}-\dot{{\lambda }_i^1}=-({\lambda }_i^2-{\lambda }_i^1)\sum _{j=1}^N({\lambda }_j^2-{\lambda }_j^1{)}^2{\sin }^2(2({\theta }_j-{\theta }_i))$$

$$\dot{{\theta }_i}=\sum _{j=1}^N\frac{1}{4}({\lambda }_j^2-{\lambda }_j^1{)}^2\sin (4({\theta }_j-{\theta }_i))$$

$$\dot{{\lambda }_i^1}+\dot{{\lambda }_i^2}=-{\mu }_i({\lambda }_i^1+{\lambda }_i^2)$$

$$\dot{{\lambda }_i^2}-\dot{{\lambda }_i^1}=-({\lambda }_i^2-{\lambda }_i^1)\sum _{j=1}^N({\lambda }_j^2-{\lambda }_j^1{)}^2{\sin }^2(2({\theta }_j-{\theta }_i))-\frac{1}{2}{\mu }_i({\lambda }_i^2-{\lambda }_i^1)$$

$$\dot{{\theta }_i}={\omega }_i+\sum _{j=1}^N\frac{1}{4}({\lambda }_j^2-{\lambda }_j^1{)}^2\sin (4({\theta }_j-{\theta }_i))$$

$$\frac{d{\theta }_j}{dt}=f_j({\theta }_1,{\theta }_2,\dots ,{\theta }_N)$$

$$\frac{d\alpha }{dt}$$

$$=-{\alpha }^3\sum _{k=1}^N\frac{1}{2}\left[1-\cos \left(4\left(\frac{\pi n(k-1)}{2N}-\frac{\pi n(j-1)}{2N}\right)\right)\right]$$

$$=-{\alpha }^3\sum _{k=1}^N\frac{1}{2}\left[1-\cos \left(\frac{2\pi n(k-1)}{2N}-\frac{2\pi n(j-1)}{2N}\right)\right].$$

$$\alpha (t)=\frac{\alpha (0)}{\sqrt{1+N{\alpha }^2(0)t}},$$

Full text

1 Introduction

The Kuramoto model,

[TABLE]

where $\omega_{i}\in\mathbb{R}$ and $\gamma$ is the coupling constant, is a simplified model for many processes involving the sychronization of nonlinear coupled oscillators. This model has been used in a wide range of biological and physical applications, including modeling the phenomena of jet lag[16] and circadian rhythms[2], and describing electrical networks[6, 19]. Numerous variations of the Kuramoto model have been developed and studied, including a non-Abelian generalization of the Kuramoto model known as the Lohe model. Lohe actually proposed several related models: the first was a flow on unitary matrices

[TABLE]

where $j=1,...,N$, $c_{kj}\in\mathbb{R}$, $k\geqslant 0$ is the coupling strength, ${\bf U}_{j}$ is a $d\times d$ unitary matrix, ${\bf U}_{j}^{*}$ is the Hermitian conjugate of ${\bf U}_{j}$, and ${\bf H}_{j}$ is a $d\times d$ Hermitian matrix whose eigenvalues correspond to the natural frequencies of the Lohe oscillator at node $j$. This is Equation (2) in Lohe’s 2009 paper [14]. In this work Lohe demonstrated that in the one dimensional case, where $U_{i}=e^{-i\theta_{i}}$, an element of the Abelian unitary group $U(1)$, his model reduces to standard Kuramoto (1.1) model. In the same paper Lohe also proposes several other models based on other matrix groups. Additionally, other non-Abelian variations of the Kuramoto model have been developed which reduce to (1.1) or variants when an Abelian Lie group is chosen, including the quantum Kuramoto model[5]. This model is of the form

[TABLE]

where $f$ is a real analytic function, ${\bf X_{i}}$ is an element of some Lie group (satisfying a certain mild structural condition) and ${\bf\Omega}_{i}$ is an element of the Lie algebra.

In a 2010 paper[15] Lohe proposes a somewhat different model, Equation (10), in the form

[TABLE]

Here each wavefunction $|\psi_{i}\rangle$ evolves according to a Hamiltonian ${\bf H}_{i}$, playing the role of the frequency in the standard Kuramoto model, and the wavefunctions interact through a nonlinear coupling. Our model is most closely related to this model, though the evolution takes place at the level of the operators, not the wavefunctions, so it is worthwhile to spend a little time discussing the structure of these models. In particular these models take the general form of a conservative (Hamiltonian) part plus a constrained gradient descent.

[TABLE]

Here $\delta_{u}$ represents the gradient in the finite dimensional case, and the Euler-Lagrange operator in the infinite dimensional case, ${\bf J}$ is some skew-adjoint operator representing the Hamiltonian structure, $E$ is some functional to be minimized, attainment of the minimum of $E$ representing synchronization, $C_{i}$ representing some functionals to be conserved, and $\mu_{i}$ representing Lagrange multipliers to enforce the constraints that $C_{i}(u)$ be constant. The Hamiltonian flow is, obviously, the ${\bf H}_{i}|\psi_{i}\rangle$ term in Equation (1.4). The remaining terms are a constrained gradient flow: these terms act to decrease the Dirichlet-like energy

[TABLE]

subject to the $N$ constraints that the wavefunctions $|\psi_{i}\rangle$ be properly normalized: $\langle\psi_{i}|\psi_{i}\rangle=\|\psi_{i}\|^{2}=1$. To see this we note that the gradient flow on the functional

[TABLE]

Here the first term, the Dirichlet energy, measures the degree of alignment of the wavefuntions $|\psi_{i}\rangle$, while the scalar quantities $\mu_{i}$ are Lagrange multipliers to maintain the constraint $\langle\psi_{i}|\psi_{i}\rangle=1$. Physically this can be interpreted as conserving the particle number of the $i^{th}$ species, and so $\mu_{i}$ is essentially the chemical potential associated with species $i$.

If we take the dynamics to be

[TABLE]

with Hamiltonian $H=\sum_{i}\langle\psi_{i}|H_{i}|\psi_{i}\rangle$ and symplectic form $J=-\frac{i}{\hbar}$ then we find the dynamics

[TABLE]

The Lagrange multipliers $\mu_{i}$ are determined by the condition that the particle numbers of individual species $\langle\psi_{i}|\psi_{i}\rangle$ be conserved: $\frac{\partial}{\partial t}\langle\psi_{i}|\psi_{i}\rangle=0$. This gives

[TABLE]

which is equivalent to (1.4) under a trivial rescaling of $|\psi_{i}\rangle$. Similarly to derive Equation (20) in Lohe’s 2009 paper we take $u=\bigoplus_{i=1}^{N}{\bf x}_{i}$, $H=\frac{1}{2}\sum_{i}\|{\bf x}_{i}\|^{2}=\|u\|^{2}$, $E=\frac{\kappa}{4N}\sum_{i,j}a_{ij}\|{\bf x}_{i}-{\bf x}_{j}\|^{2}$, ${\bf J}=\bigoplus_{i=1}^{N}{\bf\Omega}_{i}$ and the constraints $C_{i}$ to be $C_{i}({\bf x}_{i})=\frac{1}{2}\|{\bf x}_{i}\|^{2}=\text{constant}$. Substituting these into Equation (1.5) gives

[TABLE]

In this spirit we consider a different generalization of the Kuramoto flow, where the phase space is $N$ $k\times k$ real symmetric matrices. We assume that these matrices evolve according to a gradient flow, where the energy functional is given by the sum of the Hilbert-Schmidt norms of the pairwise matrix commutators. Since this is a gradient flow matrices evolve to minimize energy, whose global minimum is easily seen to be when the matrices pairwise commute.

Throughout this paper all matrices ${\bf M}_{i}$ are assumed to be real symmetric, while all matrices ${\bf\Omega}_{i}$ are assumed to be skew-symmetric. In the remainder of Section 1 we further motivate the current model. We begin by deriving the matrix analog of the Kuramoto model with all frequencies equal. We finish the introduction by establishing a concrete connection to a variant of the standard Kuramoto model. In particular, in Section 1.2 we demonstrate that for the case of $2\times 2$ matrices, this system can be described by a set of equations which is similar to the classical Kuramoto model with Hebbian interactions. In Section 2 we show that in the $2\times 2$ case, this model has a set of unstable fixed points associated with the case in which the matrices $\{{\bf M}_{1},{\bf M}_{2}\ldots{\bf M}_{n}\}$ have eigenvectors pointing in directions which are evenly spaced about the unit circle. This is related to twist states in the classical Kuramoto model. Lastly, in Section 3 we show that for the set of pairwise commuting matrices of any size, this matrix-valued model has conditionally stable fixed points.

1.1 Equations of motion and basic structure

As argued earlier all of the models proposed by Lohe take the same abstract form: they consist of a conservative (Hamiltonian) part representing some abstract rotation, together with a constrained gradient flow part. We follow this same philosophy to derive a non-commuting symmetric matrix valued flow. We first define a gradient flow on a set of symmetric matrices. We will let $[\cdot,\cdot]$ denote the usual matrix commutator $[{\bf A},{\bf B}]={\bf A}{\bf B}-{\bf B}{\bf A}$ and recall that for real symmetric matrices the Hilbert-Schmidt norm is $\|{\bf A}\|^{2}=\operatorname{Tr}{\bf A}^{2}$. We now define an energy function $E$ to be a weighted sum of the Hilbert-Schmidt norms of all pair-wise commutators,

[TABLE]

If we assume that $a_{ij}>0$ and symmetric then $E\geqslant 0$, with $E=0$ if and only if the matrices $\{{\bf M}_{j}\}_{j=1}^{N}$ all commute. The dynamics of the matrices ${\bf M}_{i}$ will be gradient flow on the function $E$.

[TABLE]

The notation here is that if $E$ is a scalar valued function of a $k\times k$ matrix ${\bf M}$ then $\frac{\partial f}{\partial{\bf M}}$ is a $k\times k$ matrix whose entries are the derivatives of $E$ with respect to the corresponding entry of ${\bf M}$, $\left(\frac{\partial f}{\partial{\bf M}}\right)_{ij}=\frac{\partial f}{\partial M_{ij}}.$ A straightforward calculation – see Appendix A for details – gives the equations of motion as

[TABLE]

These equations are not obviously related to the standard Kuramoto model, being cubic in the matrix entries instead of trigonometric, but we next show that, at least in the case of $2\times 2$ matrices the equations can be put into a form that is strongly reminiscent of the classical Kuramoto model. Specifically we show that in terms of the eigenvalues and the angles defining the eigenvectors the flow can be written as a Kuramoto flow for the angles together with a somewhat unusual Hebbian-type evolution for the eigenvalues, which play the role of a coupling strength.

More generally we can consider additional terms to this pure gradient flow representing a conservative/Hamiltonian dynamics and/or constraints. From the point of view of matrices perhaps the most natural assumption is to constrain the Hilbert-Schmidt or Euclidean norms of the operators: $\|{\bf M}_{i}\|^{2}=\operatorname{Tr}({\bf M}_{i}^{2})=\text{\rm constant}$. This leads to a more general system of the form

[TABLE]

where the Lagrange multiplier $\mu_{i},$ given by $\mu_{i}=\frac{\sum_{j}a_{ij}\operatorname{Tr}({\bf M}_{i}[{\bf M}_{j},[{\bf M}_{i},{\bf M}_{j}]]))}{\operatorname{Tr}({\bf M}_{i}^{2})}=\frac{\sum_{j}2a_{ij}\operatorname{Tr}(({\bf M}_{i}{\bf M}_{j})^{2}-({\bf M}_{i}{\bf M}_{j})^{2})}{\operatorname{Tr}{\bf M}_{i}^{2}},$ conserves the Hilbert-Schmidt norm $\operatorname{Tr}({\bf M}_{i}^{2})$ and ${\bf\Omega}_{i}$ is a real skew-symmetric matrix. The conservative term $[{\bf\Omega},{\bf M}_{i}]$ is the standard sort of term one sees in the evolution equations for an operator in the interaction representation.

We note a few invariance properties of Equation (1.7) that generalize the analogous properties of the Kuramoto model. In the classical Kuramoto model the dynamics is invariant under the common rotation $\theta_{i}\mapsto\theta_{i}+\alpha$. Here this extends to orthogonal invariance: if one conjugates all ${\bf M}_{i}$ by the same orthogonal matrix ${\bf O}$, $\tilde{\bf M}_{i}={\bf O}{\bf M}_{i}{\bf O}^{T}$ then the new variables $\tilde{\bf M}_{i}$ satisfy (1.7) as well. Similarly in the classical Kuramoto if one moves to a uniformly rotating coordinate system $\theta_{i}\mapsto\theta_{i}+\bar{\omega}t$ then this is the same as subtracting a mean from all of the frequencies: $\omega_{i}\mapsto\omega_{i}-\bar{\omega}$. Here we have a similar invariance. If we move to a coordinate system that rotates uniformly: $\tilde{\bf M}_{i}=e^{\tilde{\bf\Omega}t}{\bf M}_{i}e^{-\tilde{\bf\Omega}t}$, where $\tilde{\bf\Omega}$ is a skew-symmetric matrix, then (1.7) becomes

[TABLE]

Note that the first invariance in particular has important implications for the spectrum of the linearization about a fixed point. The fact that the dynamics is invariant under conjugation by an element of the orthogonal group ${\bf O}(n)$ implies that the kernel of the linearized operator has dimension at least $\frac{n(n-1)}{2}$. We will discuss this further in a later section.

In most of this paper we will consider ${\bf\Omega}_{i}=0$, so we are considering strict synchronization rather than phase-locking. In our analysis of the twist states the constrained version of the equations arises, with the Lagrange multipliers arising from a change of variable rather than an explicit constraint.

1.2 Spectral coordinates for the $2\times 2$ case: reduction to Kuramoto form.

In order to build intuition we first consider the lowest dimensional case of possible interest, $2\times 2$ real symmetric matrices, under the assumptions that $a_{ij}=1,\mu_{i}=0,{\bf\Omega}_{i}=0$ – no constraints, no external frequency forcing and equally weighted all-to-all coupling. In the case of the Lohe model when one restricts to the Abelian case the dynamics reduces to that of the classical Kuramoto. In our case the dynamics in the $2\times 2$ case is slightly more complicated, and can be expressed as a Kuramoto type evolution for the eigenangle, with the eigenvalues entering as an effective coupling strength between angles, together coupled with an evolution for the eigenvalues.

Let $\{{\bf M}_{i}\}$ be a collection of $N$ symmetric matrices. Using the fact that symmetric matrices are orthogonally diagnonalizable, we write each ${\bf M}_{i}$ as

[TABLE]

Here $\lambda_{i}^{1}$ and $\lambda_{i}^{2}$ represent the two eigenvalues for the $i$-th matrix and $\theta_{i}$ represents the angle of rotation of the orthogonal eigenframe from the standard axes.

We now substitute these representations of ${\bf M}_{i}$ into (1.6). It is somewhat more convenient to work with the equivalent form

[TABLE]

For a fixed $i$, we take the derivative of ${\bf M}_{i}$ and get

[TABLE]

After some algebra we find the following expression for the nested commutator $[{\bf M_{j}},[{\bf M}_{i},{\bf M}_{j}]]$:

[TABLE]

Summing (1.9) over all $j$, yields the representation of the right-hand side of (1.6) under this decomposition. The advantage of this representation is that we have passed the matrix dynamics purely onto the time-dependent eigenvalues and eigenangles. In order to more clearly see this, we again take advantage of the structure of the symmetric matrices. Note that the matrices ${\bf E}_{1}=\begin{bmatrix}1&0\\ 0&1\end{bmatrix}$, ${\bf E}_{2}=\begin{bmatrix}-\cos(2\theta_{i})&\sin(2\theta_{i})\\ \sin(2\theta_{i})&\cos(2\theta_{i})\end{bmatrix}$, and ${\bf E}_{3}=\begin{bmatrix}\sin(2\theta_{i})&\cos(2\theta_{i})\\ \cos(2\theta_{i})&-\sin(2\theta_{i})\end{bmatrix}$ form an orthogonal basis for all $2\times 2$ symmetric matrices under the inner product $\langle A,B\rangle=\operatorname{Tr}(A^{T}B)$. Using this basis, it is straightforward to find that the $N$ equations of (1.6) are equivalent to the following evolution equations on the $2N$ eigenvalues $\lambda_{i}^{1}$, $\lambda_{i}^{2}$ and $N$ eigenangles $\theta_{i}:$

[TABLE]

[TABLE]

[TABLE]

More generally if one carries out the same calculation for the model with external frequency forcing and constraints, Equation (1.7), one finds

[TABLE]

[TABLE]

[TABLE]

where \Omega_{i}=\left[\begin{array}[]{cc}0&\omega_{i}\\ -\omega_{i}&0\end{array}\right] and $\mu_{i}=\sum_{j}\frac{\operatorname{Tr}({\bf M}_{i}[{\bf M}_{j},[{\bf M}_{i},{\bf M}_{j}]])}{\operatorname{Tr}({\bf M}_{i}^{2})}=-\frac{1}{2}\sum_{j}\frac{(\lambda_{j}^{1}-\lambda_{j}^{2})^{2}(\lambda_{i}^{1}-\lambda_{i}^{2})^{2}\sin^{2}2(\theta_{i}-\theta_{j})}{(\lambda_{i}^{1})^{2}+(\lambda_{i}^{2})^{2}}$. Note that $\mu_{i}$ is always negative so it tends to increase the size of the eigenvalues.

Together, Equations (1.11) and (1.12) takes a form similar to that of classical Kuramoto, but with a dynamic coupling parameter that evolves in response to the phase-differences. Such models are very natural in many applications, including neuroscience, and have received considerable attention in the literature [1, 10, 13]. In the theoretical neuroscience literature, such systems are called Hebbian if the coupling strength increases if the oscillators are in phase, and decreases if the oscillators are out of phase, a behavior suggested in the classical neuroscience work of Hebb[12]. The behavior here is similar: in the absence of the constraint the coupling strength $(\lambda_{j}^{2}-\lambda_{j}^{1})^{2}$ is always decreasing, but the rate of decrease goes to zero as the alignment difference $\theta_{i}-\theta_{j}$ approaches zero. Enforcing the constraint (in other words including the $\mu_{i}$ terms) gives a more classically Hebbian behavior – the coupling strength is increasing in time for oscillators that are in phase and decreasing in time for oscillators that are out of phase.

That said, one obvious difference between (1.12) and the classical Kuramoto model is a rescaling of the angle difference by 4: the righthand side of (1.12) is a periodic function of $\theta_{i}$ with period $\pi/2$ rather than period $2\pi,$ as in Kuramoto. However, there is a clear geometric interpretation of this factor of four, since we are considering alignment of an orthogonal frame rather than individual vectors. From the point of view of commutation it is not important which eigenvectors of ${\bf M}_{i}$ align with which eigenvectors of ${\bf M}_{j}$, so one expects a invariance under rotations through angle $\pi/2$, whence the factor of four. Another interesting observation is that, since the coupling between frames is a function of the eigenvalue difference $(\lambda_{j}^{2}-\lambda_{j}^{1})$ there are two ways in which the system can move toward a fixed point: either the frames can align $\theta_{i}=\theta_{j}$ or the eigenvalues can degenerate $\lambda_{j}^{2}=\lambda_{j}^{1}$ in which case ${\bf M}_{j}$ decouples from the rest of the matrices, as it is a multiple of the identity and thus commutes with all the matrices. One fixed point of this system is the case in which all the angles, $\theta_{i}$, are equal. This leads to a set of matrices which all commute for any $\Delta\lambda_{i}$. Hence, the right hand side of (1.6) is equal to zero and we are at a fixed point. We analyze the stability of these fixed points in Section 3.

Figure 1 demonstrates the dynamics of $\eqref{eq:MatrixMoto}$ through a numerical example for the case of three $2\times 2$ matrices. In figure 1 we show a plot of the Hilbert-Schmidt norm of the pairwise commutators of each of the matrices as a function of time together with a plot of the angles between the eigenvector frames as a function of time. It is clear that the matrices evolve to their minimum energy, which is the case in which all of the matrices commute pairwise.

The next figure, (2), depicts the evolution of the eigenframes for eight $2\times 2$ matrices. Each matrix is represented by a pair of lines. The lines point in the direction of the eigenvector, and the length is proportional to the corresponding eigenvalue, so it is possible to determine the eigenvectors and the magnitude (though not the sign) of the eigenvalues from such a plot. The eigenframes are depicted at times $t=0,.001,.002,.003,.004,.005$. Note that the evolution occurs relatively quickly in this case, as we have not scaled the interaction with $\frac{1}{N}$ as is customary with the Kuramoto model. One can see that the lengths of the eigenvalues changes as the matrices evolve, and that there is a strong tendency for the eigenvalues of a given matrix to get closer together. (This is apparent from the equations, as the trace of each matrix is conserved, while the eigenvalue difference is a decreasing function.) We also note that the matrices with the largest eigenvalues tend to dominate the dynamics – they “pull” the small matrices towards them. This is evident from the equations of motion, where the coupling strength between matrices $i$ and $j$ is proportional to the square of the difference of the eigenvalues of ${\bf M}_{j}$.

2 Twist States

It is well-known that, in addition to synchronous or phase-locked states the Kuramoto and related models admit twist or splay states, those where $\theta_{j}=\frac{2\pi(j-1)n}{N}$, with the integer $n$ representing the winding number of the solution. Mirollo [18] gave a very general topological existence for proof of such states for models of the form

[TABLE]

where the right shift operator $\sigma:\theta_{j}\mapsto\theta_{j+1\mod N}$ acts equivariantly on the time evolution operator $F_{t}$: $F_{t}\circ\sigma=\sigma F_{t}$. In a very different vein Ferguson [7] and Delbays, Coletta and Jacquod [4] constructed twist-type solutions for the Kuramoto model on very general graphs (which typically do not satisfy Mirollo’s equivariance criteria.) Medvedev and Tang studied the stability of such states on Cayley and random graphs [17].

The properties of these depend in an interesting way on the topology of the graph. In the case of the complete graph topology (so-called all-to-all coupling) the twist states are always unstable, whereas if the graph is a single cycle the twist states are stable as long as the winding number $n$ is small enough.

In this section we construct the analogs of the twist states for the Kuramoto flow described by Equation (1.6). Interestingly the natural analog of the twist states is not a fixed point of the motion, but rather is a dynamic state on a submanifold along which all of the matrices exhibit slow (algebraic) decay to a multiple of the identity matrix. We analyze the stability of this slow decay: by a dynamic rescaling we can express these solutions as fixed points of related equations, in essence moving along an orbit of the scaling group with the solution. In this co-scaling coordinate system we analyze the stability of the fixed point and show that it is unstable.

Definition 2.1**.**

We define a twist state in $\mathbb{R}^{2}$ as one for which the matrices $\{{\bf M}_{1},{\bf M}_{2},...,{\bf M}_{N}\}$ have eigenvectors pointing in directions which are evenly spaced around the unit circle, and the eigenvalue difference (the effective coupling strength) is constant: in other words $\Delta\lambda_{j}:=\lambda_{j}^{2}-\lambda_{j}^{1}$ is independent of $j$ and $\theta_{j}=\frac{\pi(j-1)n}{2N}$.

Again we note the factor of $4$ in the definition of the twist states for the matrix model versus the classical Kuramoto model: as before this reflects the alignment of orthogonal frames in ${\mathbb{R}}^{2}$ as opposed to angles. We see that the righthand sides of equations $\eqref{eq:2x2EigenvalueSumDynamics}\text{ and }\eqref{eq:2x2Eigen-angleDynamics}$ are equal to zero for this choice of $\Delta\lambda$ and $\boldsymbol{\theta}$, implying that for this initial data the eigenvectors of ${\bf M}_{j}$ and the traces of ${\bf M}_{j}$ do not evolve in time, and thus we get a single equation $\eqref{eq:2x2CouplingDynamics}$ governing the evolution of the eigenvalue difference $\lambda_{j}^{2}-\lambda_{j}^{1}$. It is clear that the right-hand side of equation $\eqref{eq:2x2CouplingDynamics}$ is negative, so the eigenvalue difference decays in time. The decay is clearly algebraic, and given the fact that the righthand side of equation $\eqref{eq:2x2CouplingDynamics}$ is cubic in $\Delta\lambda$ it is simple to compute that the eigenvalue difference decays like $t^{-\frac{1}{2}}$. Thus we have an invariant submanifold where $\theta_{j}=\frac{\pi(j-1)n}{2N}$ and $\Delta\lambda_{j}$ is independent of $j$, and along this submanifold all of the matrices decay in time to multiples of the identity matrix: ${\bf M}_{j}\rightarrow\frac{\operatorname{Tr}{\bf M}_{j}}{2}I.$

Our goal is to understand the stability of this simple decay along a submanifold to a general perturbation. To begin with we again note that the dynamics of the Equation (1.6) is invariant under the $N$ parameter family of maps ${\bf M}_{j}\mapsto{\bf M}_{j}+\alpha_{j}I$, where $I$ is the identity matrix. Because of this invariance we can always assume that $\operatorname{Tr}{\bf M}_{j}=0.$ Our next step is to do a dynamic rescaling to turn these decaying solutions into fixed points of a modified equation. More specifically we let $(\lambda_{j}^{2}-\lambda_{j}^{1})=\alpha(t)$. Then, $\eqref{eq:2x2CouplingDynamics}$ becomes

[TABLE]

An easy calculation shows that the righthand side sums to $-\frac{N\alpha^{3}}{2},$ leading to the separable equation $\frac{d\alpha}{dt}=-\alpha^{3}\frac{N}{2}.$ Solving for $\alpha$, we get

[TABLE]

and thus we see that the eigenvalues decay algebraically. In order to study the stability of this invariant submanifold we re-scale (1.11) in order to obtain fixed points. To do this, we let ${\bf M}_{j}(t)=\alpha(t)\tilde{{\bf M}}_{j}(t)$ Then,

[TABLE]

[TABLE]

Using the fact that $\alpha$ satisfies $\frac{d\alpha}{dt}=-\alpha^{3}\frac{N}{2}$, we get

[TABLE]

In order to remove the algebraic decay from this system, we rescale time via $\alpha^{2}(t)dt=ds$. That is in the rescaled time

[TABLE]

a twist solution is a fixed point of

[TABLE]

Note that this is exactly the form of the contrained gradient flow, if we constrain the Hilbert-Schmidt norm of each matrix ${\tilde{\bf M}}_{i}$ to remain constant.

As we did in Section 1.2, using the spectral decomposition of $\tilde{\bf{\bf M}}_{j}={\bf R}(\theta_{j}){\bf\Lambda}_{j}{\bf R}(-\theta_{j})$ yields the equations of motion for the rescaled system in spectral coordinates

[TABLE]

To analyze the stability of the twist states we compute the Jacobian. Recall that, due to the invariance under the map ${\bf M}_{i}\mapsto{\bf M}_{i}+\alpha_{i}{\bf I}$, the dynamics is trivial in the $\operatorname{Tr}({\bf M}_{i})$ directions. We can see from Equation (2.2) that the Jacobian has at least $N$ elements in the kernel. What is more interesting is the remaining $2N$ degrees of freedom represented by $\theta_{i},\Delta\lambda_{i}:=\lambda_{i}^{2}-\lambda_{i}^{1}$, which evolve non-trivially. Thus we need to analyze the $2N\times 2N$ matrix that results from the Equations (2.3) and (2.4). Defining $\Delta\lambda_{j}:=(\lambda_{j}^{2}-\lambda_{j}^{1})$, the equations of motion become

[TABLE]

where

[TABLE]

and

[TABLE]

Thus we find that the Jacobian ${\bf J}$ can be written in the block form

[TABLE]

where the partial derivatives of each non-zero block are of the form

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

To compute the eigenvalues of ${\bf J}$, we take advantage of the structure of the individual blocks. Using the evenness of the cosine function, the oddness of the sine and the equivalent forms of polar angles, we have that

• •

$\frac{\partial g}{\partial\theta}$ and $\frac{\partial f}{\partial\Delta\lambda}$ are real symmetric matrices,

• •

$\frac{\partial g}{\partial\Delta\lambda}$ and $\frac{\partial f}{\partial\theta}$ are antisymmetric,

• •

$\frac{\partial g}{\partial\Delta\lambda}$ and $\frac{\partial f}{\partial\theta}$ are algebraically related by $\frac{\partial f}{\partial\theta}=4\frac{\partial g}{\partial\Delta\lambda}^{T}$, and

• •

Each block is individually a circulant matrix.

Here we remind the reader of a standard fact about the spectrum of circulant matrices: see, for example, the text of Terras[20] for details.

Lemma 2.2**.**

Suppose that ${\bf A}$ is an $n\times n$ circulant matrix: in other words that ${\bf A}_{ij}=c_{(i-j)\mod n}$ for some vector ${\bf c}\in{\mathbb{R}}^{n}$. Then the eigenvectors are given by ${\bf v}^{(\alpha)}$ with ${\bf v}^{(\alpha)}_{j}=\exp(\frac{2\pi i(j-1)\alpha}{n})$ and the eigenvalues are given by $\lambda^{\alpha}=\sum_{j=1}^{n}c_{j}\exp(\frac{2\pi i(j-1)\alpha}{n})$. In other words ${\bf A}$ is diagonalized by the discrete Fourier transform, and the eigenvalues are the discrete Fourier transform of the vector ${\bf c}$.

The fact that the Jacobian is of block-circulant form allows us to compute the eigenvalues of the Jacobian, leading to the following stability theorem.

Theorem 2.3**.**

Any collection $\{{\bf M}_{1},{\bf M}_{2},...,{\bf M}_{N}\}$ of $2\times 2$ matrices in a twist state (Definition 2.1) is an unstable orbit of Equation (1.6). More precisely, after re-scaling (2.1) for the algebraic decay in the eigenvalue difference $\Delta\lambda_{j}$, the Jacobian (2.5) evaluated at a twist state has the following eigenvalues

[TABLE]

Proof.

Given the block-circulant nature of the matrix the theorem follows from fairly straightforward computations. First note that since the matrix is block circulant the individual blocks commute (all circulant matrices have a common eigenbasis) and the eigenvectors may be assumed to be a direct sum of a vector in ${\mathbb{R}}^{2}$ with the eigenvector of a circulant matrix

[TABLE]

where, as previously, ${\bf v}^{(\alpha)}_{j}=\exp\left(\frac{2\pi i(j-1)\alpha}{n}\right)$. Let $a_{\alpha}$ denote the eigenvalue for block $\frac{\partial g}{\partial\theta}$ associated with the eigenvector ${\bf v}^{(\alpha)}$, $b_{\alpha}$ for block $\frac{\partial g}{\partial\Delta\lambda}$, and $c_{\alpha}$ for block $\frac{\partial f}{\partial\Delta\lambda}$. Note that the eigenvalues of the fourth block are $4\bar{b}_{\alpha}$, by the algebraic relationship of the off diagonal blocks. Then, focusing solely on the non-zero block of ${\bf J}$, we have

[TABLE]

or equivalently

[TABLE]

We see that for each $\alpha$, our larger eigenvalue problem reduces to the associated $2\times 2$ eigenvalue problem for the matrix defined by each block’s individual eigenvalue. Doing this for all $\alpha$ will yield $2N$ eigenvalues for the Jacobian matrix ${\bf J}$. A straightforward though somewhat lengthy computation – see Appendix B for details – gives the following expressions for $a_{\alpha},b_{\alpha},c_{\alpha}$

[TABLE]

To finish our analysis, we return to Equation (2.8) and the coefficient matrix $\begin{bmatrix}a_{\alpha}&b_{\alpha}\\ 4\bar{b}_{\alpha}&c_{\alpha}\end{bmatrix}$. For each $\alpha$ the above matrix gives two eigenvalues; the spectrum of the linearization being the union over all $\alpha\in\{0,1,\ldots,N-1\}.$ The only values of $\alpha$ that are important are $\alpha=0,1,N-1$, since the remaining $\alpha$ values give the zero matrix. When $\alpha=0$ the coefficient matrix is $\begin{bmatrix}0&0\\ 0&-N\end{bmatrix}$ with eigenvalues [math] and $-N$ and associated eigenvectors $[{\bf 0},{\bf v}^{(0)},{\bf 0}]^{T}$ and $[{\bf 0},{\bf 0},{\bf v}^{(0)}]^{T}$ of ${\bf J}$, respectively. Note that ${\bf v}^{(0)}$ is the vector $[1,1,1,\ldots,1]^{T}$. When $\alpha=1$ the coefficient matrix is $\begin{bmatrix}N/2&-Ni/4\\ Ni&N/2\end{bmatrix}$ with eigenvalues [math] and $N$. The associated eigenvectors of ${\bf J}$ are $[{\bf 0},{\bf v}^{(1)},-2i{\bf v}^{(1)}]^{T}$ and $[{\bf 0},(-Ni/4){\bf v}^{(1)},(3N/2){\bf v}^{(1)}]^{T}$, respectively. Finally for $\alpha=N-1$, we get $\begin{bmatrix}N/2&Ni/4\\ -Ni&N/2\end{bmatrix}$ which again has eigenvalues [math] and $N$. For all other $\alpha$ the coefficient matrix is the zero matrix ${\bf 0}$. This yields two 0 eigenvalues for each $\alpha\notin\{0,1,N-1\}$. Recall that we earlier saw that there are an additional $N$ zero eigenvalues corresponding to the coordinates $\lambda_{i}^{1}+\lambda_{i}^{2}=\operatorname{Tr}({\bf M}_{i}).$ Thus the $3N$ eigenvalues of ${\bf J}$ are given by

[TABLE]

∎

Figures 3 – 6 demonstrate the dynamics of a twist state and the instability of a twist state (2.1) to perturbation. We look at three collections of matrices; the first with the exact initial conditions for a twist state, the second being an example where the eigenvalue difference $\Delta\lambda$ is not constant throughout, and the third where the rotation angles are not quite equally spaced.

For the collection of matrices in Figure 3, the eigenvalue frames are all $\pi/6$ apart and the eigenvalue difference is $\Delta\lambda=2$. We numerically integrate the matrix Kuramoto flow with initial conditions satisfying the conditions of Definition 2.1 to study the dynamics. Figure 3 demonstrates there is an exact solution where the eigenframe does not change but where the difference in the eigenvalues decays algebraically. The left-most plot depicts the three eigenframes and associated eigenvalues. For each matrix we plot the eigenvectors in the $x-y$ plane, scaled so that the length is the same as the magnitude of the associated eigenvalue. One of the matrices has eigenvalues $\lambda_{1}=5$ and $\lambda_{2}=3$ with associated eigenvectors $(1,0)$ and $(0,1)$. This is represented in the plot by a line segment of length $5$ in the $(1,0)$ direction and a segment of length $3$ in the $(0,1)$ direction. The other two matrices have eigenvalues $\lambda_{1}=4,\lambda_{2}=2$ and $\lambda_{1}=1,\lambda_{2}=3$ with associated eigenvectors $v_{1}=[\cos\frac{\pi}{6},-\sin\frac{\pi}{6}]^{t}$ and $v_{2}=[\sin\frac{\pi}{6},\cos\frac{\pi}{6}]^{t}$ and $v_{1}=[\cos\frac{\pi}{3},-\sin\frac{\pi}{3}]^{t}$ and $v_{2}=[\sin\frac{\pi}{3},\cos\frac{\pi}{3}]^{t}$ respectively. The next two frames show the solutions at times $t=1.0$ and $t=8.0$ respectively. It is apparent that the eigenvectors have not changed, but that the eigenvalue difference for each matrix has converged to zero. In the third frame the axes are of approximately equal length indicating that each matrix is close to a multiple of the identity. This clearly shows the dynamics on the twist manifold.

Figure 4 demonstrates an example where the eigenvalue difference is not the same for all three matrices. Again the rotational difference is $\pi/6$, but two of the matrices have a $\Delta\lambda=2$ and one has $\Delta\lambda=3$. The initial condition here is similar looking, but the dynamics is very different. With these perturbed initial conditions the eigenframes rotate and quickly align, dominated by the eigenframe with the largest eigenvalues. Note that the time-scale is much shorter here than in the previous example: the three plots depict times $t=0,0.4,1.0$, as compared with $t=0,1.0,8.0$ in the previous figure. This is because the instability is exponential, while the decay of the non-identity portions of the matrices on the twist manifold is algebraic. Figure 5 shows a similar experiment, where the variation of $\Delta\lambda_{i}$ is much smaller: the first and third matrices have $\Delta\lambda=2$ while the second matrix has $\Delta\lambda=2.01$. We see a similar phenomenon as in the prior figure, though on a somewhat longer time scale. The eigenframes rotate in order to align.

Figure 6 shows one final experiment to illustrate the instability of the twist state. In this example the eigenvalue differences are again $2$ for each matrix, as required for a twist state, but we have perturbed one of the angles slightly. The first and third eigenframes are at angle [math] and $\frac{\pi}{3}$ but the second eigenframe is at angle $11\pi/60$ off the axis. Again we see that the eigenframes do not remain fixed, but rather they rapidly rotate to coincide.

3 Stability of Commuting Matrices Fixed Points

It is clear that any set of $N$ pairwise commuting $n\times n$ matrices $\{{\bf M}_{i}\}_{i=1}^{N}$ gives a fixed point of (1.6) or (more correctly) a point on a fixed manifold. Our goal in this section is to show that this point is conditionally stable - the linearization of the flow about this point is negative semi-definite, with the dimension of the null-space equal to the dimension of the manifold of fixed points. To begin with we first count the dimension of the fixed manifold. In addition to setting our expectations with respect to the dimension of the kernel of the linearized operator it will help to motivate the construction of the basis that diagonalizes the linearization. First of all we claim that the manifold of commuting fixed points has dimension $nN+\frac{n(n-1)}{2}.$ To see this note that a family of symmetric matrices is commuting if and only if they can be simultaneously diagonalized, so each family is parameterized by $n$ eigenvalues for each of the $N$ matrices together with a common eigenbasis given by a matrix in $SO(n).$ Thus we expect that the linearized operator will have a null-space of dimension $nN+\frac{n(n-1)}{2}$. We are interested in non-zero eigenvalues of the linearization, as these determine the stability of the fixed point. We begin by presenting a small lemma, which gives a useful representation for the tangent space to the space of symmetric matrices. Of course the vector space of symmetric $n\times n$ matrices is isomorphic to ${\mathbb{R}}^{\frac{n(n+1)}{2}}$, and the tangent space to ${\mathbb{R}}^{k}$ is ${\mathbb{R}}^{k}$, but it is convenient to give an explicit decomposition into commuting and non-commuting directions.

Lemma 3.1**.**

Given a real symmetric $n\times n$ matrix ${\bf M}$ with distinct eigenvalues $\lambda_{i}$ and orthonormal eigenvectors ${\bf v}_{i}$. There is an orthonormal decomposition of the space of symmetric $n\times n$ matrices into commuting and non-commuting components as follows: there is an $n$ dimensional subspace of matrices which commute with ${\bf M}$ and an $\frac{n(n-1)}{2}$ dimensional subspace of matrices which do not commute with ${\bf M}$ spanned respectively by $\{{\bf C}_{\alpha}\}_{\alpha=1}^{n}$ and $\{{\bf N}_{\alpha\beta}\}_{\alpha,\beta~{}\alpha>\beta}^{n}$

[TABLE]

This basis is orthonormal under the usual inner product on symmetric matrices $\langle A,B\rangle=\operatorname{Tr}(A^{T}B).$

Proof.

It is easy to see the ${\bf M}{\bf C}_{\alpha}={\bf M}{\bf v}_{\alpha}\otimes{\bf v}_{\alpha}=\lambda_{\alpha}{\bf v}_{\alpha}\otimes{\bf v}_{\alpha}={\bf C}_{\alpha}{\bf M}$ so $[{\bf C}_{\alpha},{\bf M}]=0$ so the matrices commute. Similarly we have that

[TABLE]

which is clearly non-zero if $\lambda_{\alpha}$ and $\lambda_{\beta}$ are distinct. (In the non-generic case where ${\bf M}$ has repeated eigenvalues there are “extra” directions that commute with ${\bf M}$.) As the eigenvectors $\{{\bf v}_{i}\}$ can be chosen to be orthonormal the basis vectors ${\bf C}_{\alpha},{\bf N}_{\alpha\beta}$ inherit this property. ∎

Next we consider the linearized flow. We assume that we have a family of $N$ commuting $n\times n$ matrices $\{{\bf M}^{0}_{i}\}_{i=1}^{N}$ with $[{\bf M}_{i}^{0},{\bf M}_{j}^{0}]=0$ for any fixed $i,j$, which forms a fixed point of the matrix Kuramoto flow. We will let $\tilde{{\bf M}}_{i}$ represent the perturbation of ${\bf M}_{i}$: in other words ${\bf M}_{i}={\bf M}_{i}^{0}+\tilde{\bf M}_{i}.$ With this notation it is simple to compute that $\tilde{{\bf M}}_{i}$ satisfies the following linearized flow

[TABLE]

which we can write more compactly as

[TABLE]

where ${\mathcal{L}}$ is a linear operator on the space of $N$ symmetric $n\times n$ matrices (isomorphic to ${\mathbb{R}}^{N\frac{n(n+1)}{2}})$. Interestingly we can compute the spectrum of this linear operator quite explicitly – this is the content of the next proposition.

Proposition 3.2**.**

For a given set of commuting matrices $\{{\bf M}_{i}^{0}\}_{i=1}^{N}$ the spectrum of $\boldsymbol{\mathcal{L}}$ as an operator from ${\mathbb{R}}^{N\frac{n(n+1)}{2}}$ to ${\mathbb{R}}^{N\frac{n(n+1)}{2}}$ is given as follows. Suppose that ${\bf M}_{i}^{0}$ has eigenvalues $\lambda_{i}^{\alpha}$ for $\alpha\in\{1,2,\ldots,n\}$. The eigenvalues of ${\mathcal{L}}$ are given by

[TABLE]

So generically $\boldsymbol{\mathcal{L}}$ has a kernel of dimension $nN+\frac{n(n-1)}{2}$ and $\frac{n(n-1)}{2}$ negative eigenvalues, each of multiplicity $N-1$. Of course for non-generic situations we could get additional degeneracies or additional elements of the kernel.

Proof.

Commuting matrices share a common set of eigenvectors, a fact we use to construct a particularly convenient orthonormal basis for the set of $n\times n$ symmetric matrices. Given a set $\{{\bf v}_{\alpha}\}_{\alpha=1}^{n}$ of orthonormal eigenvectors of the commuting family of matrices $\{{\bf M}_{i}^{0}\}_{i=1}^{N}$ we define a basis for for $n\times n$ symmetric matrices as follows.

[TABLE]

It is straightforward to check that (assuming that $\{{\bf v}_{\alpha}\}_{\alpha=1}^{n}$ are orthogonal) this basis is orthonormal under the usual matrix inner product $\langle A,B\rangle=\operatorname{Tr}(A^{T}B)$. Using this basis, $\tilde{{\bf M}}_{i}$ can be written as

[TABLE]

and we have the following commutation rules:

[TABLE]

We can see that the vectors $\{{\bf v}_{\alpha}\otimes{\bf v}_{\alpha}\}_{\alpha=1}^{n}$ lie in the kernel of the operator $[{\bf M}_{1}^{0},\cdot]$: they represent the “commuting” directions, and will be responsible for the $nN$ dimensional part of the kernel. The remaining directions $\{{\bf v}_{\alpha}\otimes{\bf v}_{\beta}+{\bf v}_{\beta}\otimes{\bf v}_{\alpha}\}_{\alpha>\beta}^{n}$ behave non-trivially under commutation with ${\bf M}_{i}^{0}$. Each of these will give a block, on which the linearized operator looks like a graph Laplacian. There will be one zero eigenvalue in each of these blocks, giving the additional $\binom{n}{2}$ directions in the kernel.

Given the above commutation relations we next compute the linearization using this representation of $\tilde{{\bf M}}_{i}$ for a fixed $i$, $j$.

[TABLE]

In a similar fashion we find that

[TABLE]

and finally

[TABLE]

Thus, (3.2) can be written as

[TABLE]

where the coefficients $a_{\alpha}$ and $b_{\alpha\beta}$ evolve according to

[TABLE]

Thus we have block-diagonalized the linearized equations of motion. The evolution in $a_{\alpha}$, the commuting direction, is trivial. We can represent the $\frac{d{\bf b}_{\alpha_{\beta}}}{dt}$ by

[TABLE]

where ${\bf b}_{\alpha\beta}=[b_{\alpha\beta}^{1},b_{\alpha\beta}^{2},b_{\alpha\beta}^{3},\ldots b_{\alpha\beta}^{N}]^{T}$ and the matrix ${\bf L}_{\alpha\beta}={\bf v}\otimes{\bf v}-\mu{\bf I}$ where

[TABLE]

and $\mu=\sum\limits_{j=1}^{N}(\lambda_{j}^{\alpha}-\lambda_{j}^{\beta})^{2}$. It is easy to see that the eigenvalues of a matrix ${\bf v}\otimes{\bf v}-\mu{\bf I}$ are $\|{\bf v}\|^{2}-\mu$, with multiplicity one (the corresponding eigenvector if ${\bf v}$) and $-\mu$ with multiplicity $N-1$, the eigenvectors being any orthogonal basis for the $(N-1)$ dimensional subspace orthogonal to ${\bf v}$.

Thus the eigenvalues and eigenvectors of ${\bf L}$ are:

[TABLE]

with corresponding eigenvectors

[TABLE]

∎

3.1 Numerical Simulations

We conclude this section with some numerical experiments. It appears that for generic initial conditions in the absence of external forcing (${\bf\Omega}_{i}=0$) the solutions always converge to a family of commuting matrices. We illustrate this with a couple of numerical simulations. The first numerical simulation depicts the evolution of a set of three $3\times 3$ matrices that are initially pair-wise non-commuting. Since the eigenframes are difficult to visualize in three and higher dimensions we instead choose to plot the Hilbert-Schmidt norms of the pairwise commutators, $\|[{\bf M}_{i},{\bf M}_{j}]\|^{2}$ as functions of time. The results of one such expermient are depicted in Figure 7. We see that there is a very rapid decay of the Hilbert-Schmidt norms to zero, indicating that the matrices are asymptotically commuting. Figure 8 depicts a similar experiment with three $4\times 4$ matrices. Again we see a similar behavior – there is a rapid decy in the Hilbert-Schmidt norms as the eigenframes move to align. We have repeated these experiments a number of times, varying both the number of matrices and the size of the matrices. In all cases we have observed very rpid convergence to a commuting set of matrices.

4 Conclusion

In this paper we proposed a new, purely real-valued, non-Abelian generalization of the Kuramoto model to symmetric matrix-valued variables, and proved three results for this model. We have shown that in the case of $2\times 2$ matrices the system can be described by three equations. One of these equations is similar to the classical Kuramoto model with dynamic coupling. However, the angle difference is multiplied by 4 as a result of our consideration of the aligning of eigenvector frames as opposed to angles. We have also defined twist states for the matrix-valued model in $\mathbb{R}^{2}$ to be the case in which the matrices have a constant eigenvalue difference and have eigenvectors pointing in directions which are equally spaced around the unit circle. We have proven the twist states are dynamically unstable. Additionally, we have proven that fixed points associated with families of commuting matrices of any size are conditionally stable.

Thus far we have only considered twist states in the case of $2\times 2$ matrices. In future research, it would be interesting to investigate twist states for larger matrices. In addition, throughout this paper we have assumed all oscillators to have the same natural frequency and the absence of a constrained gradient flow. In other words, analysis of Equation (1.7) should yield some interesting results.

5 Acknowledgments

J.C.B. and S.E.S. would like to acknowledge support under NSF grant NSF- DMS 1615418.

T.E.C. would like to acknowledge support from Caterpillar Fellowship Grant at Bradley University.

Appendix A Derivation of the Gradient Flow

Recall the energy function defined in Section 1,

[TABLE]

The energy $E$ is a scalar valued function of $N$, $k\times k$ matrices ${\bf M}_{i}$ and $\frac{\partial E}{\partial{\bf M}_{i}}$ is a $k\times k$ matrix whose entries are the derivatives of $E$ with respect to the corresponding entry of ${\bf M}_{i}$. Individual components of a matrix will be denoted $M_{i;\alpha\beta}$ with Greek indices for the matrix components. As in the main text the derivatives are defined

[TABLE]

where $M_{i;\alpha\beta}$ is the $(\alpha,\beta)$ entry of the matrix ${\bf M}_{i}$.

We consider $\frac{1}{4}\operatorname{Tr}[{\bf M}_{i},{\bf M}_{j}]^{2}$, a single term in the sum of the energy function and want to compute the derivative with respect to $M_{i;\alpha\beta}$. Using the fact that the trace operator is linear and invariant under cyclic permutations of the matrix product, we have that

[TABLE]

Applying the identity that $\frac{\partial M_{j;\tau\kappa}}{\partial M_{i;\alpha\beta}}=\delta_{ij}\left(\delta_{\alpha\tau}\delta_{\beta\kappa}+\delta_{\alpha\kappa}\delta_{\beta\tau}\right)$

[TABLE]

Note that the resultant here is the $\alpha\beta$ entry of the matrix product

[TABLE]

Then, for fixed $j$,

[TABLE]

Then, summing over all the terms yields the matrix Kuramoto flow (1.6)

[TABLE]

Appendix B Computation of the eigenvalues and eigenvectors of the linearization about a twist state.

In order to apply Equation (2.8), for each $\alpha\in\{0,1,2,...,N-1\}$ we need to determine the associated eigenvalue for each block of ${\bf J}$, Equation (2.5). Since the first entry of ${\bf v}^{(\alpha)}$ is $1$ the associated eigenvalue is given by the dot product of the first row of each block with the vector ${\bf v}^{(\alpha)}$. The following properties of the discrete Fourier transform will be useful here.

[TABLE]

We begin by looking at the $\frac{\partial g}{\partial\theta}$ block. The top row of this block is

[TABLE]

To find the $k^{th}$ eigenvalue $a_{k}$ we take the dot product of this row by ${\bf v}^{(k)}$ yields

[TABLE]

By the given identities this sum is equal to $N/2$ if $k=1$ or $k=N-1$ and is zero otherwise.

The eigenvalues for the other two blocks follow in the same way. The first row of the $\frac{\partial g}{\partial\Delta\lambda}$ block is given by

[TABLE]

It is clear that since the entries of the row are given by $\sin\frac{2\pi j}{N}$ the only non-zero Fourier coefficients will be $k=1$ and $k=N-1$ giving $b_{1}=-iN/4$ and $b_{N-1}=iN/4$, with all of the remaining eigenvalues equal to zero.

Finally we consider the block $\frac{\partial f}{\partial\Delta\theta}$. The first row of this block looks like

[TABLE]

Since this term is the sum of a constant and a cosine term we will have three non-zero Fourier coefficients and thus three non-zero eigenvalues. Applying the above identities gives

[TABLE]

with all of the remaining eigenvalues zero.