$\displaystyle\left.\frac{\partial\psi}{\partial\mathbf{x}}\right|_{\mathbf{x}}=\left.\frac{\partial\psi}{\partial\mathbf{x}}\right|_{\mathbf{x}^{\gamma}}+\left.\frac{\partial^{2}\psi}{\partial\mathbf{x}^{2}}\right|_{\mathbf{x}^{\gamma}}\Delta\mathbf{x}+\cdots$
(99)

and hence expanding $\Delta\mathbf{x}$ as in (83) and comparing powers of
$\psi$ with (11c) we see that
$I^{(0)}=\left.\frac{\partial\psi}{\partial\mathbf{x}}\right|_{\mathbf{x}^{\gamma}}$
and
$I^{(1)}=\left.\frac{\partial^{2}\psi}{\partial\mathbf{x}^{2}}\right|_{\mathbf{x}^{\gamma}}g^{(1)}$.
Then $\frac{{\rm d}I^{(0)}}{{\rm
d}t}=\left.\frac{\partial^{2}\psi}{\partial\mathbf{x}^{2}}\right|_{\mathbf{x}^{\gamma}(\theta)}^{T}\mathbf{F}(\mathbf{x}^{\gamma}(\theta))$.
Therefore (98) becomes

$\displaystyle\kappa$
$\displaystyle=\left(\left.\frac{\partial^{2}\psi}{\partial\mathbf{x}^{2}}\right|_{\mathbf{x}^{\gamma}}g^{(1)}\right)\cdot\mathbf{F}(\mathbf{x}^{\gamma}(\theta))+I^{(0)}(\theta)\cdot
Jg^{(1)}(\theta)$
$\displaystyle=g^{(1)}(\theta)\cdot\left.\frac{\partial^{2}\psi}{\partial\mathbf{x}^{2}}\right|_{\mathbf{x}^{\gamma}}^{T}\mathbf{F}(\mathbf{x}^{\gamma}(\theta))+g^{(1)}(\theta)\cdot
J^{T}I^{(0)}(\theta)$ $\displaystyle=g^{(1)}(\theta)\cdot\left(\frac{{\rm
d}I^{(0)}}{{\rm d}t}+J^{T}I^{(0)}(\theta)\right)=\kappa g^{(1)}(\theta)\cdot
I^{(0)}(\theta),$

where (12c) is used in the final step. Thus the required normalisation is
satisfied if $I^{(0)}(0)\cdot g^{(1)}(0)=1$.

##  Appendix B Stability of phase-locked states

###  B.1 General phase-locked states

For a phase locked state $\Phi=(\phi_{1},\ldots,\phi_{N})$ the linearisation
of (20) has $2N\times 2N$ block matrix Jacobian

$\displaystyle\mathcal{J}=\left[\begin{array}[]{@{}c|c@{}}\mathcal{H}^{(1)}(\Phi)&\mathcal{H}^{(2)}(\Phi)\\\
\hline\cr\mathcal{H}^{(3)}(\Phi)&\mathcal{H}^{(4)}(\Phi)\end{array}\right]$

where

$\displaystyle\begin{split}\mathcal{H}_{ij}^{(1)}(\Phi)=&\ \epsilon
w_{ij}\Xi^{(\theta)}_{ij}-\epsilon\delta_{ij}\sum_{k=1}^{N}w_{ik}\Xi^{(\theta)}_{ik},\end{split}$
(100a) $\displaystyle\mathcal{H}_{ij}^{(2)}(\Phi)=$ $\displaystyle\ \epsilon
w_{ij}H_{3}(\phi_{j}-\phi_{i})+\epsilon\delta_{ij}\sum_{k=1}^{N}w_{ik}H_{2}(\phi_{k}-\phi_{i}),$
(100b) $\displaystyle\begin{split}\mathcal{H}_{ij}^{(3)}(\Phi)=&\ \epsilon
w_{ij}\Xi^{(\psi)}_{ij}-\epsilon\delta_{ij}\sum_{k=1}^{N}w_{ik}\Xi^{(\psi)}_{ik},\end{split}$
(100c) $\displaystyle\begin{split}\mathcal{H}_{ij}^{(4)}(\Phi)=&\ \epsilon
w_{ij}H_{6}(\phi_{j}-\phi_{i})+\delta_{ij}\left(\kappa+\varepsilon\sum_{k=1}^{N}w_{ik}H_{5}(\phi_{k}-\phi_{i})\right)\end{split}$
(100d)

for

$\displaystyle\Xi^{(\theta)}_{ij}$
$\displaystyle=H^{\prime}_{1}(\phi_{j}-\phi_{i})+\Psi_{i}H^{\prime}_{2}(\phi_{j}-\phi_{i})+\Psi_{j}H^{\prime}_{3}(\phi_{j}-\phi_{i}),$
$\displaystyle\Xi^{(\psi)}_{ij}$
$\displaystyle=H^{\prime}_{4}(\phi_{j}-\phi_{i})+\Psi_{i}H^{\prime}_{5}(\phi_{j}-\phi_{i})+\Psi_{j}H^{\prime}_{6}(\phi_{j}-\phi_{i}).$

###  B.2 Two-cluster states

For two cluster states where the two clusters have distinct orbits, the
linearisation about $\Phi=(0,\ldots,0,\chi,\ldots,\chi)$ has $2N\times 2N$
block matrix Jacobian (31) where each of the blocks $\mathcal{H}^{(q)}$ are
themselves block matrices of the form (61). The coefficients in the definition
of the blocks (62a)–(62d) can be rearranged into arrays and are found to be
given by

$\displaystyle\begin{split}\mathcal{A}&=\begin{bmatrix}a_{1}&a_{2}\\\
a_{3}&a_{4}\end{bmatrix}=\frac{\epsilon}{N}\begin{bmatrix}H_{1}^{\prime}(0)+{\Psi}_{A}(H_{2}^{\prime}(0)+H_{3}^{\prime}(0))&H_{3}(0)\\\
H_{4}^{\prime}(0)+{\Psi}_{A}(H_{5}^{\prime}(0)+H_{6}^{\prime}(0))&H_{6}(0)\end{bmatrix},\end{split}$
(101a) $\displaystyle\begin{split}\mathcal{B}&=\begin{bmatrix}b_{1}&b_{2}\\\
b_{3}&b_{4}\end{bmatrix}=\frac{\epsilon}{N}\begin{bmatrix}H_{1}^{\prime}(\chi)+{\Psi}_{A}H_{2}^{\prime}(\chi)+{\Psi}_{B}H_{3}^{\prime}(\chi)&H_{3}(\chi)\\\
H_{4}^{\prime}(\chi)+{\Psi}_{A}H_{5}^{\prime}(\chi)+{\Psi}_{B}H_{6}^{\prime}(\chi)&H_{6}(\chi)\end{bmatrix},\end{split}$
(101b) $\displaystyle\begin{split}\mathcal{C}&=\begin{bmatrix}c_{1}&c_{2}\\\
c_{3}&c_{4}\end{bmatrix}=\frac{\epsilon}{N}\begin{bmatrix}H_{1}^{\prime}(-\chi)+{\Psi}_{B}H_{2}^{\prime}(-\chi)+{\Psi}_{A}H_{3}^{\prime}(-\chi)&H_{3}(-\chi)\\\
H_{4}^{\prime}(-\chi)+{\Psi}_{B}H_{5}^{\prime}(-\chi)+{\Psi}_{A}H_{6}^{\prime}(-\chi)&H_{6}(-\chi)\end{bmatrix},\end{split}$
(101c) $\displaystyle\begin{split}\mathcal{D}&=\begin{bmatrix}d_{1}&d_{2}\\\
d_{3}&d_{4}\end{bmatrix}=\frac{\epsilon}{N}\begin{bmatrix}H_{1}^{\prime}(0)+{\Psi}_{B}(H_{2}^{\prime}(0)+H_{3}^{\prime}(0))&H_{3}(0)\\\
H_{4}^{\prime}(0)+{\Psi}_{B}(H_{5}^{\prime}(0)+H_{6}^{\prime}(0))&H_{6}(0)\end{bmatrix},\end{split}$
(101d)
$\displaystyle\begin{split}\mathcal{S^{A}}&=\begin{bmatrix}s^{A}_{1}&s^{A}_{2}\\\
s^{A}_{3}&s^{A}_{4}\end{bmatrix}=\begin{bmatrix}-N_{A}a_{1}-N_{B}b_{1}&\epsilon(N_{A}H_{2}(0)+N_{B}H_{2}(\chi))/N\\\
-N_{A}c_{1}-N_{B}d_{1}&\kappa+\epsilon(N_{A}H_{5}(0)+N_{B}H_{5}(\chi))/N\end{bmatrix},\end{split}$
(101e)
$\displaystyle\begin{split}\mathcal{S^{B}}&=\begin{bmatrix}s^{B}_{1}&s^{B}_{2}\\\
s^{B}_{3}&s^{B}_{4}\end{bmatrix}=\begin{bmatrix}-N_{A}a_{3}-N_{B}b_{3}&\epsilon(N_{A}H_{2}(-\chi)+N_{B}H_{2}(0))/N\\\
-N_{A}c_{3}-N_{B}d_{3}&\kappa+\epsilon(N_{A}H_{5}(-\chi)+N_{B}H_{5}(0))/N\end{bmatrix}.\end{split}$
(101f)

We can then rewrite

$\displaystyle|\mathcal{J}-\mu
I_{2N}|=\begin{vmatrix}\mathcal{A}\otimes\mathbbm{1}_{N_{A}\times
N_{A}}+(\mathcal{S}^{A}-\mu I_{2})\otimes
I_{N_{A}}&\mathcal{B}\otimes\mathbbm{1}_{N_{A}\times N_{B}}\\\
\mathcal{C}\otimes\mathbbm{1}_{N_{B}\times
N_{A}}&\mathcal{D}\otimes\mathbbm{1}_{N_{B}\times N_{B}}+(\mathcal{S}^{B}-\mu
I_{2})\otimes I_{N_{B}}\end{vmatrix}.$ (102)

A generalisation of the matrix determinant lemma Ding and Zhou (2007) says
that for invertible $n\times n$ matrix $P$ and invertible $m\times m$ matrix
$Q$, if $U$ and $V$ are $n\times m$ matrices then

$\displaystyle\det(P+UQV^{T})=\det(Q^{-1}+V^{T}P^{-1}U)\det(Q)\det(P).$ (103)

Taking $P=(\mathcal{S}^{A}-\mu I_{2})\otimes I_{N_{A}}$,
$Q=(\mathcal{S}^{A}-\mu I_{2})^{-1}$ and letting $U=(\mathcal{S}^{A}-\mu
I_{2})\otimes\mathbbm{1}_{N_{A}\times 1}$ and
$V=\mathcal{A}^{T}\otimes\mathbbm{1}_{N_{A}\times 1}$ we observe that since
$UQV^{T}=\mathcal{A}\otimes\mathbbm{1}_{N_{A}\times N_{A}}$ and
$V^{T}P^{-1}U=N_{A}\mathcal{A}$,

$\displaystyle\begin{split}&\det(\mathcal{A}\otimes\mathbbm{1}_{N_{A}\times
N_{A}}+(\mathcal{S}^{A}-\mu I_{2})\otimes I_{N_{A}})\\\ &\ \
=\det(\mathcal{S}^{A}-\mu
I_{2})^{N_{A}-1}\det(\mathcal{S}^{A}+N_{A}\mathcal{A}-\mu I_{2}).\end{split}$
(104)

This reveals that eigenvalues corresponding to intracluster perturbations are
given by the eigenvalues of $\mathcal{S}^{A}$ each with multiplicity $N_{A}-1$
and the eigenvalues of $\mathcal{S}^{B}$ each with multiplicity $N_{B}-1$. The
remaining four eigenvalues can be found by considering intercluster
perturbations
$(\phi_{A},\psi_{A},\phi_{B},\psi_{B})=(\overline{\phi}_{A},{\Psi}_{A},\overline{\phi}_{B},{\Psi}_{B})+(\widetilde{\phi}_{A},\widetilde{\psi}_{A},\widetilde{\phi}_{B},\widetilde{\psi}_{B})$
where $\overline{\phi}_{A}$ and $\overline{\phi}_{B}$ are the phases of the
two clusters in the phase locked state. Then

$\displaystyle\frac{{\rm d}}{{\rm d}t}\begin{bmatrix}\widetilde{\phi}_{A}\\\
\widetilde{\psi}_{A}\\\ \widetilde{\phi}_{B}\\\
\widetilde{\psi}_{B}\end{bmatrix}=\begin{bmatrix}\mathcal{S}^{A}+N_{A}\mathcal{A}&N_{B}\mathcal{B}\\\
N_{A}\mathcal{C}&\mathcal{S}^{B}+N_{B}\mathcal{D}\end{bmatrix}\begin{bmatrix}\widetilde{\phi}_{A}\\\
\widetilde{\psi}_{A}\\\ \widetilde{\phi}_{B}\\\
\widetilde{\psi}_{B}\end{bmatrix}.$ (105)

Therefore the intercluster eigenvalues are the eigenvalues of the Jacobian
$\mathcal{J}_{M}$ in (105) which is also given in (63). One of these is zero
corresponding to pure rotations, observing that the first and third columns of
$\mathcal{J}_{M}$ are identical. The others correspond to changes in the phase
difference between clusters and changes in the isostable values for the two
synchronised clusters. In conclusion,

$\displaystyle\begin{split}\det(\mathcal{J}-\mu
I_{2N})=&\det(\mathcal{S}^{A}-\mu I_{2})^{N_{A}-1}\det(\mathcal{S}^{B}-\mu
I_{2})^{N_{B}-1}\\\ &\qquad\times\det(\mathcal{J}_{M}-\mu I_{4}).\end{split}$
(106)

##  Appendix C Higher order phase equation derivation

Here we follow and correct Park and Wilson Park and Wilson (2021) to derive
the higher order phase equation (112) explicitly including terms up to cubic
order in the interaction strength $\epsilon$. Park and Wilson Park and Wilson
(2021) take the isostable coordinate as an $O(\epsilon)$ term which may be
expressed as $\psi(t)=\epsilon p^{(1)}(t)+\epsilon^{2}p^{(2)}(t)+\cdots$ where
$p^{(k)}(t)$ are $O(1)$. Using this additional expansion in (10b) and matching
terms at different powers of $\epsilon$ results in a hierarchy of linear first
order differential equations for the $p^{(k)}(t)$ Park and Wilson (2021). The
equations at $O(\epsilon)$ and $O(\epsilon^{2})$ are

$\displaystyle\frac{{\rm d}p_{i}^{(1)}}{{\rm d}t}$ $\displaystyle=\kappa
p_{i}^{(1)}+\sum_{j=1}^{N}w_{ij}h_{4}(\theta_{i},\theta_{j})$ (107)
$\displaystyle\frac{{\rm d}p_{i}^{(2)}}{{\rm d}t}$ $\displaystyle=\kappa
p_{i}^{(2)}+\sum_{j=1}^{N}w_{ij}\left[p_{i}^{(1)}h_{5}(\theta_{i},\theta_{j})+p_{j}^{(1)}h_{6}(\theta_{i},\theta_{j})\right],$
(108)

where $h_{4},h_{5},h_{6}$ are given by (19d)–(19f). The forcing terms of the
equation at a given order depend only on the solutions of lower order
equations and hence each can be solved explicitly in turn using the
integrating factor $e^{-\kappa t}$. Following Park and Wilson (2021); Wilson
and Ermentrout (2019b) (and also correcting notational errors in Park and
Wilson (2021)) we observe that

$\displaystyle p_{i}^{(1)}(t)$
$\displaystyle=\sum_{k=1}^{N}w_{ik}q^{(1)}(\theta_{i},\theta_{k})$ (109)
$\displaystyle p_{i}^{(2)}(t)$
$\displaystyle=\sum_{k=1}^{N}\sum_{l=1}^{N}\left[w_{ik}w_{il}q^{(2)}(\theta_{i},\theta_{k},\theta_{l})+w_{ik}w_{kl}q^{(3)}(\theta_{i},\theta_{k},\theta_{l})\right],$
(110)

where

$\displaystyle q^{(1)}(\theta_{i},\theta_{k})$
$\displaystyle=\int_{0}^{\infty}e^{\kappa s}h_{4}(\theta_{i}-\omega
s,\theta_{k}-\omega s)\ {\rm d}s,$ (111a)
$\displaystyle\begin{split}q^{(2)}(\theta_{i},\theta_{k},\theta_{l})&=\int_{0}^{\infty}e^{\kappa
s}q^{(1)}(\theta_{i}-\omega s,\theta_{l}-\omega s)\\\ &\qquad\times
h_{5}(\theta_{i}-\omega s,\theta_{k}-\omega s)\ {\rm d}s,\end{split}$ (111b)
$\displaystyle\begin{split}q^{(3)}(\theta_{i},\theta_{k},\theta_{l})&=\int_{0}^{\infty}e^{\kappa
s}q^{(1)}(\theta_{k}-\omega s,\theta_{l}-\omega s)\\\ &\qquad\times
h_{6}(\theta_{i}-\omega s,\theta_{k}-\omega s)\ {\rm d}s.\end{split}$ (111c)

Higher order terms $p^{(k)}$, $k>2$ may also be determined (see Park and
Wilson (2021)) but we do not require these here. Substituting all expansions
into (10a) we obtain the higher order phase equation

$\displaystyle\begin{split}\frac{{\rm d}\theta_{i}}{{\rm d}t}=\omega\ +\
&\epsilon\sum_{j=1}^{N}w_{ij}\overline{h}_{1}(\theta_{i},\theta_{j})+\epsilon^{2}\sum_{j,k=1}^{N}\left[w_{ij}w_{ik}\overline{h}_{2}(\theta_{i},\theta_{j},\theta_{k})+w_{ij}w_{jk}\overline{h}_{3}(\theta_{i},\theta_{j},\theta_{k})\right]\\\
+\
&\epsilon^{3}\sum_{j,k,l=1}^{N}\Bigl{[}w_{ij}w_{ik}w_{il}\overline{h}_{4}(\theta_{i},\theta_{j},\theta_{k},\theta_{l})+w_{ij}w_{ik}w_{kl}\overline{h}_{5}(\theta_{i},\theta_{j},\theta_{k},\theta_{l})+w_{ij}w_{jk}w_{jl}\overline{h}_{6}(\theta_{i},\theta_{j},\theta_{k},\theta_{l})\\\
&\qquad\qquad\qquad\qquad+w_{ij}w_{jk}w_{kl}\overline{h}_{7}(\theta_{i},\theta_{j},\theta_{k},\theta_{l})+w_{ij}w_{ik}w_{jl}\overline{h}_{8}(\theta_{i},\theta_{j},\theta_{k},\theta_{l})\Bigr{]},\end{split}$
(112)

where

$\displaystyle\overline{h}_{1}(\theta_{i},\theta_{j})$
$\displaystyle=h_{1}(\theta_{i},\theta_{j}),$ (113a)
$\displaystyle\overline{h}_{2}(\theta_{i},\theta_{j},\theta_{k})$
$\displaystyle=q^{(1)}(\theta_{i},\theta_{k})h_{2}(\theta_{i},\theta_{j}),$
(113b) $\displaystyle\overline{h}_{3}(\theta_{i},\theta_{j},\theta_{k})$
$\displaystyle=q^{(1)}(\theta_{j},\theta_{k})h_{3}(\theta_{i},\theta_{j}),$
(113c)
$\displaystyle\begin{split}\overline{h}_{4}(\theta_{i},\theta_{j},\theta_{k},\theta_{l})&=q^{(2)}(\theta_{i},\theta_{k},\theta_{l})h_{2}(\theta_{i},\theta_{j})\\\
&\quad+q^{(1)}(\theta_{i},\theta_{k})q^{(1)}(\theta_{i},\theta_{l})h_{7}(\theta_{i},\theta_{j}),\end{split}$
(113d)
$\displaystyle\overline{h}_{5}(\theta_{i},\theta_{j},\theta_{k},\theta_{l})$
$\displaystyle=q^{(3)}(\theta_{i},\theta_{k},\theta_{l})h_{2}(\theta_{i},\theta_{j}),$
(113e)
$\displaystyle\begin{split}\overline{h}_{6}(\theta_{i},\theta_{j},\theta_{k},\theta_{l})&=q^{(2)}(\theta_{j},\theta_{k},\theta_{l})h_{3}(\theta_{i},\theta_{j})\\\
&\quad+q^{(1)}(\theta_{j},\theta_{k})q^{(1)}(\theta_{j},\theta_{l})h_{8}(\theta_{i},\theta_{j}),\end{split}$
(113f)
$\displaystyle\overline{h}_{7}(\theta_{i},\theta_{j},\theta_{k},\theta_{l})$
$\displaystyle=q^{(3)}(\theta_{j},\theta_{k},\theta_{l})h_{3}(\theta_{i},\theta_{j}),$
(113g)
$\displaystyle\overline{h}_{8}(\theta_{i},\theta_{j},\theta_{k},\theta_{l})$
$\displaystyle=q^{(1)}(\theta_{i},\theta_{k})q^{(1)}(\theta_{j},\theta_{l})h_{9}(\theta_{i},\theta_{j}),$
(113h)

for $h_{1},\ldots,h_{6}$ as defined in (19), and for

$\displaystyle\begin{split}h_{7}(\theta_{i},\theta_{j})&=Z^{(0)}(\theta_{i})\cdot
K_{1}(\theta_{i},\theta_{j})+Z^{(1)}(\theta_{i})\cdot
J_{1}g^{(1)}(\theta_{i})\\\
&\quad+Z^{(2)}(\theta_{i})\cdot\mathbf{G}(\mathbf{x}^{\gamma}(\theta_{i}),\mathbf{x}^{\gamma}(\theta_{j})),\end{split}$
(114a) $\displaystyle h_{8}(\theta_{i},\theta_{j})$
$\displaystyle=Z^{(0)}(\theta_{i})\cdot K_{2}(\theta_{i},\theta_{j}),$ (114b)
$\displaystyle h_{9}(\theta_{i},\theta_{j})$
$\displaystyle=Z^{(0)}(\theta_{i})\cdot
L(\theta_{i},\theta_{j})+Z^{(1)}(\theta_{i})\cdot J_{2}g^{(1)}(\theta_{j}).$
(114c)

Moving to a rotating reference frame, applying first-order averaging and
returning to the original variables we can obtain the autonomous approximation
to the higher order phase reduced equation (64) where

$\displaystyle\overline{H}_{1}(\chi)$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}\overline{h}_{1}(u,u+\chi)\,{\rm
d}u,$ (115a) $\displaystyle\overline{H}_{m}(\chi,\eta)$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}\overline{h}_{m}(u,u+\chi,u+\eta)\,{\rm
d}u,\quad m=2,3,$ (115b)
$\displaystyle\begin{split}\overline{H}_{m}(\chi,\eta,\xi)&=\frac{1}{2\pi}\int_{0}^{2\pi}\overline{h}_{m}(u,u+\chi,u+\eta,u+\xi)\,{\rm
d}u,\\\ &\qquad\qquad\qquad\qquad m=4,\ldots,8.\end{split}$ (115c)

##  Appendix D Splay state linear stability boundary for $N\geq 3$

In the large $N$ limit, from (50) we calculate

$\displaystyle\lambda_{q}^{(1)}=$ $\displaystyle\begin{cases}-\epsilon
i\left(i(1+c_{1}c_{2})+c_{2}-c_{1}\right)/2&\text{for }q=1,\\\
0&\text{otherwise, }\end{cases}$ (116) $\displaystyle\lambda_{q}^{(2)}=$
$\displaystyle\begin{cases}A(i-c_{1})(1+c_{2}^{2})\epsilon/2&\text{for
}q=1,\\\ 0&\text{otherwise, }\end{cases}$ (117)
$\displaystyle\lambda_{q}^{(3)}=$ $\displaystyle\begin{cases}\frac{\epsilon
i(ic_{1}+1)}{2A(\epsilon-1)}&\text{for }q=1,\\\ 0&\text{otherwise,
}\end{cases}$ (118) $\displaystyle\lambda_{q}^{(4)}=$
$\displaystyle\begin{cases}\kappa+\epsilon\left((i-c_{1})c_{2}+ic_{1}+5\right)/2&\text{for
}q=1,\\\ 2(\epsilon-1)&\text{otherwise. }\end{cases}$ (119)

Therefore using (47), the eigenvalues of the Jacobian for the splay state are
$0$ and $2(\epsilon-1)$, each of multiplicity $N-1$, and the eigenvalues of
$\Lambda_{1}$. It can be shown that there is a pair of eigenvalues $\mu=\pm
i\sigma$ with

$\displaystyle\sigma=\frac{\epsilon\left(c_{2}^{2}c_{1}-2c_{2}+3c_{1})\epsilon^{2}+2(c_{2}-c1)(2\epsilon-1)\right)}{2(\epsilon-1)(3\epsilon-2)},$
(120)

when $\epsilon=\epsilon_{0,PI}$ satisfies (69). Equation (69) defines a
stability boundary for the splay state in the large $N$ limit. It can
similarly be shown that for finite $N\geq 3$ the eigenvalues of the Jacobian
are $0$, $2(\epsilon-1)$ each of multiplicity $N-2$ and the remaining
eigenvalues are the eigenvalues of $\Lambda_{1}$ and their complex conjugates.
This gives the identical equation (69) for the stability boundary.

##  Appendix E The Morris-Lecar model

The Morris-Lecar model is a planar model of neural activity Morris and Lecar
(1981) with two ionic currents; an outward non-inactivating potassium current
and an inward, non-inactivating calcium current. Assuming that the calcium
dynamics operate on a much faster time scale than the potassium, the dynamics
of an isolated node are given by

$\displaystyle\begin{split}C_{m}\frac{{\rm d}v_{i}}{{\rm
d}t}&=I_{b}-g_{L}(v_{i}-E_{L})-g_{K}w_{i}(v_{i}-E_{K})\\\ &\quad-
g_{Ca}m_{\infty}(v_{i})(v_{i}-E_{Ca})\end{split}$ (121a)
$\displaystyle\frac{{\rm d}w_{i}}{{\rm d}t}$
$\displaystyle=\phi(w_{\infty}(v_{i})-w_{i})\lambda(v_{i}),$ (121b)

where

$\displaystyle m_{\infty}(v)$
$\displaystyle=0.5\left(1+\tanh\left(\left(v-V_{1}\right)/V_{2}\right)\right)$
$\displaystyle w_{\infty}(v)$
$\displaystyle=0.5\left(1+\tanh\left(\left(v-V_{3}\right)/V_{4}\right)\right)$
$\displaystyle\lambda(v)$
$\displaystyle=\cosh\left(\left(v-V_{3}\right)/\left(2V_{4}\right)\right).$

and $V_{1},\ldots,V_{4}$ and $\phi$ are constants. Here $w_{i}$ represents the
fraction of K+ channels open at node $i$ and the Ca2+ channels respond to
$v_{i}$ so rapidly that we assume instant activation. Here $g_{L}$ is the
leakage conductance, $g_{K}$ and $g_{Ca}$ are the potassium and calcium
conductances, $E_{L}$, $E_{K}$, $E_{Ca}$ are corresponding reversal
potentials, $m_{\infty}(v)$, $w_{\infty}(v)$ are voltage dependent gating
functions, $\lambda(v)$ is a voltage-dependent rate, $I_{b}$ is the applied
current and $C_{m}$ is the cell capacitance. Unless otherwise stated, the
default model parameters used throughout are as in Han, Kurrer, and Kuramoto
(1995): $\phi=1.15$, $g_{Ca}=1$, $g_{K}=2$, $g_{L}=0.5$, $E_{Ca}=1$,
$E_{K}=-0.7$, $E_{L}=-0.5$, $V_{1}=-0.01$, $V_{2}=0.15$, $V_{3}=0.1$,
$V_{4}=0.145$ and $C_{m}=1$. For these parameter values a limit cycle of the
system arises at $I_{b}=0.0730$ through a homoclinic connection. We choose
$I_{b}=0.075$ to place the system near to the homoclinic bifurcation. For
these parameter values we find that the periodic orbit has Floquet exponent
$\kappa=-0.4094$ and period $T=8.1654$.

## References

  * Park and Wilson (2021) Y. Park and D. D. Wilson, “High-order accuracy computation of coupling functions for strongly coupled oscillators,” SIAM Journal on Applied Dynamical Systems 20, 1464–1484 (2021).
  * León and Pazó (2019) I. León and D. Pazó, “Phase-reduction beyond the first order: The case of the mean-field complex Ginzburg-Landau equation,” Physical Review E 100, 012211 (2019).
  * Kuramoto (1984) Y. Kuramoto,  _Chemical Oscillations, Waves, and Turbulence_ , Springer Series in Synergetics (Springer-Verlag, Berlin, 1984).
  * Winfree (2001) A. Winfree,  _The Geometry of Biological Time_ , 2nd ed., Springer Study Edition (Springer-Verlag, New York, 2001).
  * Ermentrout and Terman (2010) G. B. Ermentrout and D. H. Terman,  _Mathematical Foundations of Neuroscience_ , Interdisciplinary Applied Mathematics, Vol. 35 (Springer, New York, 2010).
  * Brown, Moehlis, and Holmes (2004) E. Brown, J. Moehlis, and P. Holmes, “On the phase reduction and response dynamics of neural oscillator populations,” Neural Computation 16, 673–715 (2004).
  * Izhikevich (2007) E. Izhikevich,  _Dynamical systems in neuroscience: the geometry of excitability and bursting_ (MIT Press, London, 2007).
  * Pietras and Daffertshofer (2019) B. Pietras and A. Daffertshofer, “Network dynamics of coupled oscillators and phase reduction techniques,” Physics Reports 819, 1–105 (2019), network dynamics of coupled oscillators and phase reduction techniques.
  * Nakao (2016) H. Nakao, “Phase reduction approach to synchronisation of nonlinear oscillators,” Contemporary Physics 57, 188–214 (2016).
  * Monga  _et al._ (2019) B. Monga, D. Wilson, T. Matchen, and J. Moehlis, “Phase reduction and phase-based optimal control for biological systems: a tutorial,” Biological Cybernetics 113, 11–46 (2019).
  * Winfree (1967) A. T. Winfree, “Biological rhythms and the behavior of populations of coupled oscillators,” Journal of Theoretical Biology 16, 15–42 (1967).
  * Guckenheimer (1975) J. Guckenheimer, “Isochrons and phaseless sets,” Journal of Mathematical Biology 1, 259–273 (1975).
  * Dörfler, Chertkov, and Bullo (2013) F. Dörfler, M. Chertkov, and F. Bullo, “Synchronization in complex oscillator networks and smart grids,” Proceedings of the National Academy of Sciences 110, 2005–2010 (2013).
  * Dörfler and Bullo (2014) F. Dörfler and F. Bullo, “Synchronization in complex networks of phase oscillators: A survey,” Automatica 50, 1539–1564 (2014).
  * Abrams and Strogatz (2004) D. M. Abrams and S. H. Strogatz, “Chimera states for coupled oscillators,” Phys. Rev. Lett. 93, 174102 (2004).
  * Hoppensteadt and Izhikevich (1997) F. Hoppensteadt and E. Izhikevich,  _Weakly Connected Neural Networks_ , Applied Mathematical Sciences (NY: Springer-Verlag, 1997).
  * Wedgwood  _et al._ (2013) K. C. A. Wedgwood, K. K. Lin, R. Thul, and S. Coombes, “Phase-amplitude descriptions of neural oscillator models,” Journal of Mathematical Neuroscience 3, 2 (2013).
  * Letson and Rubin (2018) B. Letson and J. E. Rubin, “A new frame for an old (phase) portrait: Finding rivers and other flow features in the plane,” SIAM Journal on Applied Dynamical Systems 17, 2414–2445 (2018).
  * Guillamon and Huguet (2009) A. Guillamon and G. Huguet, “A computational and geometric approach to phase resetting curves and surfaces,” SIAM Journal on Applied Dynamical Systems 8, 1005–1042 (2009).
  * Mauroy, Mezić, and Moehlis (2013) A. Mauroy, I. Mezić, and J. Moehlis, “Isostables, isochrons, and koopman spectrum for the action–angle representation of stable fixed point dynamics,” Physica D: Nonlinear Phenomena 261, 19–30 (2013).
  * Wilson and Moehlis (2016) D. Wilson and J. Moehlis, “Isostable reduction of periodic orbits,” Phys. Rev. E 94, 052213 (2016).
  * Budišić, Mohr, and Mezić (2012) M. Budišić, R. Mohr, and I. Mezić, “Applied koopmanism,” Chaos 22, 047510 (2012).
  * Castejón, Guillamon, and Huguet (2013) O. Castejón, A. Guillamon, and G. Huguet, “Phase-amplitude response functions for transient-state stimuli,” Journal of Mathematical Neuroscience 3, 13 (2013).
  * Pérez-Cervera, M-Seara, and Huguet (2020) A. Pérez-Cervera, T. M-Seara, and G. Huguet, “Global phase-amplitude description of oscillatory dynamics via the parameterization method,” Chaos 30, 083117 (2020).
  * Mauroy and Mezić (2018) A. Mauroy and I. Mezić, “Global computation of phase-amplitude reduction for limit-cycle dynamics,” Chaos 28, 073108 (2018).
  * Shirasaka, Kurebayashi, and Nakao (2017) S. Shirasaka, W. Kurebayashi, and H. Nakao, “Phase-amplitude reduction of transient dynamics far from attractors for limit-cycling systems,” Chaos 27, 023119 (2017).
  * Wilson and Ermentrout (2018) D. Wilson and B. Ermentrout, “Greater accuracy and broadened applicability of phase reduction using isostable coordinates,” Journal of Mathematical Biology 76, 37–66 (2018).
  * Wilson and Moehlis (2022) D. Wilson and J. Moehlis, “Recent advances in the analysis and control of large populations of neural oscillators,” Annual Reviews in Control 54, 327–351 (2022).
  * Wilson and Ermentrout (2019a) D. Wilson and B. Ermentrout, “Augmented phase reduction of (not so) weakly perturbed coupled oscillators,” SIAM Review 61, 277–315 (2019a).
  * Wilson (2019) D. Wilson, “Isostable reduction of oscillators with piecewise smooth dynamics and complex floquet multipliers,” Phys. Rev. E 99, 022210 (2019).
  * Wilson (2020a) D. Wilson, “Phase-amplitude reduction far beyond the weakly perturbed paradigm,” Phys. Rev. E 101, 022220 (2020a).
  * Wilson and Ermentrout (2019b) D. Wilson and B. Ermentrout, “Phase models beyond weak coupling,” Phys. Rev. Lett. 123, 164101 (2019b).
  * Ermentrout, Park, and Wilson (2019) B. Ermentrout, Y. Park, and D. Wilson, “Recent advances in coupled oscillator theory,” Philosophical Transactions of the Royal Society A 377, 20190092 (2019).
  * Coombes  _et al._ (2023) S. Coombes, M. Şayli, R. Thul, R. Nicks, M. A. Porter, and Y. M. Lai, “Oscillatory networks: Insights from piecewise-linear modelling,”  (2023).
  * Park and Wilson (2023) Y. Park and D. Wilson, “N-body oscillator interactions of higher-order coupling functions,”  (2023).
  * Mau, Rosenblum, and Pikovsky (2023) E. T. K. Mau, M. Rosenblum, and A. Pikovsky, “High-order phase reduction for coupled 2d oscillators,”  (2023).
  * Battiston  _et al._ (2020) F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, and G. Petri, “Networks beyond pairwise interactions: Structure and dynamics,” Physics Reports 874, 1–92 (2020), networks beyond pairwise interactions: Structure and dynamics.
  * Pecora and Carroll (1998) L. M. Pecora and T. L. Carroll, “Master stability functions for synchronized coupled systems,” Phys. Rev. Lett. 80, 2109–2112 (1998).
  * Morris and Lecar (1981) C. Morris and H. Lecar, “Voltage oscillations in the barnacle giant muscle fibre,” Biophysical Journal 35, 193–213 (1981).
  * Han, Kurrer, and Kuramoto (1995) S. K. Han, C. Kurrer, and Y. Kuramoto, “Dephasing and Bursting in Coupled Neural Oscillators,” Physical Review Letters 75, 3190–3193 (1995).
  * Ermentrout and Kopell (1991) G. Ermentrout and N. Kopell, “Multiple pulse interactions and averaging in systems of coupled neural oscillators,” Journal of Mathemarical Biology 29, 195–217 (1991).
  * Guckenheimer and Holmes (1990) J. Guckenheimer and P. Holmes,  _Nonlinear oscillations, dynamical systems, and bifurcations of vector fields_ , Applied Mathematical Sciences, Vol. 42 (Springer-Verlag, 1990).
  * Ashwin and Swift (1992) P. Ashwin and J. Swift, “The dynamics of n weakly coupled identical oscillators.” Journal of Nonlinear Science 2, 69–108 (1992).
  * Kim and Lee (2000) S. Kim and S. G. Lee, “Phase dynamics in the biological neural networks,” Physica A: Statistical Mechanics and its Applications 288, 380–396 (2000), dynamics Days Asia-Pacific: First International Conference on NonLinear Science.
  * Watanabe and Swift (1997) S. Watanabe and J. Swift, “Stability of periodic solutions in series arrays of josephson junctions with internal capacitance,” Journal of Nonlinear Science 7, 503–536 (1997).
  * Ermentrout (1992) G. B. Ermentrout, “Stable periodic solutions to discrete and continuum arrays of weakly coupled nonlinear oscillators,” SIAM Journal on Applied Mathematics 52, 1665–1687 (1992).
  * Coombes and Wedgwood (2023) S. Coombes and K. C. A. Wedgwood,  _Neurodynamics: An Applied Mathematics Perspective_ , Texts in Applied Mathematics, Vol. 75 (Springer, 2023).
  * Kuramoto and Nakao (2019) Y. Kuramoto and H. Nakao, “On the concept of dynamical reduction: the case of coupled oscillators,” Philosophical Transactions of the Royal Society A 377, 2160 (2019).
  * Cumin and Unsworth (2007) D. Cumin and C. Unsworth, “Generalising the kuramoto model for the study of neuronal synchronisation in the brain,” Physica D: Nonlinear Phenomena 226, 181–196 (2007).
  * Forrester (2015) D. Forrester, “On the concept of dynamical reduction: the case of coupled oscillators,” Scientific Reports 5, 16994 (2015).
  * Wiesenfeld, Colet, and Strogatz (1998) K. Wiesenfeld, P. Colet, and S. H. Strogatz, “Frequency locking in josephson arrays: Connection with the kuramoto model,” Phys. Rev. E 57, 1563–1569 (1998).
  * Kvalheim and Revzen (2021) M. D. Kvalheim and S. Revzen, “Existence and uniqueness of global koopman eigenfunctions for stable fixed points and periodic orbits,” Physica D: Nonlinear Phenomena 425, 132959 (2021).
  * Sanders, Verhulst, and Murdock (2007) J. A. Sanders, F. Verhulst, and J. Murdock,  _Averaging methods in nonlinear dynamical systems, 2nd edn._ (New York, NY: Springer, 2007).
  * Ashwin, Coombes, and Nicks (2016) P. Ashwin, S. Coombes, and R. Nicks, “Mathematical Frameworks for Oscillatory Network Dynamics in Neuroscience,” Journal of Mathematical Neuroscience 2, 6 (2016).
  * Silvester (2000) J. R. Silvester, “Determinants of block matrices,” The Mathematical Gazette 84, 460 – 467 (2000).
  * Tee (2007) G. J. Tee, “Eigenvectors of block circulant and alternating circulant matrices,” New Zealand Journal of Mathematics 36, 195–211 (2007).
  * Kra and Simanca (2012) I. Kra and S. R. Simanca, “On circulant matrices,” Notices of the AMS 59, 368–377 (2012).
  * Ding and Zhou (2007) J. Ding and A. Zhou, “Eigenvalues of rank-one updated matrices with some applications,” Applied Mathematics Letters 20, 1223–1226 (2007).
  * Brown, Holmes, and Moehlis (2003) E. Brown, P. Holmes, and J. Moehlis, “Globally coupled oscillator networks,” in  _Perspectives and Problems in Nolinear Science: A Celebratory Volume in Honor of Lawrence Sirovich_ , edited by E. Kaplan, J. E. Marsden, and K. R. Sreenivasan (Springer New York, New York, NY, 2003) pp. 183–215.
  * Hakim and Rappel (1992) V. Hakim and W.-J. Rappel, “Dynamics of the globally coupled complex ginzburg-landau equation,” Phys. Rev. A 46, R7347–R7350 (1992).
  * Nakagawa and Kuramoto (1993) N. Nakagawa and Y. Kuramoto, “Collective Chaos in a Population of Globally Coupled Oscillators,” Progress of Theoretical Physics 89, 313–323 (1993).
  * Ermentrout (2002) G. B. Ermentrout,  _Simulating, analyzing, and animating dynamical systems - a guide to XPPAUT for researchers and students_ (SIAM, Philadelphia, 2002).
  * Mauroy and Mezić (2012) A. Mauroy and I. Mezić, “On the use of Fourier averages to compute the global isochrons of (quasi)periodic dynamics,” Chaos: An Interdisciplinary Journal of Nonlinear Science 22 (2012), 10.1063/1.4736859, 033112.
  * Bick, Böhle, and Kuehn (2023) C. Bick, T. Böhle, and C. Kuehn, “Higher-order interactions in phase oscillator networks through phase reductions of oscillators with phase dependent amplitude,”  (2023).
  * Wilson (2022) D. Wilson, “An adaptive phase-amplitude reduction framework without $\mathcal{O}(\epsilon)$ constraints on inputs,” SIAM Journal on Applied Dynamical Systems 21, 204–230 (2022).
  * Wilson (2020b) D. Wilson, “Optimal open-loop desynchronization of neural oscillator populations,” Journal of Mathematical Biology 81, 25–64 (2020b).
  * Nicks, Chambon, and Coombes (2018) R. Nicks, L. Chambon, and S. Coombes, “Clusters in nonsmooth oscillator networks,” Physical Review E 97, 032213 (2018).
  * Medvedev (2019) G. S. Medvedev, “The continuum limit of the Kuramoto model on sparse random graphs,” Communications in Mathematical Sciences 17, 883–898 (2019).
  * Ihara and Yagasaki (2023) R. Ihara and K. Yagasaki, “Continuum Limits of Coupled Oscillator Networks Depending on Multiple Sparse Graphs,” Journal of Nonlinear Science 33, 62 (2023).
  * Kotani  _et al._ (2012) K. Kotani, I. Yamaguchi, Y. Ogawa, Y. Jimbo, H. Nakao, and G. B. Ermentrout, “Adjoint Method Provides Phase Response Functions for Delay-Induced Oscillations,” Physical Review Letters 109, 044101 (2012).
  * Novičenko and Pyragas (2012) V. Novičenko and K. Pyragas, “Phase reduction of weakly perturbed limit cycle oscillations in time-delay systems,” Physica D: Nonlinear Phenomena 241, 1090–1098 (2012).
  * Magnus and Neudecker (2007) J. R. Magnus and H. Neudecker,  _Matrix Differential Calculus with Applications in Statistics and Econometrics_ (John Wiley and Sons, Chichester, 2007).