# Emergent hypernetworks in weakly coupled oscillators

Eddie Nijholt1, Jorge Luis Ocampo-Espindola2, Deniz Eroglu3,
István Z. Kiss2, Tiago Pereira1,4∗

1Instituto de Ciências Matemáticas e Computação, Universidade de São Paulo,
São Carlos, Brazil
2Department of Chemistry, Saint Louis University, St. Louis, USA
3Faculty of Engineering and Natural Sciences, Kadir Has University, Istanbul,
Turkey
4Department of Mathematics, Imperial College London, London, UK

∗To whom correspondence should be addressed; E-mail:
<EMAIL_ADDRESS>
###### Abstract

Networks of weakly coupled oscillators had a profound impact on our
understanding of complex systems. Studies on model reconstruction from data
have shown prevalent contributions from hypernetworks with triplet and higher
interactions among oscillators, in spite that such models were originally
defined as oscillator networks with pairwise interactions. Here, we show that
hypernetworks can spontaneously emerge even in the presence of pairwise albeit
nonlinear coupling given certain triplet frequency resonance conditions. The
results are demonstrated in experiments with electrochemical oscillators and
in simulations with integrate-and-fire neurons. By developing a comprehensive
theory, we uncover the mechanism for emergent hypernetworks by identifying
appearing and forbidden frequency resonant conditions. Furthermore, it is
shown that microscopic linear (difference) coupling among units results in
coupled mean fields, which have sufficient nonlinearity to facilitate
hypernetworks. Our findings shed light on the apparent abundance of
hypernetworks and provide a constructive way to predict and engineer their
emergence.

Published in Nature Communications 13, 4849 (2022)

## Introduction

Networks of weakly coupled oscillators are prolific models for a variety of
natural systems ranging from biology [1, 2] and chemistry [3, 4] to
neuroscience [5, 6] via ecology [7] to engineering [8]. Such networks serve as
stepping stones to understand collective dynamics [9, 10, 11, 12] and other
emergent phenomena in networks [13, 14]. In these models, the interactions are
described in a pairwise manner and the collective dynamics of a network can be
predicted by the superposition of such pairwise interactions.

Recent work, however, suggests that many networks described as pairwise
interactions can be better described in terms of hypernetworks with triplet
and quadruplet interactions among nodes [15, 16, 17, 18]. In fact,
hypernetworks appear as suitable representations of certain dynamical
processes found in physics [19, 20], chemistry [21] and neuroscience [22, 23].
This has ignited research aimed at understanding the impact of higher-order
interactions on the dynamical behavior of complex systems [24, 25, 26, 27].
Moreover, besides considering hypernetworks as a good description of such
models, we observed that hypernetworks could be revealed in data-driven model
reconstructions when the original model is a network. Therefore, a major
puzzle is why hypernetworks emerge as the fitting description of actual
network data.

Here, we show that hypernetworks can describe experimental data of networks of
electrochemical oscillators with nonlinear coupling. We uncover a mechanism
that generates higher-order interactions as a model to describe oscillator
networks from data. First, we show that sparse model recovery from data
reveals higher-order interactions. We then develop a theory for the emergence
of such higher-order interactions when the isolated system is close to a Hopf
bifurcation. We provide an algorithm to reveal emergent hypernetwork and its
emergent coupling functions for any network in disciplines ranging from
neuroscience to chemistry.  The emergent hypernetworks provide a dimension
reduction that allows the characterization of critical transitions.

## Results

### Emergent hypernetworks in electrochemical experiments

We designed an experimental system with four oscillatory chemical reactions
coupled with nonlinear feedback and delay arranged in a ring network (see Fig.
1 (a)). The set-up consists of a multichannel potentiostat interfaced with a
real-time controller and connected to a Pt counter, a Hg/Hg2SO4 sat K2SO4
reference, and four Ni working electrodes in 3.0 M sulfuric acid electrolyte.
At a constant circuit potential ($V_{0}$=1100 mV with respect to the reference
electrode) and with an external resistance ($R_{ind}$=1.0 kohm) attached to
each nickel wire, the electrochemical dissolution of nickel exhibits periodic
current and electrode potential oscillations with a natural frequency of 0.385
Hz.

Without coupling, we adjusted the natural frequency of each oscillator to have
a ratio with respect to oscillator 1 as $\omega_{2}/\omega_{1}$ = 2.53
($\approx 2.5$), $\omega_{3}/\omega_{1}$ =1.56 ($\approx 1.5$) and
$\omega_{4}/\omega_{1}$ = 2.53 ($\approx 2.5$) with a set of resistors and
capacitors ($C_{ind}$), see Supplementary Note 1.) The natural frequencies
create opportunities for triplet resonances, as there are small detunings for
$\omega_{1}-\omega_{2}+\omega_{3}$ and $\omega_{1}-\omega_{4}+\omega_{3}$,  as
well as pairwise resonances $\omega_{2}\approx\omega_{4}$.

The individual electrode potentials ($E_{k}$) were recorded and rescaled and
offset corrected

$\tilde{E}_{k}=O_{k}[E_{k}-o_{k}],$ (1)

where $o_{k}$ and $O_{k}$ are the time-averaged electrode potential and
amplitude rescaling factor, respectively. (The rescaling factors,
$O_{k}={0.5,1,0.5,1}$ V were applied to counter the different amplitudes of
the slow oscillators.) A ring-coupling can be introduced with external
feedback (see Fig. 1 (b,c)) according to

$\displaystyle
V_{k}(t)=V_{0,k}+K\sum_{\ell=1}^{4}A_{k\ell}h[\tilde{E}_{k}(t),\tilde{E}_{\ell}(t-\tau)],$
(2)

where $V_{k}(t)$ and $V_{0,k}$ are the applied and the offset circuit
potential of the $k$th electrode, respectively, $K$ is the coupling strength,
$A_{k\ell}$ is the adjacency matrix, $\tau$ is a time delay, and

$\displaystyle
h[\tilde{E}_{k}(t),\tilde{E}_{\ell}(t-\tau)]=(\tilde{E}_{k}(t)+\tilde{E}_{k}(t)^{2})\tilde{E}_{\ell}(t-\tau).$
(3)

This delayed nonlinear feedback modulates the impact of the coupled units with
a bias towards positive values (similar to a diode operation in the $(-1,1)$
interval). Note that this form of feedback is fundamentally different from
previously applied nonlinear schemes [4] in that it does not produce obvious
synchronization patterns, for example, one and multi-cluster states.

Figure 1 (d) shows the time series of the electrode potential for $K$=5.2 and
$\tau$ =1.65 s. The slow oscillators (1 and 3) have larger amplitudes and the
time series exhibit nonlinear waveform modulations without any obvious
synchronization pattern (one-cluster state).

Figure 1: Emergent hypernetworks in an electrochemical network experiment. (a)
Experimental setup. (b) Schematic illustration of the electrochemical
experiment with the nonlinear feedback. The blue, orange, yellow, and green
lines represent the elements 1 to 4, respectively. The electrode potential
signals ($E_{k}$) of the four (nearly) isolated electrodes are nonlinearly
modulated and fed back with a delay $\tau$ to the corresponding circuit
potential ($V_{k}$), which drives the metal dissolution. (The delay is
implemented by storing the past data in the memory of the computer.)  (c)
Representation of the in a ring network topology used in the experiment. (d)
Electrode potential time series. (e) Filtered and fitted (dark red line)
instantaneous frequency using LASSO for hypernetwork reconstruction
corresponding from top to bottom to oscillators 1 to 4, respectively. (f)
Experimental recovery of the phase interactions given by a hypernetwork.

From the potentials $\tilde{E}_{k}$ we extract the frequencies
$\dot{\theta}_{k}$ and apply a first-order Savitzky-Golay filter with a time
window of 45 s to remove the in-cycle and short-range phase fluctuation, as
shown in Figure 1 (e) (solid line). For each oscillator, a slow variation is
seen as the oscillators slow down and speed up on a timescale of about 100 s
(or 40 cycles); notably, the elements 1 and 3 exhibit similar
$\dot{\theta}_{k}$ oscillations, which are different from those in elements 2
and 4.

To describe the nature of the phase dynamics, we consider the slow triplet
phase differences

$\displaystyle\begin{split}\phi_{1}&=\theta_{1}-\theta_{2}+\theta_{3}\\\
\phi_{2}&=\theta_{1}-\theta_{4}+\theta_{3},\end{split}$ (4)

which correspond to the triplet frequency detunings.

The impact of triplet interactions on the dynamics can be extracted with a
LASSO fit to

$\dot{\theta}_{k}=\hat{\omega}_{k}(t)+\sum_{j=1}^{2}C_{j}^{k}\sin(\phi_{j})+D_{j}^{k}\cos(\phi_{j})$
(5)

where
$\hat{\omega}_{k}(t)=\hat{\omega}^{0}_{k}+\hat{\omega}^{1}_{k}t+\hat{\omega}^{2}_{k}t^{2}$
is the fitted, slowly drifting (up to quadratic variation in time) natural
frequency, and $C_{j}^{k}$ and $D_{j}^{k}$ are the amplitudes of the sin and
cos phase coupling functions corresponding to the appropriate triplet phase
differences. The strength of the triplet interactions $j=1,2$ (for $\phi_{j}$)
on oscillator $k$ is given by the amplitudes
$H_{j}^{k}=\sqrt{(C^{k}_{j})^{2}+(D^{k}_{j})^{2}}$.

The dynamics of oscillators 1 and 3 are impacted by both triplet interactions;
$\phi_{1}$ impacts oscillators 1 and 3 with amplitudes
$4.9\text{\times}{10}^{-3}$ and $4.4\text{\times}{10}^{-3}$, and $\phi_{2}$
with $2.3\text{\times}{10}^{-3}$ and $3.2\text{\times}{10}^{-3}$,
respectively. However, the dynamics of oscillators 2 and 4 are only impacted
by triplet interactions $\phi_{1}$ (with amplitude
$1.33\text{\times}{10}^{-2}$ ) and $\phi_{2}$ ($1.7\text{\times}{10}^{-2}$),
respectively. These triplet interactions describe phase fluctuations over the
long time scale (red curves in Fig. 1 (e)). Therefore, we can conclude that
the phase dynamics of the oscillators coupled in a ring can be described by a
hypernetwork shown in Fig. 1 (f).

The fact that model recovery provides triplets as the best description is
rather puzzling. Also given that the resonant behavior
$\omega_{2}\approx\omega_{4}$ did not appear in the model recovery from data.
This suggests an interplay between the resonant frequencies and the network
topology. The question arises, which resonances/triplet interactions emerge
from a large number of possibilities in a given network, natural frequencies,
and nonlinear coupling?  An outstanding question is what is the origin of
these triplet interactions that were generated by pairwise physical coupling?

### A theory for emergent higher-order interactions

To answer these questions, we develop a theory that captures the important
characteristics of the experiments: nonlinear coupling and triplet resonance
conditions. We consider the networks

$\dot{z}_{k}=f_{k}(z_{k})+\alpha\sum_{\ell=1}^{n}A_{k\ell}h_{k}(z_{k},z_{\ell})$
(6)

where $z_{k}\in\mathbb{C}$ is the state of the $k$th oscillator,
$h_{k}:\mathbb{C}\times\mathbb{C}\rightarrow\mathbb{C}$ is the pairwise
coupling function, $A_{k\ell}$ is the adjacency matrix, and $\alpha>0$ is the
coupling strength. When the isolated system is close to a Hopf bifurcation,
the dynamics is described by
$f_{k}(z_{k})=\gamma_{k}z_{k}+\beta_{k}z_{k}|z_{k}|^{2}$ [28]. The Hopf
bifurcation is a common route to oscillations in nonlinear systems and
describes the appearance of oscillations in applications [2, 5, 6, 8, 3]. Our
proofs are valid for $\gamma_{k}=\lambda+i\omega_{k}$ with small $\lambda$ and
$\omega_{k}$ satisfying resonance conditions. We fix $\beta_{k}=-1$, but this
value is immaterial. We develop a normal form theory to eliminate unnecessary
terms of $h(z_{k},z_{\ell})$ and to expose higher-order ones that predict the
dynamics. To a network of the form of Eq. (S161) we associate non-resonance
conditions that allow us to get rid of the leading interaction terms in
$\alpha$.

Since $h(z_{k},z_{\ell})$ is a linear combination of monomials and the theory
can be applied to each monomial independently, we assume first that
$h(z_{k},z_{\ell})$ is a single monomial of the form

$h(z_{k},z_{\ell})=z_{k}^{d_{1}}\bar{z}_{k}^{d_{2}}z_{\ell}^{d_{3}}\bar{z}_{\ell}^{d_{4}}$
(7)

for non-negative numbers $d_{1},\dots,d_{4}$.  Our major theoretical result is
a formulation of a non-resonance condition given by

$(d_{1}-d_{2}-1)\omega_{k}+(d_{3}-d_{4})\omega_{\ell}\not=0.$ (8)

This condition shows up naturally in our approach, as a monomial Equation (7)
can only be eliminated by a transformation that divides by the left-hand side
of Equation (8). Hence, an interaction term in the coupling function $h$ given
by Equation (7) can only be removed if the non-resonance condition is
satisfied. The non-resonance condition is defined as the union over all non-
resonance conditions of its monomial terms. The network non-resonance
conditions are given by the union over all non-resonance conditions of
$h(z_{k},z_{\ell})$ for which $A_{{k\ell}}\not=0$. Our result is the
following:

In Methods, we show that given Eq. (S161) with
$h:\mathbb{C}\times\mathbb{C}\rightarrow\mathbb{C}$ a smooth map with
vanishing constant terms, under the network non-resonance conditions, there is
a coordinate transformation that eliminates pairwise interaction terms and
reveals the higher-order interactions. The proof consists of two main steps:

i) Existence of a polynomial change of variables. Consider

$u_{k}=z_{k}-\alpha P_{k}$ (9)

for some polynomials $P_{k}$. The goal is to design $P_{k}$ such that in the
variables $u_{k}$ interaction terms linear in $\alpha$ vanish. We obtain
higher-order interactions of order $\alpha^{2}$. For Eq. (S161) we use

$P_{k}(z)=\sum_{\ell=1}^{n}A_{k\ell}\tilde{h}_{k\ell}(z_{k},z_{\ell})\,,$ (10)

where $\tilde{h}_{k\ell}(z,w)$ is the function obtained from $h(z,w)$ by
transforming each monomial according to the following replacement rule:

$z^{d_{1}}\bar{z}^{d_{2}}w^{d_{3}}\bar{w}^{d_{4}}\mapsto\frac{z^{d_{1}}\bar{z}^{d_{2}}w^{d_{3}}\bar{w}^{d_{4}}}{(d_{1}-1)\gamma_{k}+d_{2}\bar{\gamma}_{k}+d_{3}\gamma_{\ell}+d_{4}\bar{\gamma}_{\ell}}\,$
(11)

Note that the imaginary part of the denominator in Equation (11) is precisely
the left-hand side of Equation (8). While bringing the equations to the new
form, we face a major challenge to understand the combinatorial behavior of
the Taylor coefficients during the transformation. We define a bracket on the
space of polynomials to track these coefficients.

ii) Dealing with transformed isolated dynamics. The second major challenge
lies in the fact that another coordinate transformation is needed to eliminate
terms coming from the isolated dynamics $f_{k}$. Indeed, as we eliminate
coupling terms linear in $\alpha$, other terms linear in $\alpha$ appear due
to the isolated dynamics. A remarkable fact is that the same non-resonance
conditions also ensure that the second transformation exists.

Our theorem is applicable to a much broader class of coupling functions and
network formalisms than what is described by Eq. (S161). A rich variety of new
interaction rules can emerge, depending on the specifics of the set-up (see
Supplementary Note 2).

Applying the replacement rule Eq. (11) we obtain

$\displaystyle\dot{u}_{k}=f_{k}(u_{k})-\alpha^{2}\left\\{\sum_{\ell=1}^{n}\sum_{p=1}^{n}\left[A_{k\ell}A_{kp}\prescript{1}{}{G}_{k}^{\ell
p}(u_{k},u_{\ell},u_{p})-A_{k\ell}A_{\ell p}\prescript{2}{}{G}_{k}^{\ell
p}(u_{k},u_{\ell},u_{p})\right]\right\\},$ (12)

up to higher-order terms in $\alpha$ and $u$. In Methods, we discuss the new
coupling functions $\prescript{1}{}{G}_{k}$ and $\prescript{2}{}{G}_{k}$ some
their properties. The coupling is now $\alpha^{2}$ explaining anomalous
synchronization transitions that appears in networks (see Supplementary Note
3).

### Emergent hypernetworks explain experimental data

Similar to the experiments we consider a ring of four oscillators with
coupling function

$h(z,w)=z\bar{w}+z^{2}\bar{w}.$ (13)

Instead of delay, the oscillators are coupled through a conjugate variable
that enables a streamlined theoretical treatment. Close to a Hopf bifurcation,
the delay would have an effect of advancing the oscillations over half a
period. As before, we consider $\omega_{1}-\omega_{2}+\omega_{3}$ and
$\omega_{1}-\omega_{4}+\omega_{3}$ to be close to zero, so, capturing the
triplet resonance in the experiments. We apply our theory to this case to
unravel how higher-order interactions appear in the data.

The coupling function is a combination of $z\bar{w}$ and $z^{2}\bar{w}$,
providing $d_{1}=1$ and $d_{4}=1$ for the first monomial and $d_{1}=2$ and
$d_{4}=1$ for the latter. The resonance condition Eq. (8) is satisfied for
both. Using the replacement rule Eq. (11), we find

$\displaystyle
u_{k}=z_{k}+\alpha\left(\frac{z_{k-1}z_{k}}{\bar{\gamma}_{k-1}}+\frac{z_{k}z_{k+1}}{\bar{\gamma}_{k+1}}+\frac{z_{k-1}^{2}\bar{z}_{k}}{\gamma_{k-1}+\bar{\gamma}_{k}}+\frac{z_{k}^{2}\bar{z}_{k+1}}{\gamma_{k}+\bar{\gamma}_{k+1}}\right)$
(14)

Each node equation contains $16$ interaction terms as in Eq. (12). We discuss
some of these terms for the first node. ${}^{2}G^{23}_{1}$ appears as node $1$
is connected to node $2$ and $2$ to $3$. This interaction is resonant, see
Figure 2 (a). ${}^{2}G^{43}_{1}$ appears because node $1$ is connected to $4$
and node $4$ to $3$. This term is also resonant, see Figure 2 (b).
${}^{1}G^{24}_{1}$ is nonzero and nonresonant. This term appear as $1$ is
directed connected to $2$ and $4$, see Figure 2 (c). Finally, the term
${}^{2}G^{24}_{1}$ is a forbidden, the term would appear from an interaction
of $1$ to $2$ and from $2$ to $4$, however, in the original network the later
interaction is absent, see Figure 2 (d). Remarkably, not all interactions are
relevant when the goal is to describe slow oscillations in the phases.

Figure 2: Emergent higher-order interactions from the original ring network.
Coupling functions appearing in Eq. (12) of node $1$. Colors correspond to
signs in the phase combination with blue standing for positive and orange for
negative. (a) resonant interaction term appearing as ${}^{2}G^{23}_{1}$. (b)
resonant interaction term appearing as ${}^{2}G^{43}_{1}$. Finally, (c) is a
nonresonant term and (d) ${}^{2}G^{24}_{1}$ is a forbidden term (it does not
appear). These new interaction terms can be predicted from the combinatorics
of the original network and coupling function.

Indeed, once we analyse the phases in the new equations, the coupling term
coming from ${}^{2}G_{1}^{23}$ will lead to oscillations with frequency close
to $\omega_{1}-\omega_{2}+\omega_{3}$ while the term coming from
${}^{2}G_{1}^{43}$ leads to a frequency close to
$\omega_{1}-\omega_{4}+\omega_{3}$. This implies that both terms are slowly
varying. In contrast, the term coming from ${}^{2}G_{1}^{24}$ leads to
oscillations with frequency
$\omega_{1}-\omega_{2}+\omega_{4}\approx\omega_{1}$ and is fast oscillating in
comparison to the slow terms with small frequencies. In virtue of the
averaging theory, such fast oscillating terms can be neglected. In fact, only
resonant terms connected by local trees in the original graph will survive
such as the resonant ones involving $\omega_{1}-\omega_{2}+\omega_{3}$ and
$\omega_{1}-\omega_{4}+\omega_{3}$. This yields

$\displaystyle\begin{split}\dot{u}_{1}&=f_{1}(u_{1})-\alpha^{2}\eta_{12}u_{1}^{2}\bar{u}_{2}u_{3}-\alpha^{2}\eta_{14}u_{1}^{2}\bar{u}_{4}u_{3}\\\
\dot{u}_{2}&=f_{2}(u_{2})-\alpha^{2}\zeta_{231}u_{2}^{2}\bar{u}_{1}\bar{u}_{3}\\\
\dot{u}_{3}&=f_{3}(u_{3})-\alpha^{2}\eta_{32}u_{3}^{2}\bar{u}_{2}u_{1}-\alpha^{2}\eta_{34}u_{3}^{2}\bar{u}_{4}u_{1}\\\
\dot{u}_{4}&=f_{4}(u_{4})-\alpha^{2}\zeta_{431}u_{4}^{2}\bar{u}_{1}\bar{u}_{3}\end{split}$
(15)

where $\eta_{pq}=\frac{1}{\gamma_{p}+\bar{\gamma}_{q}}$ and
$\zeta_{pqr}=\frac{2}{\gamma_{p}+\bar{\gamma}_{q}}+\frac{2}{\gamma_{p}+\bar{\gamma}_{r}}+\frac{1}{\bar{\gamma}_{q}}+\frac{1}{\bar{\gamma}_{r}}$.
Writing $u=re^{i\theta}$ we obtain equations for the phases $\theta$. The
averaging theorem gives

$\displaystyle\begin{split}\dot{\theta}_{1}&=\omega_{1}-\alpha^{2}r_{0}^{3}\left[\rho_{12}(\phi_{1})+\rho_{14}(\phi_{2})\right],\\\
\dot{\theta}_{2}&=\omega_{2}-\alpha^{2}r_{0}^{3}\sigma_{231}(\phi_{1})\\\
\dot{\theta}_{3}&=\omega_{3}-\alpha^{2}r_{0}^{3}\left[\rho_{32}(\phi_{1})+\rho_{34}(\phi_{2})\right],\\\
\dot{\theta}_{4}&=\omega_{4}-\alpha^{2}r_{0}^{3}\sigma_{431}(\phi_{2}),\end{split}$
(16)

where the phases $\phi_{1}$ and $\phi_{2}$ are given in Eq. (S165). The
functions $\rho$ and $\sigma$ are provided in the Supplementary Note 4. The
emergent hypernetwork explains the experimental fitting found in Eq. (S28).
These functions represent hyperlinks as shown in Figure 1 (f).

The phase triplets $\phi_{1}$ and $\phi_{2}$ are revealed from phase reduction
in the normal form and they are not obvious from the original Eq. (S161). We
confirm these predictions by direct simulations of Eq. (S161) (Supplementary
Note 5). We present examples for a three-node path in Supplementary Note 6 and
a six-node network in Supplementary Note 7.

### Predicting the slow phase interactions in experiments

In Supplementary Note 3, we show that the experimental recovery of a
hypernetwork is not an artifact. Rather, we prove that imposing sparsity
unavoidably leads to the recovery of the normal form instead. Indeed, as the
recovery allows for a small least square deviation between the data and the
model, the recovery finds the hypernetwork as a simpler description of the
system. So, by measuring the original variables and attempting a model
recovery while imposing sparsity, model recovery learns only the higher-order
interactions. We now use the emergent network prediction for the ring network
with the corresponding resonance conditions as in the experiment to explain
the slow phase dynamics.

From the data we extract the slow phases $\phi_{1}$ and $\phi_{2}$ as shown in
Figure 3 in solid lines. Using our theory, from Eq. (16), we obtain that

$\displaystyle\dot{\phi}_{i}=\Omega_{i}+\sum_{j=1}^{2}a_{ij}\cos\phi_{j}+b_{ij}\sin\phi_{j}$
(17)

where $a$’s and $b$’s are given in terms of the functions $\sigma$ and $\rho$
in Eq. (16) see Supplementary Note 5. We treat $a$’s and $b$’s as fitting
parameters from the vector field in Eq. (17) obtained from first principles,
since the corresponding coupling parameter and amplitudes are unknown. The
resulting solutions agree with the experimental data as seen in Figure 3.  Our
findings are not strictly limited to electrochemical oscillators. As shown in
Supplementary Note 9, we detected the same hypernetworks in nonlinearly
coupled integrate-and-fire neuron models.

Figure 3: Normal form theory explains the experimental results. We show the
time series of the slow phase $\phi_{1}$ and $\phi_{2}$ from experimental data
(solid) and the prediction of the emergent hypernetwork (dashed) capturing
higher-order interactions. The vector field describing the phase interaction
is obtained from first principles. The coefficients of the vector field are
obtained by least-square minimization.

### Emergent hypernetworks among network modules coupled through mean-fields

The requirement of a nonlinear coupling, at first sight, seems to be a
limitation for practical applications. However, here we analyze how
hypernetworks emerge in modular networks with microscopic pairwise coupling
through phase differences.

We consider four subpopulations of $N$ interacting Kuramoto oscillators [13].
Nodes in each subpopulation interact strongly among themselves with coupling
strength $\mu$ and weakly between subgroups with coupling strength $\alpha$,
see Figure 4. As we will show at the macroscopic mean-field level, the
interaction is nonlinear. According to our theory, although the mean-fields
have a pairwise interaction, their model recovery will be in terms of
hypernetworks. We first consider the microscopic description; each oscillator
is described by

$\dot{\psi}_{km}=\omega_{km}+\frac{\mu}{N}\sum_{n=1}^{N}\sin(\psi_{kn}-\psi_{km})+\sum_{\ell=1}^{4}A_{kl}\left(\frac{\alpha}{N}\sum_{n=1}^{N}\sin(\psi_{ln}-\psi_{km})\right)$
(18)

or in terms of mean-fields $\dot{\psi}_{km}=\omega_{km}+{\rm Im}\left(\mu
z_{k}+\alpha\sum A_{kl}z_{l}\right)e^{-i\psi_{km}}$ where

$z_{k}=\frac{1}{N}\sum_{m=1}^{N}e^{i\psi_{km}}$ (19)

is the mean-field of the subpopulation $k$. The frequencies $\omega_{km}$ are
distributed according to a Lorenzian $\rho(\omega,\Omega_{k},\sigma_{k})$
where $\Omega_{k}$ is the mean subpopulation frequency and $\sigma_{k}$ is the
frequency dispersion. Applying the Ott-Antonsen ansatz [15], we obtain the
macroscopic equations describing the mean-fields in the limit
$N\rightarrow\infty$ as

$\dot{z}_{k}=f_{k}(z_{k})+\sum_{l=1}^{4}A_{kl}h(z_{k},z_{l})$ (20)

where $f_{k}$ is the Hopf normal form with constants
$\gamma_{k}=(i\Omega_{k}+\mu-\sigma_{k})$ and $\beta_{k}=-\mu$ and

$h(z_{k},z_{l})=\alpha z_{l}+\alpha\bar{z}_{l}z_{k}^{2},$ (21)

thus, in the macroscopic description the coupling is nonlinear. We interpret
$\alpha$ as a bifurcation parameter and deal with $\alpha z_{l}$ as a
nonlinear term as in bifurcation theory. We consider the ensemble frequencies
to satisfy the resonance conditions $\Omega_{1}+\Omega_{3}\approx 2\Omega_{2}$
and $\Omega_{2}+\Omega_{4}\approx 2\Omega_{1}$. At $\alpha=0$ each
subpopulation will have an order parameter behaving as
$z_{k}(t)=r_{k}e^{i\theta_{k}(t)}$ where
$r_{k}=\sqrt{\frac{\mu-\sigma_{k}}{\mu}}$ and $\dot{\theta}_{k}=\Omega_{k}$.
To obtain the phase model, we bring the network to its normal form and apply
the phase reduction. In Supplementary Note 10, we perform the calculations of
such resonance conditions to obtain the new normal form equations. After
discarding nonresonant terms the phase equations of the mean-fields read as

$\displaystyle\begin{split}\dot{\theta}_{1,3}&=\Omega_{1,3}+F_{1,3}(\varphi_{1})\\\
\dot{\theta}_{2,4}&=\Omega_{2,4}+F_{2,4}(\varphi_{2})\end{split}$ (22)

where $F_{i}$ is a linear combination of sine and cosine.

Next, we fix the ensemble frequencies as
$\Omega_{1}=2,\Omega_{2}=3,\Omega_{3}=4$ and $\Omega_{4}=1$ as well as the
coupling strengths $\mu=0.5$, $\sigma_{k}=0.48$ yielding $r_{k}=0.15$ and
$\alpha=0.1$ for all subpopulations. We numerically integrate the mean-field
equations and obtain the complex fields $z_{1}(t),z_{2}(t),z_{3}(t)$ and
$z_{4}(t)$ which enables us to extract the phase dynamics
$\theta_{1}(t),\theta_{2}(t),\theta_{3}(t)$ and $\theta_{4}(t)$. Performing a
Lasso regression we recover the vector fields of Eq. (22) confirming the
theoretical prediction of higher order interactions, see Supplementary Note
10.

Figure 4: Interacting subpopulations lead to higher order interaction of mean-
fields. a) The original network of coupled subpopulations (with four distinct
colours, namely, red, yellow, blue and orange). Oscillators are interacting by
an internal coupling constant $\mu$ and inter-subpopulations coupling constant
$\alpha$. b) Higher order phase interaction of the mean-fields represented
with the same colors as in a) (red, yellow, blue and orange). Applying our
approach we uncover that the phase interaction between the mean-fields is
described by a hypernetwork. c) The mean-field slow phase variables
$\varphi_{1}$ (green) and $\varphi_{2}$ (purple) were computed from the data
collected from the simulations of mean fields on the associated network. The
dashed curve is the simulation of the vector field of the slow phases
$\varphi_{1,2}$ reconstructed from data using the Lasso method.

As before, we introduce the slow phases

$\displaystyle\begin{split}\varphi_{1}=\theta_{1}-2\theta_{2}+\theta_{3},\\\
\varphi_{2}=\theta_{2}-2\theta_{1}+\theta_{4}.\end{split}$ (23)

The theory predicts the higher order interaction between the slow phases as
$\dot{\varphi}_{k}=\varepsilon_{k}+G_{k}(\varphi_{1},\varphi_{2})$, as shown
in Supplementary Note 10. The fitting the predicted vector field of $\varphi$
to the data is excellent as can be observed in Figure 4 c).

For these four subpopulation on a ring, the condition on the frequencies is
close to the subspace

$V_{res}=\\{\Omega_{1}+\Omega_{3}=2\Omega_{2},\,\Omega_{2}+\Omega_{4}=2\Omega_{1}\\}\,,$
forming a co-dimension $2$ resonance surface. That is, the emergence of
hypernetworks is generic in a two parameter family of frequencies.

## Discussion

We have uncovered a mechanism by which nonlinear pairwise interactions with
triplet resonance conditions result in nontrivial phase dynamics on a
hypernetwork.  Such interactions traditionally were attributed in brain
dynamics to synaptic transmission between two neurons mediated by chemical
messengers from a third neuron (heterosynaptic plasticity) [29]. Our findings
provide an alternative mechanism. On one hand, this finding shows that phase
dynamics can be mediated through ‘virtual’ interactions not physically present
in the system. On the other hand, such a mechanism could be leveraged to
design interactions between remote components not directly connected but
instead having correlations in natural frequencies.

The experimental system with a generic network motif with a ring of four
electrochemical oscillators presented here was an example, where a relatively
simple nonlinear modulation of the coupling induced a hypernetwork driven
phase dynamics. Networks with a ring topology are selected for the experiment
since they are common for many network based complex systems, e.g., in lasers,
biological systems, neuronal dynamics and many disciplines [30, 31]. Such
nonlinear modulation of the coupling can be quite general in gene expressions;
for example, it was used to describe the coupling among circadian cells
through Michaelis-Menten mechanism where coupling from one cell modulated the
maximum gene expression rate in the other [32].

Strikingly, we showed that the coupling resulting in mean-field coupling among
network modules has sufficient nonlinearity to facilitate hypernetwork
interactions.  In particular, event related modulation of spectral responses
of magnetoencephalogram (MEG) recordings (i.e., modulation of frequency-
specific oscillations in the motor network established by a handgrip task)
have shown very strong evidence for nonlinear, between-frequency coupling of
remote brain regions[33]. Our results strongly suggest that in these MEG
recordings, given the appropriate resonances and nonlinearities, hypernetwork
description could facilitate the long-range modulation of frequencies. In
conclusion, the findings open new avenues for hypernetwork based description
and engineering of complex systems with heterogeneous frequencies and
nonlinear interactions.

## Methods

Our results give an algorithmic procedure for obtaining a hypernetwork that
accurately describes the observed behavior of the original system. This
emergent higher order system depends on details of the given network, the
original coupling function and the resonance relations among the phases.

### Normal form calculations

In Supplementary Note 2, we consider ODEs of the general form

$\displaystyle\dot{z}_{k}$
$\displaystyle=\gamma_{k}z_{k}-\beta_{k}z_{k}|z_{k}|^{2}+\alpha
H_{k}(z_{1},\dots,z_{n})\,,\quad k\in\\{1,\dots,n\\}\,,$ (24)

with $z_{k}\in\mathbb{C}$ and $\alpha\in\mathbb{R}$. The numbers
$\beta_{k},\gamma_{k}\in\mathbb{C}$ are assumed non-zero, and we furthermore
write $\gamma_{k}=\lambda+i\omega_{k}$. Here $\lambda\in\mathbb{R}$ is seen as
the bifurcation parameter for a Hopf bifurcation, and we assume the
interaction functions $H_{k}:\mathbb{C}^{n}\rightarrow\mathbb{C}$ to be smooth
(i.e. $C^{\infty}$) for convenience. Moreover, we initially assume each
$H_{k}$ satisfies $H_{k}(0)=0$ and $DH_{k}(0)=0$, though the condition on its
derivative is later dropped.

Our main result shows that the ODE (24) can be put in a normal form that
allows us to predict the phase dynamics of the oscillators. We do this by
using two successive transformations:

$\displaystyle w_{k}$ $\displaystyle=z_{k}-\alpha P_{k}(z);$ (25)
$\displaystyle u_{k}$ $\displaystyle=w_{k}-\alpha Q_{k}(w),\,$ (26)

for some appropriately chosen polynomials $P_{k}$ and $Q_{k}$. The first of
these coordinate transformations is used to remove the term $\alpha H_{k}(z)$
from the Equation (24). This will generate additional terms in $\alpha^{2}$
that may be expressed in the coefficients of $H_{k}$ and $P_{k}$ following
certain combinatorial rules. We manage this combinatorial behavior by
introducing a special bracket $[\bullet||\bullet]$ on the space of
polynomials. In addition to these new interaction terms, the transformation
will also produce terms in $\alpha$ involving $P_{k}$ and
$\beta_{k}z_{k}|z_{k}|^{2}$, which obscure an interpretation of the system as
a (hyper) network. We therefore remove these additional terms using the second
coordinate transformation. A crucial observation here is that the non-
resonance conditions needed for the first transformation are sufficient to
ensure the second. We are able to prove this using the precise bookkeeping
enabled by the aforementioned bracket.

When dealing with the case where $DH_{k}(0)\not=0$, we instead remove only the
non-linear terms in $H_{k}$ using the transformations (25) and (26). This
reveals higher order terms as before. Even though $DH_{k}(0)$ accounts only
for nonresonant terms by assumption, this linear term will nevertheless cause
an overall frequency shift that has to be accounted for. More precisely, if we
denote by $\Omega$ the diagonal matrix with entries the frequencies
$\omega_{1},\dots,\omega_{n}$, then the natural frequencies in the coupled
case will be given by the imaginary part of the eigenvalues of $i\Omega+\alpha
DH(0)$. Here we have set $H=(H_{1},\dots,H_{n})$. These new frequencies can be
approximated by standard eigenvalue perturbation techniques.

### Properties of the coupling functions $\prescript{1}{}{G}_{k}^{\ell p}$ and
$\prescript{2}{}{G}_{k}^{\ell p}$

Applying the transformation of the theorem to Eq. (S161) yields a new system
of the form Eq. (12). In Supplementary Note 2, we show that

$\displaystyle\ \begin{split}\prescript{1}{}{G}_{k}^{\ell
p}(u_{k},u_{\ell},u_{p})&=\frac{\partial\tilde{h}_{k\ell}(u_{k},u_{\ell})}{\partial
u_{k}}h(u_{k},u_{p})+\frac{\partial\tilde{h}_{k\ell}(u_{k},u_{\ell})}{\partial\bar{u}_{k}}\overline{h(u_{k},u_{p})}\,\\\
\prescript{2}{}{G}_{k}^{\ell
p}(u_{k},u_{\ell},u_{p})&=\frac{\partial\tilde{h}_{k\ell}(u_{k},u_{\ell})}{\partial
u_{\ell}}h(u_{\ell},u_{p})+\frac{\partial\tilde{h}_{k\ell}(u_{k},u_{\ell})}{\partial\bar{u}_{\ell}}\overline{h(u_{\ell},u_{p})}\,.\end{split}$
(27)

In Eq. (27) a term of degree $d$ in $h$ and a term of degree $\tilde{d}$ in
$\tilde{h}_{k\ell}$ combine to form a term of degree $d+\tilde{d}-1$ in
$\prescript{1}{}{G}_{k}^{\ell p}$. As both $h$ and $\tilde{h}_{k\ell}$ have
terms of degree 2 and higher, we see that $\prescript{1}{}{G}_{k}^{\ell p}$
only has terms of degree 3 and higher. The same holds true for
$\prescript{2}{}{G}_{k}^{\ell p}$, which means that a classical network
description involving directed edges is no longer possible.

The third order terms are moreover easily found by replacing $h$ and
$\tilde{h}_{k\ell}$ in Eq. (27) by their quadratic terms. Likewise, the fourth
order terms are found by replacing $h$ by its quadratic terms and
$\tilde{h}_{k\ell}$ by its cubic terms and vice versa in Eq. (27). We may also
argue that these higher order terms in $\prescript{1}{}{G}_{k}^{\ell p}$ and
$\prescript{2}{}{G}_{k}^{\ell p}$ are non-vanishing in general. Indeed, the
coefficients in front of these terms are rational functions of $\gamma_{k}$
and the coefficients of $h$. Such functions are either identical to the zero
function (which Eq. (27) excludes) or non-vanishing on an open dense set.

New terms emerge that have an interpretation as higher-order interactions. The
two double sums in Eq. (12) have a combinatorial interpretation. The first
double sum counts all pairs of nodes $(\ell,p)$ that both influenced node $k$
in the original network. The second double sum counts all pairs $(\ell,p)$
where $\ell$ influenced $k$ and $p$ influenced $\ell$ and $p$ need not
influence $k$ directly in the old network, so that new node-dependency is
formed.

### An explicit algorithm for predicting the emergent hypernetwork

We present an algorithm for obtaining an emergent hypernetwork from a given
network system. Its input consists of the adjacency matrix $A$, the function
$h$ and the phases $\omega_{1}$ through $\omega_{n}$, and we assume the
nonresonance conditions of the theorem to hold. The algorithm is as follows:

Algorithm 1 Emergent Hypernetworks

Input: Adjacency matrix $A$, coupling function $h$, frequencies and amplitudes
$\gamma_{i}$’s
Output: Hypernetwork and Coupling functions

1:for each $k\in\mathcal{S}$ do

2: for each $\ell\in\mathcal{S}$ do

3: if $A_{k\ell}\not=0$ then

4: form the polynomials $\tilde{h}_{k\ell}(u_{k},u_{\ell})$ by the replacement
rule
$z^{d_{1}}\bar{z}^{d_{2}}w^{d_{3}}\bar{w}^{d_{4}}\mapsto\frac{z^{d_{1}}\bar{z}^{d_{2}}w^{d_{3}}\bar{w}^{d_{4}}}{(d_{1}-1)\gamma_{k}+d_{2}\bar{\gamma}_{k}+d_{3}\gamma_{\ell}+d_{4}\bar{\gamma}_{\ell}}\,$

5: for each $p\in\mathcal{S}$ do

6: if $A_{k\ell}A_{kp}\not=0$ then

7: Compute ${}^{1}G_{k}^{\ell p}$

8: if $A_{k\ell}A_{\ell p}\not=0$ then

9: Compute ${}^{2}G_{k}^{\ell p}$

10: procedure Resonant terms in the coupling functions $G$

11: for each
$u_{k}^{d_{1}}\bar{u}_{k}^{d_{2}}u_{\ell}^{d_{3}}\bar{u}_{\ell}^{d_{4}}u_{p}^{d_{5}}\bar{u}_{p}^{d_{6}}$
monomial of ${}^{1}G_{k}^{\ell p}$ and ${}^{2}G_{k}^{\ell p}$ do

12: if
$(d_{1}-d_{2}-1)\omega_{k}+(d_{3}-d_{4})\omega_{\ell}+(d_{5}-d_{6})\omega_{p}\not=0$
then

13: discard term

14: procedure Remaining monomials are the couplings of node $k$

## Data Availability

We provide the experimental time-series and the extracted phases of the
oscillations (Fig. 1) at [34].

## Code availability

The source code for reconstructing the functions representing hypernetwork
dynamics from oscillatory networks dynamics is available [35].

## References

  * [1] Watts, M., Tabak, J., Zimliki, C., Sherman, A. & Bertram, R.  Slow variable dominance and phase resetting in phantom bursting.  _J. Theor. Biol._ 276, 218–228 (2011).
  * [2] Kralemann, B. _et al._ In vivo cardiac phase response curve elucidates human respiratory heart rate variability.  _Nat. Commun._ 4, 1–9 (2013).
  * [3] Sebek, M., Tönjes, R. & Kiss, I. Z.  Complex rotating waves and long transients in a ring network of electrochemical oscillators with sparse random cross-connections.  _Phys. Rev. Lett._ 116, 068701 (2016).
  * [4] Bick, C., Sebek, M. & Kiss, I. Z.  Robust weak chimeras in oscillator networks with delayed linear and quadratic interactions.  _Phys. Rev. Lett._ 119, 168301 (2017).
  * [5] Schneidman, E., Berry, M. J., Segev, R. & Bialek, W.  Weak pairwise correlations imply strongly correlated network states in a neural population.  _Nature_ 440, 1007–1012 (2006).
  * [6] Ermentrout, G. B. & Terman, D. H.  _Mathematical foundations of neuroscience_ , vol. 35 (Springer Science & Business Media, 2010).
  * [7] Blasius, B., Huppert, A. & Stone, L.  Complex dynamics and phase synchronization in spatially extended ecological systems.  _Nature_ 399, 354–359 (1999).
  * [8] Matheny, M. H. _et al._ Exotic states in a simple network of nanoelectromechanical oscillators.  _Science_ 363 (2019).
  * [9] Smeal, R. M., Ermentrout, G. B. & White, J. A.  Phase-response curves and synchronized neural networks.  _Philos. Trans. R. Soc. Lond., B, Biol. Sci._ 365, 2407–2422 (2010).
  * [10] Omel’chenko, E. & Wolfrum, M.  Nonuniversal transitions to synchrony in the sakaguchi-kuramoto model.  _Phys. Rev. Lett._ 109, 164101 (2012).
  * [11] Hong, H. & Strogatz, S. H.  Kuramoto Model of Coupled Oscillators with Positive and Negative Coupling Parameters: An Example of Conformist and Contrarian Oscillators.  _Phys. Rev. Lett._ 106, 054102 (2011).
  * [12] Kuramoto, Y.  _Chemical oscillations, waves, and turbulence_ (Courier Corporation, 2003).
  * [13] Stankovski, T., Pereira, T., McClintock, P. V. E. & Stefanovska, A.  Coupling functions: Universal insights into dynamical interaction mechanisms.  _Rev. Mod. Phys._ 89, 045001 (2017).
  * [14] Rodrigues, F. A., Peron, T. K. D., Ji, P. & Kurths, J.  The kuramoto model in complex networks.  _Physics Reports_ 610, 1–98 (2016).
  * [15] Tönjes, R., Fiore, C. E. & Pereira, T.  Coherence resonance in influencer networks.  _Nature Communications_ 12, 1–8 (2021).
  * [16] Giusti, C., Pastalkova, E., Curto, C. & Itskov, V.  Clique topology reveals intrinsic geometric structure in neural correlations.  _Proc. Natl. Acad. Sci. U.S.A._ 112, 13455–13460 (2015).
  * [17] Reimann, M. W. _et al._ Cliques of neurons bound into cavities provide a missing link between structure and function.  _Frontiers in computational neuroscience_ 11, 48 (2017).
  * [18] Bassett, D. S., Zurn, P. & Gold, J. I.  On the nature and use of models in network neuroscience.  _Nat. Rev. Neurosci._ 19, 566–578 (2018).
  * [19] Millán, A. P., Torres, J. J. & Bianconi, G.  Explosive higher-order kuramoto dynamics on simplicial complexes.  _Phys. Rev. Lett._ 124, 218301 (2020).
  * [20] Bick, C., Ashwin, P. & Rodrigues, A.  Chaos in generically coupled phase oscillator networks with nonpairwise interactions.  _Chaos_ 26, 094814 (2016).
  * [21] Kori, H., Kuramoto, Y., Jain, S., Kiss, I. Z. & Hudson, J. L.  Clustering in globally coupled oscillators near a hopf bifurcation: theory and experiments.  _Phys. Rev. E_ 89, 062906 (2014).
  * [22] Giusti, C., Ghrist, R. & Bassett, D. S.  Two’s company, three (or more) is a simplex.  _J. Comput. Neurosci._ 41, 1–14 (2016).
  * [23] Bassett, D. S. & Sporns, O.  Network neuroscience.  _Nat. Neurosci._ 20, 353–364 (2017).
  * [24] Grilli, J., Barabás, G., Michalska-Smith, M. J. & Allesina, S.  Higher-order interactions stabilize dynamics in competitive network models.  _Nature_ 548, 210–213 (2017).
  * [25] Skardal, P. S. & Arenas, A.  Abrupt desynchronization and extensive multistability in globally coupled oscillator simplexes.  _Phys. Rev. Lett._ 122, 248301 (2019).
  * [26] Mulas, R., Kuehn, C. & Jost, J.  Coupled dynamics on hypergraphs: Master stability of steady states and synchronization.  _Phys. Rev. E_ 101, 062313 (2020).
  * [27] Bilal, S. & Ramaswamy, R.  Synchronization and amplitude death in hypernetworks.  _Phys. Rev. E_ 89, 062923 (2014).
  * [28] Shil’nikov, L., Shil’nikov, A., Turaev, D. & Chua, L.  _Methods of qualitative theory in nonlinear dynamics_ , vol. 5 (World Scientific, 2001).
  * [29] Chistiakova, M., Bannon, N. M., Bazhenov, M. & Volgushev, M.  Heterosynaptic plasticity.  _Neuroscientist_ 20, 483–498 (2014).
  * [30] Popovych, O. V., Yanchuk, S. & Tass, P. A.  Delay-and coupling-induced firing patterns in oscillatory neural loops.  _Physical review letters_ 107, 228102 (2011).
  * [31] Takamatsu, A. _et al._ Spatiotemporal symmetry in rings of coupled biological oscillators of physarum plasmodial slime mold.  _Physical Review Letters_ 87, 078102 (2001).
  * [32] Schroder, S., Herzog, E. D. & Kiss, I. Z.  Transcription-Based Oscillator Model for Light-Induced Splitting as Antiphase Circadian Gene Expression in the Suprachiasmatic Nuclei.  _J. Biol. Rhytms_ 27, 79–90 (2012).
  * [33] Chen, C. C. _et al._ Nonlinear Coupling in the Human Motor System.  _J. Neurosci._ 30, 8393–8399 (2010).
  * [34] Nijholt, E., Ocampo-Espindola, J. L., Eroglu, D., Kiss, I. Z. & Pereira, T.  Emergent hypernetworks in weakly coupled oscillators (this paper). GitHub repository which includes the data measured from the experiments. (2021).  URL https://github.com/jloespindola/Hypernetwork_data.
  * [35] Nijholt, E., Ocampo-Espindola, J. L., Eroglu, D., Kiss, I. Z. & Pereira, T.  Emergent hypernetworks in weakly coupled oscillators (this paper). Zenodo repository: code for computing the functions representing hypernetwork dynamics from the networks dynamics. (2021).  URL https://doi.org/10.5281/zenodo.5749164.
  * [36] Pikovsky, A., Kurths, J., Rosenblum, M. & Kurths, J.  _Synchronization: a universal concept in nonlinear sciences_ , vol. 12 (Cambridge university press, 2003).
  * [37] Eroglu, D., Lamb, J. S. W. & Pereira, T.  Synchronisation of chaos and its applications.  _Contemporary Physics_ 58, 207–243 (2017).
  * [38] Novaes, M., Santos, E. R. d. & Pereira, T.  Recovering sparse networks: Basis adaptation and stability under extensions.  _arXiv preprint arXiv:2104.00796_ (2021).
  * [39] Kori, H., Kiss, I. Z., Jain, S. & Hudson, J. L.  Partial synchronization of relaxation oscillators with repulsive coupling in autocatalytic integrate-and-fire model and electrochemical experiments.  _Chaos_ 28, 045111 (2018).  URL https://doi.org/10.1063/1.5022497.
  * [40] Kaptanoglu, A. A. _et al._ Pysindy: A comprehensive python package for robust sparse system identification.  _Journal of Open Source Software_ 7, 3994 (2022).  URL https://doi.org/10.21105/joss.03994.

## Acknowledgments

We thank Sajjad Bakrani, Zachary G. Nicolaou, Marcel Novaes, Edmilson Roque,
Robert Ronge and Jeroen Lamb for enlightening discussions. TP was supported in
part by FAPESP Cemeai Grant No. 2013/07375-0 and is a Newton Advanced Fellow
of the Royal Society NAF$\backslash$R1$\backslash$180236\. TP and EN were
partially supported by Serrapilheira Institute (Grant No. Serra-1709-16124).
DE was supported by TUBITAK Grant No. 118C236 and the BAGEP Award of the
Science Academy. JLO-E acknowledges financial support from CONACYT. IZK
acknowledges support from National Science Foundation (grant CHE-1900011).

## Author Contributions Statement

EN and TP designed the overall study and formulated the theory. JLO-E and IZK
designed and performed the experiments. DE implemented the numerical
simulations and analyses. All authors contributed to the writing of the
manuscript. All authors reviewed and approved the final manuscript.

## Competing Interests Statement

The authors declare no competing interests.

List of supplementary materials

Supplementary Text

Supplementary Information to
Emergent hypernetworks in weakly coupled oscillators

Eddie Nijholt1, Jorge Luis Ocampo-Espindola2, Deniz Eroglu3,

István Z. Kiss2, Tiago Pereira1,4∗

1Instituto de Ciências Matemáticas e Computação, Universidade de São Paulo,
São Carlos, Brazil

2Department of Chemistry, Saint Louis University, St. Louis, USA

3 Faculty of Engineering and Natural Sciences, Kadir Has University, Istanbul,
Turkey

4 Department of Mathematics, Imperial College London, SW7 2AZ, London, United
Kingdom

###### Supplementary Notes

  1. 1 Experimental setup and methods
    1. 1.1 Experimental setup
    2. 1.2 Behavior without coupling
    3. 1.3 Phase dynamics
  2. 2 Proof of emergent higher-order networks
    1. 2.1 Preliminaries and results
    2. 2.2 A coordinate transformation
    3. 2.3 The first reduction
    4. 2.4 The second reduction
  3. 3 Anomalous synchronization on a 4-node ring
    1. 3.1 Explanation of anomalous synchronization via normal form theory
  4. 4 Phase reduction for $h=(z^{2}+z)\bar{w}$ and resonance $\omega_{1}-\omega_{2,4}+\omega_{3}=0$
  5. 5 Slow phase dynamics for $h=(z^{2}+z)\bar{w}$ and resonance $\omega_{1}-\omega_{2,4}+\omega_{3}=0$
    1. 5.1 Emergent hypernetwork predicts data behaviour
  6. 6 Model recovery of a 3-path with coupling $h=(z^{2}+z)\bar{w}$ and resonance $\omega_{1}-\omega_{2}+\omega_{3}=0$
    1. 6.1 Emergent network explanation
  7. 7 6 nodes network examples with $h(z,w)=z\bar{w}$
  8. 8 Model recovery and normal form representation
  9. 9 Emergent hypernetworks in an integrate-and-fire model
  10. 10 Mean field interaction
    1. 10.1 Frequency shifts
    2. 10.2 Model Recovery
    3. 10.3 Normal Form Calculations

Supplementary Note

##  1 Experimental setup and methods

In this section, we describe the details about the experimental setup, the
dynamical behavior of the oscillators without coupling, the of the phase model
using LASSO.

###  1.1 Experimental setup

The experiments were carried out in a standard three-electrode electrochemical
cell.

The cell consists of a nickel-array working electrode (W), a Pt-coated Ti rod
as a counter electrode (C), and a Hg/Hg2SO4 sat. K2SO4 as a reference
electrode (R). The electrolyte was a 3.0 M sulfuric acid solution at a
constant temperature of 10 ∘C. The electrode array consisted of four 1-mm
diameter nickel wires embedded in epoxy with a spacing of 3 mm. A multichannel
potentiostat (Gill-IK64, ACM Instruments) interfaced with a real-time LabVIEW
controller measured the potential drop $[E_{k}(t)$, with respect to the
reference electrode] and set the circuit potential ($V_{0,k}$) of the working
electrodes individually at a rate of 200 Hz. The electrode potentials are
corrected for offset $o_{1}$=0.92 V, $o_{2}$=0.98 V, $o_{3}$= 0.91 V, and
$o_{4}$=0.97 V.

###  1.2 Behavior without coupling

The offset circuit potential to each oscillator was established 20 mV above
the Hopf bifurcation ($V_{0,1}$=1850 mV, $V_{0,2}$=1100 mV, $V_{0,3}$=1660 mV,
$V_{0,4}$=1103 mV). The natural frequencies [Supplementary Fig. 5 (a)] were
adjusted to have values of $\omega_{1}$=0.152 Hz, $\omega_{2}$=0.385 Hz,
$\omega_{3}$=0.237 Hz and $\omega_{4}$=0.384 Hz with a set of resistors and
capacitors $R_{ind,1}$=12.0 kohm, $C_{ind}$=440 $\mu$F, $R_{ind,2}$=1.00 kohm,
$R_{ind,3}$=12.0 kohm, $R_{ind,4}$=1.00 kohm. Without coupling, we observed
that the slow oscillators (1 and 3) have about twice the amplitude than the
fast oscillators (2 and 4). Supplementary Fig. 5 (c) shows the electrode
potential time series of each oscillator.

###  1.3 Phase dynamics

Phase definition.We used the peak-finding approach [36] to extract the phase
of each oscillator and then linear interpolation between peaks from the
experimental electrode potential time series. When there is no coupling, the
pairwise phase difference shows a linear growth [Supplementary Fig.. 5 (b)]
and the triplet phase differences, $\phi_{j}$, $j$=1, 2, do not show phase
slip behavior [Supplementary Fig. 5 (d)].

Fitting of phase dynamics As described in the main text, the impact of triplet
interactions on the dynamics can be extracted with a LASSO fit to the
$\dot{\theta}_{k}$ values according to

$\dot{\theta}_{k}=\hat{\omega}_{k}(t)+\sum_{j=1}^{2}C_{j}^{k}\sin(\phi_{j})+D_{j}^{k}\cos(\phi_{j})$
(S28)

where
$\hat{\omega}_{k}(t)=\hat{\omega}^{0}_{k}+\hat{\omega}^{1}_{k}t+\hat{\omega}^{2}_{k}t^{2}$
is the fitted, slowly drifting (up to quadratic variation in time) natural
frequency, and $C_{j}^{k}$ and $D_{j}^{k}$ are the amplitudes of the sin and
cos phase coupling functions corresponding to the appropriate triplet phase
differences.

For the fit, the instantaneous frequency, $\dot{\theta}_{k}$ was obtained with
the numerical derivative of the phase of each oscillator from the experimental
times series. The $\dot{\theta}_{k}$ was filtered by a first order Savitzky-
Golay filter for 45 s. Using $\phi_{j}$, we fitted the $\dot{\theta}_{k}$ with
LASSO method. In LASSO, the regularization parameter determines how many
parameters in the fitted model should be set to zero. We used a regularization
parameter so that the mean square error is 20% higher than the best fit (no
regularization). The fitted parameters are:

Table 1: Recovered coefficients from Eq. S28. Coefficients | Oscillator number
---|---
| 1 | 2 | 3 | 4
$\hat{\omega}_{k}^{0}$ |  0.953 |  2.368 |  1.467 |  2.383
$\hat{\omega}_{k}^{1}$ |  $2.76\times 10^{-5}$ | 0 | $-4.61\times 10^{-5}$ | 0
$\hat{\omega}_{k}^{2}$ | $-6.40\times 10^{-8}$ | $-6.63\times 10^{-8}$ | $-1.32\times 10^{-7}$ | 0
$C_{1}^{k}$ |  $4.89\times 10^{-3}$ |  $5.04\times 10^{-3}$ | $\,\,\,\,3.15\times 10^{-3}$ | 0
$D_{1}^{k}$ |  $9.27\times 10^{-5}$ |  $1.23\times 10^{-2}$ | $-3.10\times 10^{-3}$ | 0
$C_{2}^{k}$ | $-1.49\times 10^{-3}$ | 0 |  $3.16\times 10^{-3}$ |  $4.63\times 10^{-3}$
$D_{2}^{k}$ | $-1.73\times 10^{-3}$ | 0 | $-7.37\times 10^{-4}$ |  $1.64\times 10^{-2}$
Figure 5: Dynamical behavior without coupling in the experiments. a) Natural
frequency of each oscillators without coupling $\omega_{1}$=0.152 Hz,
$\omega_{2}$=0.385 Hz, $\omega_{3}$=0.237 Hz and $\omega_{4}$=0.384 Hz. The
blue, orange, yellow and green dots represent the elements 1 to 4
respectively. b) Time series of the phase difference with respect oscillator
one. The blue line: $\theta_{1}$-$\theta_{1}=0$, orange line:
$\theta_{2}$-$\theta_{1}$, yellow line: $\theta_{3}$-$\theta_{1}$ and green
line: $\theta_{4}$-$\theta_{1}$. c) Electrode potential time series. Blue,
orange, yellow, and green line corresponds to oscillator one to four
respectively. d) Time series of the slow phases, $\phi_{1}$ (blue) and
$\phi_{2}$ (green) without coupling.

##  2 Proof of emergent higher-order networks

We consider ODEs of the general form

$\displaystyle\dot{z}_{k}$
$\displaystyle=\gamma_{k}z_{k}-\beta_{k}z_{k}|z_{k}|^{2}+\alpha
H_{k}(z_{1},\dots,z_{n})\,,$ (S29)

for $k\in\\{1,\dots,n\\}$. Here, each $z_{k}$ takes values in $\mathbb{C}$ and
$\alpha\in\mathbb{R}$ denotes the coupling parameter of the interaction. We
moreover have $\beta_{k},\gamma_{k}\in\mathbb{C}$ non-zero, and write
$\gamma_{k}=\lambda+i\omega_{k}$ for all $k\in\\{1,\dots,n\\}$. Note that
$\lambda\in\mathbb{R}$ may be seen as the bifurcation parameter for a Hopf
bifurcation, which might in particular vanish. Each interaction function
$H_{k}:\mathbb{C}^{n}\rightarrow\mathbb{C}$ is assumed smooth (i.e.
$C^{\infty}$) for convenience, and satisfies $H_{k}(0)=0$ and $DH_{k}(0)=0$.

We will show that the ODE (S29) can be put in a particular normal form that
allows us to predict the dynamics of the phases of the oscillators. Our
technique for doing so involves two successive coordinate transformations:

$\displaystyle w_{k}$ $\displaystyle=z_{k}-\alpha P_{k}(z)$ (S30)
$\displaystyle u_{k}$ $\displaystyle=w_{k}-\alpha Q_{k}(w)\,$

for some appropriately chosen polynomials $P_{k}$ and $Q_{k}$. The first of
these transformations is used to remove the term $\alpha H_{k}(z)$ from the
ODE (S29). This will generate additional terms in $\alpha^{2}$ that may be
expressed in the coefficients of $H_{k}$ and $P_{k}$ following certain
combinatorial rules. In order to describe this combinatorial behavior, we
first introduce a useful bracket $[\bullet||\bullet]$ on the space of
polynomials, see Definition 2.6. The first coordinate transformation will also
produce terms in $\alpha$ involving $P_{k}$ and $\beta_{k}z_{k}|z_{k}|^{2}$.
Again our bracket allows for a precise description of these new terms, which
we then remove using the second coordinate transformation. The precise
bookkeeping enabled by the bracket will be crucial in determining what non-
resonance conditions are needed for the second transformation. In fact, it
will turn out that the non-resonance conditions needed for the first
transformation are sufficient to ensure the second.

We first present the main result, Theorem 2.4, in Subsection 2.1. There we
also develop the necessary definitions, notation and machinery needed for the
proof, which is then presented in the remaining subsections.

###  2.1 Preliminaries and results

In order to analyse the ODE (S29), it will be useful to write

$H_{k}(z)={H}^{d}_{k}(z)+\mathcal{O}(|z|^{d+1})\,,$

where ${H}^{d}_{k}(z)$ is a polynomial denoting the terms up to degree $d$ in
the Taylor expansion of $H_{k}(z)$ around the origin. We will mostly work with
the value $d=5$. Note that ${H}^{d}_{k}(z)$ is therefore a polynomial in both
the variables $z_{1},\dots,z_{n}$ and their complex conjugates
$\overline{z}_{1},\dots,\overline{z}_{n}$, with complex coefficients. In
general, whenever we talk about a polynomial we will always mean a complex
polynomial in its given complex variables and their complex conjugates. It
will also be useful to write $H(z)=(H_{k}(z))$ for the vector valued function
that captures all interaction functions $H_{k}(z)$ as its components, and
similarly set $H^{d}(z)=(H^{d}_{k}(z))$.

As is often the case with normal form calculations, we will need to assume
some conditions on the $\omega_{k}$ (or more precisely the $\gamma_{k}$).
These will depend on the coefficients of ${H}^{5}_{k}(z)$. To this end, we
define:

######  Definition 2.1.

Let

$\displaystyle R(z)=cz_{1}^{s_{1}}\dots
z_{n}^{s_{n}}\overline{z}_{1}^{t_{1}}\dots\overline{z}_{n}^{t_{n}}\,$ (S31)

be a monomial term in ${H}^{d}_{k}(z)$, where $c$ is a complex number and
$s_{1},\dots,s_{n}$, $t_{1},\dots,t_{n}$ are non-negative integers. The $k$th
_non-resonance condition_ of $R(z)$ is the condition

$\displaystyle s_{1}\omega_{1}+\dots+s_{n}\omega_{n}-t_{1}\omega_{1}-\dots-
t_{n}\omega_{n}-\omega_{k}\not=0\,.$ (S32)

Note that, as $\omega_{\ell}$ denotes the imaginary part of $\gamma_{\ell}$
for all $\ell\in\\{1,\dots,n\\}$, the $k$th non-resonance condition guarantees
in particular that:

$\displaystyle
s_{1}\gamma_{1}+\dots+s_{n}\gamma_{n}+t_{1}\overline{\gamma}_{1}+\dots+t_{n}\overline{\gamma}_{n}-\gamma_{k}\not=0\,,$
(S33)

which will play a role in much of the proofs and constructions below. In fact,
varying $\lambda$ and allowing the particular case $\lambda=0$, we see that
equations (S32) and (S33) are equivalent in general. Next, the $k$th _non-
resonance condition_ of a polynomial is defined as the union of the $k$th non-
resonance conditions of all of its monomial terms. Finally, the _non-resonance
condition_ of a polynomial map ${H}^{d}(z)=({H}^{d}_{k}(z))$ is the union over
all $k\in\\{1,\dots,n\\}$ of the $k$th non-resonance conditions of
${H}^{d}_{k}(z)$. $\triangle$

######  Example 2.2.

Suppose the interaction functions are given by the polynomials

$H_{k}(z)=\sum_{\ell=1}^{n}c_{k,\ell}z_{k}\overline{z}_{\ell}\,,$ (S34)

for some (possibly weighted) connection matrix
$c=(c_{k,\ell})\in\mathbb{C}^{n\times n}$. It follows that the $k$th non-
resonance condition of $H_{k}(z)$ is given by

$\omega_{k}-{\omega}_{\ell}-\omega_{k}=-{\omega}_{\ell}\not=0\quad\text{ for
all }{\ell}\text{ such that }c_{k,\ell}\not=0\,.$ (S35)

Hence, we see that the non-resonance condition of $H(z)=(H_{k}(z))$ is
satisfied if we simply have $\omega_{\ell}\not=0$ for all nodes
${\ell}\in\\{1,\dots,n\\}$. $\triangle$

######  Example 2.3.

Suppose the interaction functions are given by the polynomials

$H_{k}(z)=\sum_{\ell=1}^{n}c_{k,\ell}(z_{k}\overline{z}_{\ell}+z_{k}^{2}\overline{z}_{\ell})\,,$
(S36)

for some connection matrix $c=(c_{k,\ell})$. The $k$th non-resonance condition
of $H_{k}(z)$ is now given by

$\displaystyle\omega_{k}-\omega_{\ell}-\omega_{k}$
$\displaystyle=-\omega_{\ell}\not=0\quad\text{ and }$ (S37) $\displaystyle
2\omega_{k}-\omega_{\ell}-\omega_{k}$
$\displaystyle=\omega_{k}-\omega_{\ell}\not=0$

for all $\ell$ such that $c_{k,l}\not=0$. If we assume for convenience that
$c$ encodes a symmetric, connected graph, then the non-resonance condition of
$H(z)$ is satisfied if

$\displaystyle\omega_{\ell}$ $\displaystyle\not=0\text{ for all nodes
}\ell\in\\{1,\dots,n\\}\text{ and }$ (S38)
$\displaystyle\omega_{p}-\omega_{q}$ $\displaystyle\not=0\text{ for all edges
}e=[p,q]\text{ between nodes }p\text{ and }q\,.$

$\triangle$

We are now ready to formulate the main theorem. It tells us that, under the
relevant non-resonance conditions, we may transform the ODE (S29) into a
system with leading interaction terms involving only three-way “hyper-
interactions” and with coupling constant $\alpha^{2}$. See Proposition 2.9 for
an exact description of the new leading interaction terms in Theorem 2.4.

######  Theorem 2.4.

Let ${H}^{5}_{k}$ denote the fifth order Taylor expansion of the $k$th
interaction function $H_{k}:\mathbb{C}^{n}\rightarrow\mathbb{C}$. Assume the
non-resonance conditions for ${H}^{5}=({H}^{5}_{k})$ to hold. Then the ODE
(S29) is locally conjugate to

$\displaystyle\dot{u}_{k}$
$\displaystyle=\gamma_{k}u_{k}-\beta_{k}u_{k}|u_{k}|^{2}-\alpha^{2}G_{k}(u)$
(S39)
$\displaystyle+\mathcal{O}(|\alpha||u|^{6}+|\alpha|^{2}|u|^{5}+|\alpha|^{3}|u|^{4})\,,$

with $u_{k}\in\mathbb{C}$ and for some complex polynomials $G_{k}$ with only
terms of degree $3$ and higher. See Proposition 2.9 for an exact description
of the $G_{k}$.

Note that Equation (S39) gives a precise description of $\dot{u}_{k}$ up to
sixth order in $u$ and $\alpha$.

######  Remark 2.5.

Theorem 2.4 tells us that, under the relevant non-resonance conditions, the
$\alpha$-terms can be ‘pushed back’ to order $|\alpha||u|^{6}$. Technically
speaking, this means we can ignore these remainder terms only when $\alpha
u^{6}<<\alpha^{2}u^{4}$, i.e. when $u^{2}<<\alpha$. However, it will be clear
from the proof in subsections 2.2 and 2.4 that if the non-resonance conditions
hold for terms in $H_{k}$ beyond ${H}_{k}^{5}$, then we may arrange for a
remainder in $\alpha$ with higher order terms in $u$. That is, we then get the
remainder
$\mathcal{O}(|\alpha||u|^{N}+|\alpha|^{2}|u|^{5}+|\alpha|^{3}|u|^{4})$ for a
corresponding value of $N>6$. As the interaction functions $H_{k}$ are taken
to be polynomials in our examples, with the non-resonance conditions holding
for all terms, we in fact expect a remainder of the form
$\mathcal{O}(|\alpha||u|^{N}+|\alpha|^{2}|u|^{5}+|\alpha|^{3}|u|^{4})$ for
arbitrarily high value of $N$. Hence, we may neglect all terms that are first
order in $\alpha$, and obtain a new coupled system with coupling constant
$\alpha^{2}$. $\triangle$

The new interaction functions $G_{k}$ can be obtained from $H$ by a
combinatorial construction on the Taylor coefficients. This is best described
using a bracket on polynomials that we define below. We will furthermore make
extensive use of this bracket throughout the proof of Theorem 2.4.

######  Definition 2.6.

Let $R(z)$ be a complex polynomial and let $S(z)=(S_{1}(z),\dots,S_{n}(z))$ be
an $n$-tuple (i.e., a vector) of complex polynomials
$S_{1}(z),\dots,S_{n}(z)$. We let $[R||S](z)$ be the complex polynomial
obtained by (formally) assuming each variable $z_{j}$ is time-dependent (i.e.,
$z_{j}=z_{j}(t)$) and satisfies $\dot{z}_{j}=S_{j}(z)$, after which we
differentiate $R(z)$ with respect to $t$. That is, we set

$\displaystyle[R||S](z):=\frac{d}{dt}R(z)\left|{\begin{array}[]{l}\dot{z}_{j}=S_{j}(w)\\\
\forall\,j\in\\\ \\{1,\dots,n\\}\end{array}}\right.\,.$ (S43)

The reason we choose this notation, instead of one involving the Jacobian of
$R$, is to avoid confusion with the complex conjugate variables
$\overline{z}_{j}$. Because each term $\dot{\overline{z}}_{j}$ is substituted
by $\overline{S_{j}(w)}$, the expression $[R||S]$ is in general not complex
linear in $S$.

In the special case where $S(z)=(\gamma_{1}z_{1},\dots,\gamma_{n}z_{n})$ (with
$\gamma_{k}$ as in Equation (S29)), we set

$\displaystyle\Gamma
R(z):=[R||S](z)=\frac{d}{dt}R(z)\left|{\begin{array}[]{l}\dot{z}_{j}=\gamma_{j}z_{j}\\\
\forall\,j\in\\\ \\{1,\dots,n\\}\end{array}}\right.\,.$ (S47)

$\triangle$

######  Example 2.7.

It is not hard to see that the polynomial $\Gamma R(z)$ is obtained by
replacing every monomial

$cz_{1}^{s_{1}}\dots
z_{n}^{s_{n}}\overline{z}_{1}^{t_{1}}\dots\overline{z}_{n}^{t_{n}}\,,\quad
c\in\mathbb{C}$

in $R(z)$ by

$c(s_{1}\gamma_{1}+\dots+s_{n}\gamma_{n}+t_{1}\overline{\gamma}_{1}+\dots+t_{n}\overline{\gamma}_{n})z_{1}^{s_{1}}\dots
z_{n}^{s_{n}}\overline{z}_{1}^{t_{1}}\dots\overline{z}_{n}^{t_{n}}\,.$

Consider for instance the monomial
$R(z)=R(z_{1},z_{2},\overline{z}_{1},\overline{z}_{2})=z_{1}^{2}\overline{z}_{2}$.
We have

$\displaystyle\frac{d}{dt}R(z)=2z_{1}\dot{z}_{1}\overline{z}_{2}+z_{1}^{2}\dot{\overline{z}}_{2}\,.$
(S48)

Hence, we indeed find

$\displaystyle\Gamma
R(z)=\frac{d}{dt}R(z)\left|{\begin{array}[]{l}\dot{z}_{1}=\gamma_{1}z_{1}\\\
\dot{z}_{2}=\gamma_{2}z_{2}\end{array}}\right.$
$\displaystyle=2z_{1}(\gamma_{1}z_{1})\overline{z}_{2}+z_{1}^{2}(\overline{\gamma_{2}z_{2}})$
(S51)
$\displaystyle=(2\gamma_{1}+\overline{\gamma}_{2})z_{1}^{2}\overline{z}_{2}=(2\gamma_{1}+\overline{\gamma}_{2})R(z)\,.$

$\triangle$

The term $G_{k}$ in Theorem 2.4 will be given as the bracket
$[\bullet||\bullet]$ between ${H}^{5}$ and a polynomial obtained by slightly
modifying ${H}^{5}_{k}$. Hence, intuitively, $G_{k}$ should be thought of as
$[H_{k}||H]$. More precisely, we define:

######  Definition 2.8.

Let $P$ be a polynomial for which its $k$th non-resonance conditions are met.
The ($k$th) _modified polynomial_ $\widehat{P}_{k}$ is obtained from $P$ by
replacing each monomial

$cz_{1}^{s_{1}}\dots
z_{n}^{s_{n}}\overline{z}_{1}^{t_{1}}\dots\overline{z}_{n}^{t_{n}}\,,\quad
c\in\mathbb{C}$

in $P$ by

$\frac{cz_{1}^{s_{1}}\dots
z_{n}^{s_{n}}\overline{z}_{1}^{t_{1}}\dots\overline{z}_{n}^{t_{n}}}{s_{1}\gamma_{1}+\dots+s_{n}\gamma_{n}+t_{1}\overline{\gamma}_{1}+\dots+t_{n}\overline{\gamma}_{n}-\gamma_{k}}\,.$

In the special case where $P={H}^{d}_{k}$ for some $d\leq 5$, we will simply
write $\widehat{H}^{d}_{k}:=\widehat{(H^{d}_{k})}_{k}$ to denote the
corresponding modified polynomial. $\triangle$

######  Proposition 2.9.

In Theorem 2.4 the terms $G_{k}$ are given by

$G_{k}=[\widehat{H}^{3}_{k}||{H}^{3}]\,.$ (S52)

It can be shown that $[\widehat{H}^{3}_{k}||{H}^{3}]$ indeed only has terms of
degree $3$ and higher, using the assumption that each $H_{k}$ (and therefore
each $\widehat{H}^{d}_{k}$) only has terms of degree 2 and higher. See Remark
2.15 below. Note that we only care about the third and fourth order terms of
$G_{k}$, as the rest are absorbed in the remainder terms of Equation (S39). It
will be clear from Remark 2.15 that these lowest order terms do not change if
we instead define

$G_{k}=[\widehat{H}^{d}_{k}||{H}^{d}]\,,$ (S53)

for $d=4$ or $d=5$. For this reason we will often simply write

$G_{k}=[\widehat{H}_{k}||{H}]\,.$ (S54)

######  Remark 2.10.

Let $I_{k}\subset\\{1,\dots,n\\}$ denote the _input set_ of a node
$k\in\\{1,\dots,n\\}$. That is, $I_{k}$ denotes those nodes that influence
$k$, or more precisely those nodes $\ell$ for which

$\frac{\partial H_{k}(z)}{\partial z_{\ell}}\not=0\text{ or }\frac{\partial
H_{k}(z)}{\partial\overline{z}_{\ell}}\not=0\,.$

Note that $I_{k}$ might not contain $k$ itself. It follows that in general
$G_{k}(u)$ depends on variables for nodes in the set

$I^{2}_{k}:=\left(\bigcup_{{\ell}\in I_{k}}I_{\ell}\right)\cup I_{k}\,.$ (S55)

This is because $G_{k}=[\widehat{H}^{3}_{k}||{H}^{3}]$ is formed by replacing
a term $z_{\ell}$ (or $\overline{z}_{\ell}$) in $\widehat{H}^{3}_{k}$ by
${H}^{3}_{\ell}$ (or $\overline{{H}^{3}_{\ell}}$), and this is done for each
$\ell\in I_{k}$. We have also used here that $\widehat{H}^{3}_{k}$ likewise
only depends on the variables corresponding to nodes in $I_{k}$, or possibly a
strict subset thereof.

In a similar way one sees that the third order terms of $G_{k}$ (that is, its
leading order terms) are given by ‘triplet terms’
$\tilde{u}_{r}\tilde{u}_{s}\tilde{u}_{t}$, where we have $r\in I_{k}$ and
$s,t\in I_{\ell}$ for some $\ell\in I_{k}$. (Here each $\tilde{u}_{p}$ may
independently denote $u_{p}$ or its complex conjugate $\overline{u}_{p}$). See
Figure 6 for a schematic depiction of these emergent triplet terms. Of course
the specifics of $H_{k}(z)$ might put additional constraints on the terms that
can show up in $G_{k}(u)$. $\triangle$

kr$\ell$ts$I_{k}$$I_{\ell}$ Figure 6: Schematic depiction of the hidden
‘triplet terms’ that show up in $G_{k}$.

######  Example 2.11.

As in Example 2.2, let us make the particular choice for the interaction
functions

$H_{k}(z)={H}_{k}^{3}(z)=\sum_{\ell=1}^{n}c_{k,\ell}z_{k}\overline{z}_{\ell}\,.$
(S56)

Here $c=(c_{k,\ell})$ is the connectivity matrix of the network. We will
assume that $c_{k,\ell}\in\\{0,1\\}$, or more generally
$c_{k,\ell}\in\mathbb{R}$, if the network is weighted. We find

$\widehat{H}_{k}(z)=\widehat{H}^{3}_{k}(z)=\sum_{\ell=1}^{n}\frac{c_{k,\ell}}{\overline{\gamma}_{\ell}}z_{k}\overline{z}_{\ell}\,.$
(S57)

A direct calculation now shows that

$\displaystyle G_{k}(u)=[\widehat{H}_{k}||H](u)$
$\displaystyle=\sum_{\ell=1}^{n}\frac{c_{k,\ell}}{\overline{\gamma}_{\ell}}(H_{k}(u)\overline{u}_{\ell}+u_{k}\overline{H_{\ell}(u)})$
(S58)
$\displaystyle=\sum_{\ell=1}^{n}\frac{c_{k,\ell}}{\overline{\gamma}_{\ell}}\left(\sum_{p=1}^{n}c_{k,p}u_{k}\overline{u}_{p}\overline{u}_{\ell}+u_{k}\sum_{p=1}^{n}c_{\ell,p}\overline{u}_{\ell}u_{p}\right)$
$\displaystyle=\sum_{\ell=1}^{n}\sum_{p=1}^{n}\frac{c_{k,\ell}c_{k,p}}{\overline{\gamma}_{\ell}}u_{k}\overline{u}_{\ell}\overline{u}_{p}+\sum_{\ell=1}^{n}\sum_{p=1}^{n}\frac{c_{k,\ell}c_{\ell,p}}{\overline{\gamma}_{\ell}}u_{k}\overline{u}_{\ell}u_{p}\,.$

The first of these two terms counts all trees in the network that are formed
by the node $k$ and two nodes that influence node $k$. The second term counts
all trees formed by the nodes $k$, $\ell$ and $p$, where $k$ depends on $\ell$
and $\ell$ depends on $p$. $\triangle$

We will gather some more facts about the bracket $[\bullet||\bullet]$. One
important observation is given by:

######  Lemma 2.12.

The map $(R,S)\mapsto[R||S]$ is complex linear in $R$ and real linear in $S$.
In other words, given polynomials $R,R^{\prime}$ and complex numbers
$\lambda,\mu$, we have

$[\lambda R+\mu R^{\prime}||S]=\lambda[R||S]+\mu[R^{\prime}||S]\,.$ (S59)

Given tuples $S,S^{\prime}$ and real numbers $a,b$ we have

$[R||aS+bS^{\prime}]=a[R||S]+b[R||S^{\prime}]\,.$ (S60)

###### Proof.

Complex linearity in $R$ is clear from the definition:

$\displaystyle[R||S](z):=\frac{d}{dt}R(z)\left|{\begin{array}[]{l}\dot{z}_{j}=S_{j}(w)\\\
\forall\,j\in\\\ \\{1,\dots,n\\}\end{array}}\right.\,.$ (S64)

Because of this, it suffices to show real linearity in $S$ when $R$ is given
by a monomial of the form

$R(z)=z_{i_{1}}z_{i_{2}}\dots
z_{i_{m}}\overline{z}_{j_{1}}\overline{z}_{j_{2}}\dots\overline{z}_{j_{l}},$

for some (not necessarily distinct)
$i_{1},\dots,i_{m},j_{1},\dots,j_{l}\in\\{1,\dots,n\\}$. We get

$\displaystyle[R||S](z)=\sum_{s=1}^{m}\frac{R(z)}{z_{i_{s}}}S_{i_{s}}(z)+\sum_{r=1}^{l}\frac{R(z)}{\overline{z}_{i_{r}}}\overline{S_{i_{r}}(z)}\,,$
(S65)

from which real linearity in $S$ follows readily. ∎

We will also make extensive use of the following definition:

######  Definition 2.13.

Given a complex monomial

$Q(z)=cz_{1}^{s_{1}}\dots
z_{n}^{s_{n}}\overline{z}_{1}^{t_{1}}\dots\overline{z}_{n}^{t_{n}}\,,$

where $c\in\mathbb{C}$ and with $s_{1},\dots,s_{n},t_{1},\dots,t_{n}$ non-
negative integers, we define the degree of $Q(z)$ as the number
$s_{1}+\dots+s_{n}+t_{1}+\dots+t_{n}$. The degree of a polynomial $P(z)$ is
then defined as the maximum of the degrees of all the monomial terms of
$P(z)$, as is common. Similarly, we define the _lower degree_ of a polynomial
$P(z)$ as the minimum of the degrees of all of its monomial terms. $\triangle$

It turns out our bracket has a predictable effect on degrees:

######  Lemma 2.14.

Let $R$ be a polynomial and $S=(S_{1},\dots,S_{n})$ a vector of polynomials.
If $R$ has degree $p>0$ and each of the polynomial components of $S$ has
degree at most $p^{\prime}$, then $[R||S]$ has degree $p+p^{\prime}-1$ or
lower. If $R$ has lower degree $q>0$ and each of the polynomial components of
$S$ has lower degree at least $q^{\prime}$, then $[R||S]$ has lower degree
$q+q^{\prime}-1$ or higher.

###### Proof.

By linearity of the bracket (see Lemma 2.12), it suffices to show that the
degree of $[R||S]$ is $d+d^{\prime}-1$ if $R$ is a monomial of degree $d$ and
the components of $S$ are all monomials of degree $d^{\prime}$. (By
convention, we treat the zero-polynomial as a polynomial of any degree.) As in
the proof of Lemma 2.12, we write

$R(z)=z_{i_{1}}z_{i_{2}}\dots
z_{i_{m}}\overline{z}_{j_{1}}\overline{z}_{j_{2}}\dots\overline{z}_{j_{l}},$

for some (not necessarily distinct)
$i_{1},\dots,i_{m},j_{1},\dots,j_{l}\in\\{1,\dots,n\\}$. It follows that
$m+l=d$. As in the previous proof we find

$\displaystyle[R||S](z)=\sum_{s=1}^{m}\frac{R(z)}{z_{i_{s}}}S_{i_{s}}(z)+\sum_{r=1}^{l}\frac{R(z)}{\overline{z}_{i_{r}}}\overline{S_{i_{r}}(z)}\,,$
(S66)

from which it follows readily that the degree of $[R||S]$ is indeed
$d+d^{\prime}-1$. This completes the proof. ∎

######  Remark 2.15.

Lemmas 2.12 and 2.14 give us an easy way of finding the third and fourth order
terms of $G_{k}=[\widehat{H}_{k}||{H}]$ (and higher terms if necessary). It
follows that the third order terms of $G_{k}$ are given by
$[\widehat{H}^{2}_{k}||{H}^{2}]$, where we use that $H$ (and therefore each
$\widehat{H}^{2}_{k}$) has no constant and linear terms. Likewise, the fourth
order terms of $G_{k}$ are given by the bracket between the second order terms
of $\widehat{H}_{k}$ and the third order terms of $H$, plus the bracket
between the third order terms of $\widehat{H}_{k}$ and the second order terms
of $H$. Note also that $G_{k}$ need only be defined up to third and fourth
order in Theorem 2.4, as higher order terms of $\alpha^{2}G_{k}(u)$ can be
absorbed in the remainder $\mathcal{O}(|\alpha|^{2}|u|^{5})$ of Equation
(S39). For this reason the terms of degree 4 and higher in $\widehat{H}_{k}$
and ${H}$ play no role in the construction of (the relevant terms of)
$G_{k}=[\widehat{H}_{k}||{H}]$. $\triangle$

Next, we consider $G_{k}$ in the situation where $H_{k}$ describes a coupled
cell system as in the examples of the main text.

######  Remark 2.16.

Suppose each $H_{k}$ is of the special form

$H_{k}(z)=\sum_{\ell=1}^{n}c_{k,\ell}h(z_{k},z_{\ell})\,,$ (S67)

where $z=(z_{1},\dots,z_{n})$. Here $(c_{k,\ell})$ is a real adjacency matrix
and $h:\mathbb{C}^{2}\rightarrow\mathbb{C}$ has vanishing constant and linear
terms. It follows that we may write

$\widehat{H}_{k}(z)=\sum_{\ell=1}^{n}c_{k,\ell}\widehat{h}_{k,\ell}(z_{k},z_{\ell})\,,$
(S68)

where $\widehat{h}_{k,\ell}(z_{k},z_{\ell})$ is obtained from
$h(z_{k},z_{\ell})$ by applying a monomial substitution to its terms in
precisely the same way $\widehat{H}_{k}$ is obtained from $H_{k}$. By
linearity of the bracket $[\bullet||\bullet]$ in the first slot, we get

$\displaystyle
G_{k}=[\widehat{H}_{k}||{H}]=[\sum_{\ell=1}^{n}c_{k,\ell}\widehat{h}_{k,\ell}(u_{k},u_{\ell})||H]=\sum_{\ell=1}^{n}c_{k,\ell}[\widehat{h}_{k,\ell}(u_{k},u_{\ell})||H]\,.$
(S69)

Moreover, we find

$\displaystyle[\widehat{h}_{k,\ell}(u_{k},u_{\ell})||H](u)$
$\displaystyle=\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial
u_{k}}H_{k}+\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial\overline{u}_{k}}\overline{H}_{k}+\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial
u_{\ell}}H_{\ell}+\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial\overline{u}_{\ell}}\overline{H}_{\ell}$
(S70)
$\displaystyle=\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial
u_{k}}\left(\sum_{p=1}^{n}c_{k,p}h(u_{k},u_{p})\right)+\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial\overline{u}_{k}}\left(\sum_{p=1}^{n}c_{k,p}\overline{h(u_{k},u_{p})}\right)$
$\displaystyle+\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial
u_{\ell}}\left(\sum_{p=1}^{n}c_{\ell,p}h(u_{\ell},u_{p})\right)+\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial\overline{u}_{\ell}}\left(\sum_{p=1}^{n}c_{\ell,p}\overline{h(u_{\ell},u_{p})}\right)$
$\displaystyle=\sum_{p=1}^{n}c_{k,p}\left(\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial
u_{k}}h(u_{k},u_{p})+\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial\overline{u}_{k}}\overline{h(u_{k},u_{p})}\right)$
$\displaystyle+\sum_{p=1}^{n}c_{\ell,p}\left(\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial
u_{\ell}}h(u_{\ell},u_{p})+\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial\overline{u}_{\ell}}\overline{h(u_{\ell},u_{p})}\right)$
$\displaystyle=\sum_{p=1}^{n}c_{k,p}\prescript{1}{}{G}_{k}^{\ell
p}(u_{k},u_{\ell},u_{p})+\sum_{p=1}^{n}c_{\ell,p}\prescript{2}{}{G}_{k}^{\ell
p}(u_{k},u_{\ell},u_{p})\,,$

where we have set

$\displaystyle\prescript{1}{}{G}_{k}^{\ell p}(u_{k},u_{\ell},u_{p})$
$\displaystyle:=\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial
u_{k}}h(u_{k},u_{p})+\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial\overline{u}_{k}}\overline{h(u_{k},u_{p})}\quad\text{
and }$ (S71) $\displaystyle\prescript{2}{}{G}_{k}^{\ell
p}(u_{k},u_{\ell},u_{p})$
$\displaystyle:=\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial
u_{\ell}}h(u_{\ell},u_{p})+\frac{\partial\widehat{h}_{k,\ell}(u_{k},u_{\ell})}{\partial\overline{u}_{\ell}}\overline{h(u_{\ell},u_{p})}\,.$
(S72)

Combining equations (S69) through (S72), we obtain

$\displaystyle
G_{k}(u)=\sum_{\ell=1}^{n}\sum_{p=1}^{n}c_{k,\ell}c_{k,p}\prescript{1}{}{G}_{k}^{\ell
p}(u_{k},u_{\ell},u_{p})+\sum_{\ell=1}^{n}\sum_{p=1}^{n}c_{k,\ell}c_{\ell,p}\prescript{2}{}{G}_{k}^{\ell
p}(u_{k},u_{\ell},u_{p})\,.$ (S73)

We may interpret Equation (S73) as representing a new interaction structure,
one where the interaction is now encoded through certain trees in the graph
instead of links. In this regard, the emergent interaction function (S73)
looks a lot like our original response function (S67), but counting such trees
instead of links. The only way in which Equation (S73) does not generalize
Equation (S67) perfectly is by the fact that $\prescript{1}{}{G}_{k}^{\ell p}$
and $\prescript{1}{}{G}_{k}^{\ell p}$ have indices $k,\ell$ and $p$ (whereas
${h}$ does not). However, we see from equations (S71) and (S72) that there is
no dependence on $p$; this index is only there for notational purposes.
Moreover, the dependence on $k$ and $\ell$ is only through a rescaling of the
monomials. Hence, we find an emergent interaction that is in very good
agreement with a generalization of our original interaction to tree
interaction. What is more, the trees that Equation (S73) counts are easily
identified in the original graph. See Example 2.11, which describes a special
case of interaction through (S67), and the corresponding Figure 3 of the main
manuscript. $\triangle$

######  Example 2.17.

We return to Example 2.3, where the interaction functions are given by

$H_{k}(z)=\sum_{\ell=1}^{n}c_{k,\ell}(z_{k}\overline{z}_{\ell}+z_{k}^{2}\overline{z}_{\ell})\,.$
(S74)

This is of the form (S67) as discussed in Remark 2.16, with $h$ given by

$h(z_{k},z_{\ell})=(z_{k}+z_{k}^{2})\overline{z}_{\ell}\,.$ (S75)

Following the notation of Remark 2.16, we see that

$\widehat{h}_{k,\ell}(z_{k},z_{\ell})=\frac{z_{k}\overline{z}_{\ell}}{\overline{\gamma}_{\ell}}+\frac{z_{k}^{2}\overline{z}_{\ell}}{\gamma_{k}+\overline{\gamma}_{\ell}}.$
(S76)

We therefore find

$\displaystyle\prescript{1}{}{G}_{k}^{\ell p}(u_{k},u_{\ell},u_{p})$
$\displaystyle=\left(\frac{\overline{z}_{\ell}}{\overline{\gamma}_{\ell}}+\frac{2z_{k}\overline{z}_{\ell}}{\gamma_{k}+\overline{\gamma}_{\ell}}\right)(z_{k}+z_{k}^{2})\overline{z}_{p}=\frac{(z_{k}+z_{k}^{2})\overline{z}_{\ell}\overline{z}_{p}}{\overline{\gamma}_{\ell}}+\frac{2(z^{2}_{k}+z_{k}^{3})\overline{z}_{\ell}\overline{z}_{p}}{\gamma_{k}+\overline{\gamma}_{\ell}}\,,$
(S77) $\displaystyle\prescript{2}{}{G}_{k}^{\ell p}(u_{k},u_{\ell},u_{p})$
$\displaystyle=\left(\frac{z_{k}}{\overline{\gamma}_{\ell}}+\frac{z_{k}^{2}}{\gamma_{k}+\overline{\gamma}_{\ell}}\right)(\overline{z}_{\ell}+\overline{z}_{\ell}^{2}){z}_{p}=\frac{z_{k}(\overline{z}_{\ell}+\overline{z}_{\ell}^{2}){z}_{p}}{\overline{\gamma}_{\ell}}+\frac{z_{k}^{2}(\overline{z}_{\ell}+\overline{z}_{\ell}^{2}){z}_{p}}{\gamma_{k}+\overline{\gamma}_{\ell}}\,.$

As we may ignore terms of degree 5 and higher, we may also set

$\displaystyle\prescript{1}{}{G}_{k}^{\ell p}(u_{k},u_{\ell},u_{p})$
$\displaystyle=\frac{(z_{k}+z_{k}^{2})\overline{z}_{\ell}\overline{z}_{p}}{\overline{\gamma}_{\ell}}+\frac{2z^{2}_{k}\overline{z}_{\ell}\overline{z}_{p}}{\gamma_{k}+\overline{\gamma}_{\ell}}\,,$
(S78) $\displaystyle\prescript{2}{}{G}_{k}^{\ell p}(u_{k},u_{\ell},u_{p})$
$\displaystyle=\frac{z_{k}(\overline{z}_{\ell}+\overline{z}_{\ell}^{2}){z}_{p}}{\overline{\gamma}_{\ell}}+\frac{z_{k}^{2}\overline{z}_{\ell}{z}_{p}}{\gamma_{k}+\overline{\gamma}_{\ell}}\,,$

which describe the new interaction through Equation (S73). $\triangle$

###  2.2 A coordinate transformation

In this subsection and the next ones we prove Theorem 2.4 and the accompanying
Proposition 2.9. Recall that we want to transform

$\displaystyle\dot{z}_{k}$
$\displaystyle=\gamma_{k}z_{k}-\beta_{k}z_{k}|z_{k}|^{2}+\alpha H_{k}(z)\,,$
(S79)

into an ODE where the leading interaction terms are of order $\alpha^{2}$.
Recall as well that $H_{k}^{5}(z)$ denotes the Taylor expansion of $H_{k}(z)$
up to fifth order. In particular, we may write

$\displaystyle H_{k}(z)=H^{5}_{k}(z)+\mathcal{O}(|z|^{6})\,.$ (S80)

It follows that $H^{5}_{k}(z)$ is a complex polynomial of order $5$ in the
variables $z_{1},\dots,z_{n}$ and $\overline{z}_{1},\dots,\overline{z}_{n}$.
We write $z=(z_{1},\dots,z_{n})$, and similarly for other variables, and
assume implicitly that any function of $z$ may also depend on its complex
conjugate $\overline{z}=(\overline{z}_{1},\dots,\overline{z}_{n})$.

We start by rewriting the ODE (S79) using the transformation

$w_{k}=z_{k}-\alpha\widehat{H}_{k}^{5}(z)=z_{k}-\alpha P_{k}(z)\,,$ (S81)

where we have set $P_{k}:=\widehat{H}_{k}^{5}$ for convenience. Note that each
$P_{k}:\mathbb{C}^{n}\rightarrow\mathbb{C}$ is a complex polynomial of lower
degree $2$ (see Definition 2.13). It follows that Expression (S81) describes
an invertible transformation around $z=0$. The following lemma deals with its
inverse.

######  Lemma 2.18.

Suppose the variables $w=(w_{1},\dots,w_{n})$ may be expressed in
$z=(z_{1},\dots,z_{n})$ and $\alpha$ by

$w_{k}=z_{k}-\alpha P_{k}(z)\,,$ (S82)

for some polynomials $P_{k}$ of lower degree $d\geq 2$. Then $z$ can be
expressed in $w$ and $\alpha$ by the formal expression

$z_{k}=w_{k}+\alpha
P_{k}(w)+\alpha^{2}R_{k,2}(w)+\alpha^{3}R_{k,3}(w)+\dots\,.$ (S83)

Here the $R_{k,t}(w)$ are polynomials with lower degree $(d-1)t+1$ or higher.

###### Proof.

We write

$z_{k}=R_{k,0}(w)+\alpha R_{k,1}(w)+\alpha^{2}R_{k,2}(w)+\dots\,,$ (S84)

for some functions $R_{k,t}(w):\mathbb{C}^{n}\rightarrow\mathbb{C}$. To
determine these functions, we substitute the $z$ variables in Equation (S82)
by Expression (S84). We obtain

$\displaystyle w_{k}$ $\displaystyle=z_{k}-\alpha P_{k}(z_{1},\dots,z_{n})$
(S85) $\displaystyle=[R_{k,0}(w)+\alpha R_{k,1}(w)+\dots]-\alpha
P_{k}([R_{1,0}(w)+\alpha R_{1,1}(w)+\dots],\dots,[R_{n,0}(w)+\alpha
R_{n,1}(w)+\dots])\,.$

Comparing constant terms in $\alpha$ (i.e. $\alpha^{0}$), Expression (S85)
gives us

$w_{k}=R_{k,0}(w)\,.$ (S86)

This simplifies Equation (S85) to

$\displaystyle w_{k}$ $\displaystyle=[w_{k}+\alpha R_{k,1}(w)+\dots]-\alpha
P_{k}([w_{1}+\alpha R_{1,1}(w)+\dots],\dots,[w_{n}+\alpha
R_{n,1}(w)+\dots])\,.$ (S87)

Comparing $\alpha$-terms now yields

$0=R_{k,1}(w)-P_{k}(w)\,,$ (S88)

so that

$R_{k,1}(w)=P_{k}(w)\,.$ (S89)

It remains to show that the higher order terms are indeed polynomials of the
required lower degree. We will show this by induction on $t$. Note that
$R_{k,0}(w)=w_{k}$ has (lower) degree $(d-1)0+1=1$. Likewise,
$R_{k,1}(w)=P_{k}(w)$ is of lower degree $(d-1)1+1=d$. We therefore fix an
integer $T>1$ and assume that the function $R_{k,t}(w)$ is a complex
polynomial of lower degree $(d-1)t+1$ or higher for all $t<T$ and
$k\in\\{1,\dots,n\\}$. The $\alpha^{T}$ terms in Equation (S87) are given by

$\displaystyle 0$ $\displaystyle=R_{k,T}(w)-[\alpha^{T-1}]P_{k}([w_{1}+\alpha
R_{1,1}(w)+\dots],\dots,[w_{n}+\alpha R_{n,1}(w)+\dots])\,.$ (S90)

Here $[\alpha^{T-1}]F(\alpha)$ denotes the $\alpha^{T-1}$ term in the
expansion of a function $F$ in $\alpha$. As $P_{k}$ is a polynomial of lower
degree $d$, the $\alpha^{T-1}$ term in

$P_{k}([w_{1}+\alpha R_{1,1}(w)+\dots],\dots,[w_{n}+\alpha R_{n,1}(w)+\dots])$

must be a finite sum of scalar multiples of expressions of the form

$\tilde{R}_{i_{1},t_{1}}(w)\tilde{R}_{i_{2},t_{2}}(w)\dots\tilde{R}_{i_{s},t_{s}}(w)\,,$

for $s\geq d$ and for some $i_{1},\dots,i_{s}\in\\{1,\dots,n\\}$ and
$t_{1},\dots,t_{s}\in\mathbb{Z}_{\geq 0}$ satisfying $t_{1}+\dots+t_{s}=T-1$.
Each term $\tilde{R}_{i_{j},t_{j}}(w)$ may furthermore independently denote
${R}_{i_{j},t_{j}}(w)$ or its complex conjugate
$\overline{R_{i_{j},t_{j}}(w)}$. As we have $t_{1}+\dots+t_{s}=T-1$, it in
particular holds that $t_{1},\dots,t_{s}\leq T-1$. By the induction
hypothesis, we therefore know that each of the terms $\tilde{R}_{i_{j},t_{j}}$
is a polynomial of lower degree $(d-1)t_{j}+1$ or higher. This means the
expression

$\tilde{R}_{i_{1},t_{1}}(w)\tilde{R}_{i_{2},t_{2}}(w)\dots\tilde{R}_{i_{s},t_{s}}(w)\,$

is a polynomial of lower degree $D$ satisfying

$\displaystyle D$
$\displaystyle\geq[(d-1)t_{1}+1]+[(d-1)t_{2}+1]+\dots+[(d-1)t_{s}+1]$
$\displaystyle=(d-1)(t_{1}+t_{2}+\dots+t_{s})+s$
$\displaystyle=(d-1)(T-1)+s\geq(d-1)(T-1)+d$ $\displaystyle=(d-1)T+1\,.$

It follows from Equation (S90) that $R_{k,T}(w)$ is indeed a polynomial of
lower degree $(d-1)T+1$ or higher for all $k\in\\{1,\dots,n\\}$. This proves
the lemma by induction. ∎

Setting $d=2$, it follows that the inverse of Equation (S81) is given by

$z_{k}=w_{k}+\alpha P_{k}(w)+\mathcal{O}(|\alpha|^{2}|w|^{3})\,.$ (S91)

At some point later on, we will need to know Expression (S91) up to higher
order terms. To this end, we will show how our bracket $[\bullet||\bullet]$
from Definition 2.6 shows up when performing coordinate transformations.

######  Lemma 2.19.

Let $R(z)$ be a complex polynomial and suppose we may express the
$z$-variables in some new $w$-variables by

$z_{k}=w_{k}+\alpha S_{k}(w)+\mathcal{O}(|\alpha|^{2})\,,\quad
k\in\\{1,\dots,n\\}\,.$ (S92)

Here each $S_{k}$ is a complex polynomial and we have $\alpha\in\mathbb{R}$.
Then $R(z)$ is given in the $w$-variables by

$R(z)=R(w)+\alpha[R||S](w)+\mathcal{O}(|\alpha|^{2})\,,$ (S93)

where we have set $S=(S_{1},\dots,S_{n})$.

###### Proof.

We write $R(z)$ as

$R(z)=T_{0}(w)+\alpha T_{1}(w)+\mathcal{O}(|\alpha|^{2})\,,$ (S94)

where $T_{0}(w)$ and $T_{1}(w)$ are to be determined. Assume first that $R(z)$
is given by

$R(z)=z_{i_{1}}z_{i_{2}}\dots
z_{i_{m}}\overline{z}_{j_{1}}\overline{z}_{j_{2}}\dots\overline{z}_{j_{l}},$

for some (not necessarily distinct)
$i_{1},\dots,i_{m},j_{1},\dots,j_{l}\in\\{1,\dots,n\\}$. We get

$\displaystyle R(z)$ $\displaystyle=z_{i_{1}}\dots
z_{i_{m}}\overline{z}_{j_{1}}\dots\overline{z}_{j_{l}}$ (S95)
$\displaystyle=(w_{i_{1}}+\alpha S_{i_{1}}(w))\dots(w_{i_{m}}+\alpha
S_{i_{m}}(w))\overline{(w_{j_{1}}+\alpha
S_{j_{1}}(w))}\dots\overline{(w_{j_{l}}+\alpha
S_{j_{l}}(w))}+\mathcal{O}(|\alpha|^{2})$ $\displaystyle=w_{i_{1}}\dots
w_{i_{m}}\overline{w}_{j_{1}}\dots\overline{w}_{j_{l}}+\alpha\left(\sum_{s=1}^{m}\frac{R(w)}{w_{i_{s}}}S_{i_{s}}(w)+\sum_{r=1}^{l}\frac{R(w)}{\overline{w}_{i_{r}}}\overline{S_{i_{r}}(w)}\right)+\mathcal{O}(|\alpha|^{2})$
$\displaystyle=R(w)+\alpha[R||S](w)+\mathcal{O}(|\alpha|^{2})\,,$

where in the last line we have used Expression (S65) from the proof of Lemma
2.12. As $T_{0}$ and $T_{1}$ are determined linearly by $R$, we may conclude
from Lemma 2.12 that $T_{0}(w)=R(w)$ and $T_{1}(w)=[R||S](w)$ for general
polynomials $R$. This completes the proof. ∎

######  Example 2.20.

Suppose we are given the polynomial
$Q(z)=Q(z_{1},z_{2},\overline{z}_{1},\overline{z}_{2})=z_{1}^{2}+z_{1}\overline{z}_{2}$.
Equation (S92) gives

$\displaystyle Q(z)$ $\displaystyle=(w_{1}+\alpha
S_{1}(w)+\mathcal{O}(|\alpha|^{2}))^{2}+(w_{1}+\alpha
S_{1}(w)+\mathcal{O}(|\alpha|^{2})\overline{(w_{2}+\alpha
S_{2}(w)+\mathcal{O}(|\alpha|^{2}))}$ (S96) $\displaystyle=(w_{1}+\alpha
S_{1}(w))^{2}+(w_{1}+\alpha
S_{1}(w))(\overline{w}_{2}+\alpha\overline{S_{2}(w)})+\mathcal{O}(|\alpha|^{2})$
$\displaystyle=w_{1}^{2}+w_{1}\overline{w}_{2}+\alpha(2w_{1}S_{1}(w)+w_{1}\overline{S_{2}(w)}+S_{1}(w)\overline{w}_{2})+\mathcal{O}(|\alpha|^{2})$
$\displaystyle=Q(w)+\alpha\frac{d}{dt}Q(w)\left|{\begin{array}[]{l}\dot{w}_{1}=S_{1}(w)\\\
\dot{w}_{2}=S_{2}(w)\end{array}}\right.+\mathcal{O}(|\alpha|^{2})$ (S99)
$\displaystyle=Q(w)+\alpha[Q||S](w)+\mathcal{O}(|\alpha|^{2})\,,$

which is in accordance with Lemma 2.19. $\triangle$

Returning to the transformation (S81) with inverse Equation (S91), we may in
fact conclude the following:

######  Lemma 2.21.

Suppose we have a coordinate transformation of the form

$w_{k}=z_{k}-\alpha P_{k}(z)\,,$ (S100)

where each $P_{k}:\mathbb{C}^{n}\rightarrow\mathbb{C}$ is a complex polynomial
of lower degree $2$ or higher. The inverse transformation is given by

$z_{k}=w_{k}+\alpha
P_{k}(w)+\alpha^{2}[P_{k}||P](w)+\mathcal{O}(|\alpha|^{3}|w|^{4})\,,$ (S101)

where $P=(P_{1},\dots,P_{n})$.

###### Proof.

It follows from Lemma 2.18 that we may write

$z_{k}=w_{k}+\alpha
P_{k}(w)+\alpha^{2}R_{k,2}(w)+\mathcal{O}(|\alpha|^{3}|w|^{4})\,,$ (S102)

for some function $R_{k,2}(w)$. Hence, we only have to show that
$R_{k,2}(w)=[P_{k}||P](w)$. To this end, we rewrite Expression (S100) as

$z_{k}=w_{k}+\alpha P_{k}(z)\,.$ (S103)

Next, we use Equation (S102) to write

$z_{k}=w_{k}+\alpha P_{k}(w)+\mathcal{O}(|\alpha|^{2})\,.$ (S104)

Applying Lemma 2.19 to the term $P_{k}(z)$ and the transformation (S104)
yields

$\displaystyle
P_{k}(z)=P_{k}(w)+\alpha[P_{k}||P](w)+\mathcal{O}(|\alpha|^{2})\,.$ (S105)

Combined with Equation (S103), we obtain

$\displaystyle z_{k}$
$\displaystyle=w_{k}+\alpha(P_{k}(w)+\alpha[P_{k}||P](w)+\mathcal{O}(|\alpha|^{2}))$
(S106) $\displaystyle=w_{k}+\alpha
P_{k}(w)+\alpha^{2}[P_{k}||P](w)+\mathcal{O}(|\alpha|^{3})\,.$

Comparing the two expressions (S102) and (S106) for $z_{k}$, we see that
indeed

$z_{k}=w_{k}+\alpha
P_{k}(w)+\alpha^{2}[P_{k}||P](w)+\mathcal{O}(|\alpha|^{3}|w|^{4})\,.$ (S107)

This proves the lemma. ∎

Our next step is to differentiate Equation (S81) with respect to time. This
gives us

$\displaystyle\dot{w}_{k}=\dot{z}_{k}-\alpha\frac{d}{dt}P_{k}(z)\,.$ (S108)

We will first focus on the term

$\frac{d}{dt}P_{k}(z)\,,$ (S109)

and then deal with the term $\dot{z}_{k}$.

#### The term $\partial_{t}P_{k}$

We first focus on the term (S109). We start by rewriting Equation (S79) as

$\displaystyle\dot{z}_{k}$
$\displaystyle=\gamma_{k}z_{k}-\beta_{k}z_{k}|z_{k}|^{2}+\alpha H_{k}(z)$
(S110)
$\displaystyle=\gamma_{k}z_{k}-\beta_{k}z_{k}|z_{k}|^{2}+\alpha{H}^{5}_{k}(z)+\mathcal{O}(|\alpha||z|^{6})\,,$

where we recall that ${H}^{5}_{k}(z)$ denotes the Taylor expansion of
$H_{k}(z)$ up to fifth order. From Equation (S110) and Lemma 2.12 we get

$\displaystyle\frac{d}{dt}P_{k}(z)$
$\displaystyle=\frac{d}{dt}P_{k}(z)\left|{\begin{array}[]{l}\dot{z}_{j}=\gamma_{j}z_{j}-\beta_{j}z_{j}|z_{j}|^{2}+\alpha{H}_{j}(z)\\\
\forall\,j\in\\{1,\dots,n\\}\\\ \end{array}}\right.$ (S113)
$\displaystyle=\frac{d}{dt}P_{k}(z)\left|{\begin{array}[]{l}\dot{z}_{j}=\gamma_{j}z_{j}-\beta_{j}z_{j}|z_{j}|^{2}+\alpha{H}^{5}_{j}(z)\\\
\forall\,j\in\\{1,\dots,n\\}\\\
\end{array}}\right.+\quad\mathcal{O}(|\alpha||z|^{7})$ (S116)
$\displaystyle=[P_{k}||(\dots,\gamma_{j}z_{j}-\beta_{j}z_{j}|z_{j}|^{2}+\alpha{H}^{5}_{j}(z),\dots)](z)+\mathcal{O}(|\alpha||z|^{7})$
$\displaystyle=[P_{k}||(\dots,\gamma_{j}z_{j},\dots)](z)-[P_{k}||(\dots,\beta_{j}z_{j}|z_{j}|^{2},\dots)](z)$
$\displaystyle+\alpha[P_{k}||(\dots,{H}^{5}_{j}(z),\dots)](z)+\mathcal{O}(|\alpha||z|^{7})$
$\displaystyle=\Gamma
P_{k}(z)-L^{1}_{k}(z)+\alpha[P_{k}||{H}^{5}](z)+\mathcal{O}(|\alpha||z|^{7})\,,$

where we have set

$L^{1}_{k}(z):=[P_{k}||(\dots,\beta_{j}z_{j}|z_{j}|^{2},\dots)](z)$

and

${H}^{5}:=(\dots,{H}^{5}_{j}(z),\dots)\,.$

We have moreover used that $P_{k}(z)$ has lower degree at least $2$ to arrive
at the remainder term $\mathcal{O}(|\alpha||z|^{7})$, and we refer to
Definition 2.6 for the meaning of the term $\Gamma P_{k}(z)$. Note that
$L^{1}_{k}(z)$ is a polynomial of lower degree at least $4$, whereas
$[P_{k}||{H}^{5}](z)$ has lower degree $3$ or higher.

Next, we return to Equation (S91), which we recall states

$z_{k}=w_{k}+\alpha
P_{k}(w)+\mathcal{O}(|\alpha|^{2}|w|^{3})=w_{k}+\mathcal{O}(|\alpha||w|^{2})\,.$
(S117)

We obtain

$\displaystyle L^{1}_{k}(z)$
$\displaystyle=L^{1}_{k}(w)+\mathcal{O}(|\alpha||w|^{5})$ (S118)
$\displaystyle[P_{k}||{H}^{5}](z)$
$\displaystyle=[P_{k}||{H}^{5}](w)+\mathcal{O}(|\alpha||w|^{4})\,.$ (S119)

From Lemma 2.19 we furthermore get

$\displaystyle\Gamma P_{k}(z)$ $\displaystyle=\Gamma P_{k}(w)+\alpha[\Gamma
P_{k}||P](w)+\mathcal{O}(|\alpha|^{2}|w|^{4})\,,$ (S120)

where we have set $P=(P_{1},\dots,P_{n})$. The remainder term in Equation
(S120) follows from the lower degrees of $\Gamma P_{k}(z)$ and $P_{k}(z)$, and
the remainder in Equation (S117).

Note that $[\Gamma P_{k}||P](w)$ is a polynomial of lower degree $3$ or
higher. Combining equations (S113), (S118), (S119) and (S120), we arrive at:

######  Lemma 2.22.

The term

$\frac{d}{dt}P_{k}(z)$

may be expressed in the new $w$ coordinates by

$\displaystyle\frac{d}{dt}P_{k}(z)$ $\displaystyle=\Gamma
P_{k}(z)-L^{1}_{k}(z)+\alpha[P_{k}||{H}^{5}](z)+\mathcal{O}(|\alpha||z|^{7})$
(S121) $\displaystyle=\Gamma
P_{k}(w)-L^{1}_{k}(w)+\alpha[P_{k}||{H}^{5}](w)+\alpha[\Gamma
P_{k}||P](w)+\mathcal{O}(|\alpha||w|^{5}+|\alpha|^{2}|w|^{4})\,.$

#### The term $\dot{z}_{k}$

Next, we focus on the term $\dot{z}_{k}$. Again we write.

$\displaystyle\dot{z}_{k}$
$\displaystyle=\gamma_{k}z_{k}-\beta_{k}z_{k}|z_{k}|^{2}+\alpha H_{k}(z)$
(S122) $\displaystyle=\gamma_{k}z_{k}-\beta_{k}z_{k}|z_{k}|^{2}+\alpha
H_{k}^{5}(z)+\mathcal{O}(|\alpha||z|^{6})\,,$

where ${H}^{5}_{k}(z)$ denotes the Taylor expansion of $H_{k}(z)$ up to fifth
order.

Recall the result of Lemma 2.21, which tells us that

$z_{k}=w_{k}+\alpha
P_{k}(w)+\alpha^{2}[P_{k}||P](w)+\mathcal{O}(|\alpha|^{3}|w|^{4})\,.$ (S123)

Combined, and using Lemma 2.19, we get

$\displaystyle\dot{z}_{k}$
$\displaystyle=\gamma_{k}z_{k}-\beta_{k}z_{k}|z_{k}|^{2}+\alpha{H}^{5}_{k}(z)+\mathcal{O}(|\alpha||z|^{6})$
(S124) $\displaystyle=\gamma_{k}(w_{k}+\alpha
P_{k}(w)+\alpha^{2}[P_{k}||P](w))+\mathcal{O}(|\alpha|^{3}|w|^{4})$
$\displaystyle-\beta_{k}w_{k}|w_{k}|^{2}-\alpha[\beta_{k}w_{k}|w_{k}|^{2}||P](w)+\mathcal{O}(|\alpha|^{2}|w|^{5})$
$\displaystyle+\alpha{H}^{5}_{k}(w)+\alpha^{2}[{H}^{5}_{k}||P](w)+\mathcal{O}(|\alpha|^{3}|w|^{4})$
$\displaystyle+\mathcal{O}(|\alpha||w|^{6})$
$\displaystyle=\gamma_{k}w_{k}-\beta_{k}w_{k}|w_{k}|^{2}+\alpha(\gamma_{k}P_{k}(w)+{H}^{5}_{k}(w))-\alpha[\beta_{k}w_{k}|w_{k}|^{2}||P](w)$
$\displaystyle+\alpha^{2}(\gamma_{k}[P_{k}||P](w)+[{H}^{5}_{k}||P](w))+\mathcal{O}(|\alpha||w|^{6}+|\alpha|^{2}|w|^{5}+|\alpha|^{3}|w|^{4})\,.$

We will write

$\displaystyle L^{2}_{k}(w):=[\beta_{k}w_{k}|w_{k}|^{2}||P](w)\,,$ (S125)

which has lower degree $4$ or higher, to arrive at:

######  Lemma 2.23.

The term $\dot{z}_{k}$ may be expressed in the new $w$ coordinates by

$\displaystyle\dot{z}_{k}$
$\displaystyle=\gamma_{k}w_{k}-\beta_{k}w_{k}|w_{k}|^{2}$ (S126)
$\displaystyle+\alpha(\gamma_{k}P_{k}(w)+{H}^{5}_{k}(w))-\alpha L^{2}_{k}(w)$
$\displaystyle+\alpha^{2}(\gamma_{k}[P_{k}||P](w)+[{H}^{5}_{k}||P](w))$
$\displaystyle+\mathcal{O}(|\alpha||w|^{6}+|\alpha|^{2}|w|^{5}+|\alpha|^{3}|w|^{4})\,.$

###  2.3 The first reduction

We may now substitute the results of Lemma 2.22 and Lemma 2.23 into

$\displaystyle\dot{w}_{i}=\dot{z}_{i}-\alpha\frac{d}{dt}P_{i}(z)\,.$ (S127)

We obtain

$\displaystyle\dot{w}_{k}$
$\displaystyle=\gamma_{k}w_{k}-\beta_{k}w_{k}|w_{k}|^{2}$ (S128)
$\displaystyle+\alpha(\gamma_{k}P_{k}(w)+{H}^{5}_{k}(w))-\alpha L^{2}_{k}(w)$
$\displaystyle+\alpha^{2}(\gamma_{k}[P_{k}||P](w)+[{H}^{5}_{k}||P](w))$
$\displaystyle+\mathcal{O}(|\alpha||w|^{6}+|\alpha|^{2}|w|^{5}+|\alpha|^{3}|w|^{4})$
$\displaystyle-\alpha(\Gamma
P_{k}(w)-L^{1}_{k}(w)+\alpha[P_{k}||{H}^{5}](w)+\alpha[\Gamma P_{k}||P](w))$
$\displaystyle=\gamma_{k}w_{k}-\beta_{k}w_{k}|w_{k}|^{2}$
$\displaystyle+\alpha(\gamma_{k}P_{k}(w)+{H}^{5}_{k}(w)-\Gamma
P_{k}(w))+\alpha(L^{1}_{k}(w)-L^{2}_{k}(w))$
$\displaystyle+\alpha^{2}(\gamma_{k}[P_{k}||P](w)+[{H}^{5}_{k}||P](w)-[P_{k}||{H}^{5}](w)-[\Gamma
P_{k}||P](w))$
$\displaystyle+\mathcal{O}(|\alpha||w|^{6}+|\alpha|^{2}|w|^{5}+|\alpha|^{3}|w|^{4})\,.$

By Lemma 2.12 we may further write this as

$\displaystyle\dot{w}_{k}$
$\displaystyle=\gamma_{k}w_{k}-\beta_{k}w_{k}|w_{k}|^{2}$ (S129)
$\displaystyle+\alpha(\gamma_{k}P_{k}(w)+{H}^{5}_{k}(w)-\Gamma
P_{k}(w))+\alpha(L^{1}_{k}(w)-L^{2}_{k}(w))$
$\displaystyle+\alpha^{2}[\gamma_{k}P_{k}+{H}^{5}_{k}-\Gamma
P_{k}||P](w)-\alpha^{2}[P_{k}||{H}^{5}](w)$
$\displaystyle+\mathcal{O}(|\alpha||w|^{6}+|\alpha|^{2}|w|^{5}+|\alpha|^{3}|w|^{4})\,.$

Next, we claim that our choice of polynomial $P_{k}=\widehat{H}^{5}_{k}$
guarantees that the term

$\gamma_{k}P_{k}+{H}^{5}_{k}-\Gamma P_{k}$

vanishes. More precisely, we prove:

######  Lemma 2.24.

Let $Q$ be a polynomial for which the $k$th non-resonance condition is
satisfied. In particular, it follows that the corresponding modified
polynomial $\widehat{Q}_{k}$ is well-defined. We then have

$\displaystyle\gamma_{k}\widehat{Q}_{k}+Q-\Gamma\widehat{Q}_{k}=0\,.$ (S130)

###### Proof.

By definitions 2.6 and 2.8, we see that it suffices to show this when $Q$ is
given by a single monomial

$Q(z)=z_{1}^{s_{1}}\dots
z_{n}^{s_{n}}\overline{z}_{1}^{t_{1}}\dots\overline{z}_{n}^{t_{n}}\,,$

where $t_{1},\dots,t_{n},s_{1},\dots,s_{n}$ are non-negative integers. More
precisely, we use here that the maps $Q\mapsto\widehat{Q}_{k}$ and
$Q\mapsto\Gamma Q=[Q||\dots,\gamma_{j}z_{j},\dots]$ are complex linear, when
defined. By Definition 2.8 we find

$\widehat{Q}_{k}(z)=\frac{z_{1}^{s_{1}}\dots
z_{n}^{s_{n}}\overline{z}_{1}^{t_{1}}\dots\overline{z}_{n}^{t_{n}}}{s_{1}\gamma_{1}+\dots+s_{n}\gamma_{n}+t_{1}\overline{\gamma}_{1}+\dots+t_{n}\overline{\gamma}_{n}-\gamma_{k}}=\frac{Q(z)}{s_{1}\gamma_{1}+\dots+s_{n}\gamma_{n}+t_{1}\overline{\gamma}_{1}+\dots+t_{n}\overline{\gamma}_{n}-\gamma_{k}}\,.$

Example 2.7 now tells us that

$\displaystyle\Gamma\widehat{Q}_{k}(z)$
$\displaystyle=\frac{(s_{1}\gamma_{1}+\dots+s_{n}\gamma_{n}+t_{1}\overline{\gamma}_{1}+\dots+t_{n}\overline{\gamma}_{n})z_{1}^{s_{1}}\dots
z_{n}^{s_{n}}\overline{z}_{1}^{t_{1}}\dots\overline{z}_{n}^{t_{n}}}{s_{1}\gamma_{1}+\dots+s_{n}\gamma_{n}+t_{1}\overline{\gamma}_{1}+\dots+t_{n}\overline{\gamma}_{n}-\gamma_{k}}$
(S131)
$\displaystyle=\frac{(s_{1}\gamma_{1}+\dots+s_{n}\gamma_{n}+t_{1}\overline{\gamma}_{1}+\dots+t_{n}\overline{\gamma}_{n})Q(z)}{s_{1}\gamma_{1}+\dots+s_{n}\gamma_{n}+t_{1}\overline{\gamma}_{1}+\dots+t_{n}\overline{\gamma}_{n}-\gamma_{k}}\,.$

We therefore conclude that

$\displaystyle\Gamma\widehat{Q}_{k}(z)-\gamma_{k}\widehat{Q}_{k}(z)$
$\displaystyle=\frac{(s_{1}\gamma_{1}+\dots+s_{n}\gamma_{n}+t_{1}\overline{\gamma}_{1}+\dots+t_{n}\overline{\gamma}_{n})Q(z)}{s_{1}\gamma_{1}+\dots+s_{n}\gamma_{n}+t_{1}\overline{\gamma}_{1}+\dots+t_{n}\overline{\gamma}_{n}-\gamma_{k}}$
(S132)
$\displaystyle-\frac{\gamma_{k}Q(z)}{s_{1}\gamma_{1}+\dots+s_{n}\gamma_{n}+t_{1}\overline{\gamma}_{1}+\dots+t_{n}\overline{\gamma}_{n}-\gamma_{k}}$
$\displaystyle=\frac{(s_{1}\gamma_{1}+\dots+s_{n}\gamma_{n}+t_{1}\overline{\gamma}_{1}+\dots+t_{n}\overline{\gamma}_{n}-\gamma_{k})Q(z)}{s_{1}\gamma_{1}+\dots+s_{n}\gamma_{n}+t_{1}\overline{\gamma}_{1}+\dots+t_{n}\overline{\gamma}_{n}-\gamma_{k}}=Q(z)\,.$

Thus, we precisely find

$\displaystyle\gamma_{k}\widehat{Q}_{k}+Q-\Gamma\widehat{Q}_{k}=-(\Gamma\widehat{Q}_{k}-\gamma_{k}\widehat{Q}_{k})+Q=-Q+Q=0\,,$
(S133)

which completes the proof. ∎

As we have used the shorthand notation
$\widehat{H}_{k}^{5}:=\widehat{({H^{5}_{k}})}_{k}$, we see that indeed

$\displaystyle\gamma_{k}P_{k}+{H}^{5}_{k}-\Gamma
P_{k}=\gamma_{k}\widehat{H}_{k}^{5}+{H}^{5}_{k}-\Gamma\widehat{H}_{k}^{5}=0\,.$
(S134)

Returning to Equation (S129), we find that it simplifies to

$\displaystyle\dot{w}_{k}$
$\displaystyle=\gamma_{k}w_{k}-\beta_{k}w_{k}|w_{k}|^{2}+\alpha(L^{1}_{k}(w)-L^{2}_{k}(w))-\alpha^{2}[P_{k}||{H}^{5}](w)$
(S135)
$\displaystyle+\mathcal{O}(|\alpha||w|^{6}+|\alpha|^{2}|w|^{5}+|\alpha|^{3}|w|^{4})\,,$

where we recall that $L^{1}_{k}$ and $L^{2}_{k}$ are defined as

$\displaystyle L^{1}_{k}(w)$ $\displaystyle:=[P_{k}||(\dots,\beta
w_{j}|w_{j}|^{2},\dots)](w)\text{ and }$ (S136) $\displaystyle L^{2}_{k}(w)$
$\displaystyle:=[\beta_{k}w_{k}|w_{k}|^{2}||P](w)\,,$

which are both polynomials of lower degree $4$ or higher.

###  2.4 The second reduction

Next, we wish to get rid of the term $\alpha(L^{1}_{k}(w)-L^{2}_{k}(w))$ in
Equation (S135). This follows along the same lines as in the previous
reduction. We start by defining new variables

$u_{k}=w_{k}-\alpha Q_{k}(w)\,$ (S137)

where each $Q_{k}$ is a polynomial of lower degree $4$ or higher. Note that by
Lemma 2.18 we may write

$w_{k}=u_{k}+\alpha Q_{k}(u)+\mathcal{O}(|\alpha|^{2}|u|^{7})\,.$ (S138)

Using Equation (S135) we obtain

$\displaystyle\dot{u}_{k}$
$\displaystyle=\dot{w}_{k}-\alpha\frac{d}{dt}Q_{k}(w)\left|{\begin{array}[]{l}\dot{w}_{j}=\gamma_{j}w_{j}\\\
\forall\,j\in\\{1,\dots,n\\}\\\
\end{array}}\right.+\quad\mathcal{O}(|\alpha||w|^{6})$ (S141)
$\displaystyle=\dot{w}_{k}-\alpha\Gamma Q_{k}(w)+\mathcal{O}(|\alpha||w|^{6})$
$\displaystyle=\gamma_{k}w_{k}-\beta_{k}w_{k}|w_{k}|^{2}+\alpha(L^{1}_{k}(w)-L^{2}_{k}(w))-\alpha^{2}[P_{k}||{H}^{5}](w)-\alpha\Gamma
Q_{k}(w)$
$\displaystyle+\mathcal{O}(|\alpha||w|^{6}+|\alpha|^{2}|w|^{5}+|\alpha|^{3}|w|^{4})\,.$

Next, substituting $w_{k}$ by the right hand side of Equation (S207) yields

$\displaystyle\dot{u}_{k}$
$\displaystyle=\gamma_{k}w_{k}-\beta_{k}w_{k}|w_{k}|^{2}+\alpha(L^{1}_{k}(w)-L^{2}_{k}(w))-\alpha^{2}[P_{k}||{H}^{5}](w)-\alpha\Gamma
Q_{k}(w)$ (S142)
$\displaystyle+\mathcal{O}(|\alpha||w|^{6}+|\alpha|^{2}|w|^{5}+|\alpha|^{3}|w|^{4})$
$\displaystyle=\gamma_{k}u_{k}+\alpha\gamma_{k}Q_{k}(u)-\beta_{k}u_{k}|u_{k}|^{2}+\alpha(L^{1}_{k}(u)-L^{2}_{k}(u))-\alpha^{2}[P_{k}||{H}^{5}](u)-\alpha\Gamma
Q_{k}(u)$
$\displaystyle+\mathcal{O}(|\alpha||u|^{6}+|\alpha|^{2}|u|^{5}+|\alpha|^{3}|u|^{4})$
$\displaystyle=\gamma_{k}u_{k}-\beta_{k}u_{k}|u_{k}|^{2}+\alpha(L^{1}_{k}(u)-L^{2}_{k}(u)+\gamma_{k}Q_{k}(u)-\Gamma
Q_{k}(u))-\alpha^{2}[P_{k}||{H}^{5}](u)$
$\displaystyle+\mathcal{O}(|\alpha||u|^{6}+|\alpha|^{2}|u|^{5}+|\alpha|^{3}|u|^{4})\,.$

It remains to choose $Q_{k}$ such that

$L^{1}_{k}(u)-L^{2}_{k}(u)+\gamma_{k}Q_{k}(u)-\Gamma Q_{k}(u)=0\,.$ (S143)

Setting $S_{k}(u):=L^{1}_{k}(u)-L^{2}_{k}(u)$, Equation (S143) becomes

$\displaystyle\gamma_{k}Q_{k}+S_{k}-\Gamma Q_{k}=0\,,$ (S144)

which is of the same form as Equation (S130). It therefore follows from Lemma
(2.24) that a solution to Equation (S143) is given by
$Q_{k}=\widehat{S}_{k}:=\widehat{(S_{k})}_{k}$, if indeed this is well-
defined. The following lemmas show that the non-resonance conditions of
$H_{k}^{5}$ are enough to ensure $\widehat{S}_{k}$ exists.

######  Lemma 2.25.

The polynomial $S_{k}(u)=L^{1}_{k}(u)-L^{2}_{k}(u)$ may be expressed as the
sum of terms $u_{k}^{2}\overline{R(u)}$ and $|u_{j}|^{2}R(u)$ for
$j\in\\{1,\dots,n\\}$ and with $R(u)$ a monomial term appearing in $P_{k}(u)$.

###### Proof.

We start with $L^{2}_{k}(u)$. By definition, we have

$\displaystyle L^{2}_{k}(u)$
$\displaystyle=[\beta_{k}u_{k}|u_{k}|^{2}||P](u)=2\beta_{k}u_{k}\overline{u}_{k}P_{k}(u)+\beta_{k}u_{k}^{2}\overline{P_{k}(u)}$
(S145)
$\displaystyle=2\beta_{k}|u_{k}|^{2}P_{k}(u)+\beta_{k}u_{k}^{2}\overline{P_{k}(u)}\,.$

As $P_{k}(u)$ may be expressed as the sum of monomials that appear in
$P_{k}(u)$ (tautologically), we see that $L^{2}_{k}(u)$ can indeed be written
as the sum of terms $u_{k}^{2}\overline{R(u)}$ and $|u_{j}|^{2}R(u)$, with
$R(u)$ a monomial appearing in $P_{k}(u)$.

Next, recall that $L^{1}_{k}(u)$ is defined as

$L^{1}_{k}(u):=[P_{k}||(\dots,\beta_{j}u_{j}|u_{j}|^{2},\dots)](u)\,.$ (S146)

By definition of the bracket $[\bullet||\bullet]$, this means $L^{1}_{k}(u)$
is obtained from $P_{k}$ by substituting terms $u_{j}$ by
$\beta_{j}u_{j}|u_{j}|^{2}$ and terms $\overline{u}_{j}$ by
$\overline{\beta_{j}u_{j}|u_{j}|^{2}}=\overline{\beta_{j}u_{j}}|u_{j}|^{2}$.
More precisely, if $R(u)$ is a monomial term of $P_{k}(u)$ given by

$\displaystyle R(u)=u_{1}^{s_{1}}\dots
u_{n}^{s_{n}}\overline{u}_{1}^{t_{1}}\dots\overline{u}_{n}^{t_{n}}\,,$ (S147)

then we find

$\displaystyle[R(u)||(\dots,\beta_{j}u_{j}|u_{j}|^{2},\dots)](u)$ (S148)
$\displaystyle=$ $\displaystyle\sum_{j=1}^{n}s_{j}u_{1}^{s_{1}}\dots
u_{j}^{s_{j}-1}(\beta_{j}u_{j}|u_{j}|^{2})\dots
u_{n}^{s_{n}}\overline{u}_{1}^{t_{1}}\dots\overline{u}_{n}^{t_{n}}+\sum_{j=1}^{n}t_{j}u_{1}^{s_{1}}\dots
u_{n}^{s_{n}}\overline{u}_{1}^{t_{1}}\dots\overline{u}_{j}^{t_{j}-1}\overline{(\beta_{j}u_{j}|u_{j}|^{2})}\dots\overline{u}_{n}^{t_{n}}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{n}s_{j}\beta_{j}|u_{j}|^{2}u_{1}^{s_{1}}\dots
u_{j}^{s_{j}}\dots
u_{n}^{s_{n}}\overline{u}_{1}^{t_{1}}\dots\overline{u}_{n}^{t_{n}}+\sum_{j=1}^{n}t_{j}\overline{\beta}_{j}|u_{j}|^{2}u_{1}^{s_{1}}\dots
u_{n}^{s_{n}}\overline{u}_{1}^{t_{1}}\dots\overline{u}_{j}^{t_{j}}\dots\overline{u}_{n}^{t_{n}}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{n}s_{j}\beta_{j}|u_{j}|^{2}R(u)+\sum_{j=1}^{n}t_{j}\overline{\beta}_{j}|u_{j}|^{2}R(u)\,.$

Hence, by linearity of $[\bullet||\bullet]$ in the first slot (see Lemma
2.12), we see that $L^{1}_{k}(u)$ is again of the right form.

It follows that $S_{k}(u)=L^{1}_{k}(u)-L^{2}_{k}(u)$ can indeed be expressed
as a sum of the given monomials terms. This completes the proof. ∎

######  Lemma 2.26.

Let $R(u)$ be a monomial and let $k,j\in\\{1,\dots,n\\}$ be fixed indices. The
$k$th non-resonance condition of $R(u)$ is satisfied if and only if the $k$th
non-resonance condition of $|u_{j}|^{2}R(u)$ is satisfied, if and only if the
$k$th non-resonance condition of $u_{k}^{2}\overline{R(u)}$ is satisfied.

###### Proof.

We write

$\displaystyle R(u)=u_{1}^{s_{1}}\dots
u_{n}^{s_{n}}\overline{u}_{1}^{t_{1}}\dots\overline{u}_{n}^{t_{n}}\,,$ (S149)

so that the $k$th non-resonance condition of $R(u)$ is given by

$\displaystyle
s_{1}\omega_{1}+\dots+{s_{n}}\omega_{n}-t_{1}\omega_{1}-\dots-{t_{n}}\omega_{n}-\omega_{k}\not=0\,.$
(S150)

It follows that the $k$th non-resonance condition of $|u_{j}|^{2}R(u)$ is
given by

$\displaystyle\omega_{j}-\omega_{j}+$ $\displaystyle
s_{1}\omega_{1}+\dots+{s_{n}}\omega_{n}-t_{1}\omega_{1}-\dots-{t_{n}}\omega_{n}-\omega_{k}$
(S151) $\displaystyle=$ $\displaystyle
s_{1}\omega_{1}+\dots+{s_{n}}\omega_{n}-t_{1}\omega_{1}-\dots-{t_{n}}\omega_{n}-\omega_{k}\not=0\,,$

which coincides with that of $R(u)$.
Likewise, the $k$th non-resonance condition of $u_{k}^{2}\overline{R(u)}$ is
given by

$\displaystyle
2\omega_{k}-s_{1}\omega_{1}-\dots-{s_{n}}\omega_{n}+t_{1}\omega_{1}+\dots+{t_{n}}\omega_{n}-\omega_{k}$
(S152) $\displaystyle=-$ $\displaystyle
s_{1}\omega_{1}-\dots-{s_{n}}\omega_{n}+t_{1}\omega_{1}+\dots+{t_{n}}\omega_{n}+\omega_{k}$
$\displaystyle=-$
$\displaystyle(s_{1}\omega_{1}+\dots+{s_{n}}\omega_{n}-t_{1}\omega_{1}-\dots-{t_{n}}\omega_{n}-\omega_{k})\not=0\,,$

which is again equivalent to Equation (S150). This completes the proof. ∎

Lemmas 2.25 and 2.26 guarantee that the $k$th non-resonance condition of
$S_{k}(u)=L^{1}_{k}(u)-L^{2}_{k}(u)$ is satisfied if the $k$th non-resonance
condition of $P_{k}=\widehat{H}^{5}_{k}$ is satisfied. As
$\widehat{H}^{5}_{k}$ has the same monomial terms as $H^{5}_{k}$ (though
rescaled), we see that the $k$th non-resonance conditions of $S_{k}(u)$ are
indeed satisfied. Therefore, a solution to Equation (S143) exists by Lemma
(2.24) and may be given by $Q_{k}=\widehat{S}_{k}$. Note that this choice of
$Q_{k}$ has lower degree $4$ or higher, as we assumed throughout.

Returning to Equation (S142), we finally arrive at

$\displaystyle\dot{u}_{k}=\gamma_{k}u_{k}-\beta_{k}u_{k}|u_{k}|^{2}-\alpha^{2}[P_{k}||{H}^{5}](u)+\mathcal{O}(|\alpha||u|^{6}+|\alpha|^{2}|u|^{5}+|\alpha|^{3}|u|^{4})\,.$
(S153)

We have therefore shown:

###### Proof of Theorem 2.4 and Proposition 2.9.

The calculations in this section show that the successive coordinate
transformations

$\displaystyle w_{k}$ $\displaystyle=z_{k}-\alpha P_{k}(z)$ (S154)
$\displaystyle u_{k}$ $\displaystyle=w_{k}-\alpha Q_{k}(w)\,$

bring the ODE

$\dot{z}_{k}=\gamma_{k}z_{k}-\beta_{k}z_{k}|z_{k}|^{2}+\alpha H_{k}(z)$

into the form

$\dot{u}_{k}=\gamma_{k}u_{k}-\beta_{k}u_{k}|u_{k}|^{2}-\alpha^{2}[P_{k}||{H}^{5}](u)+\mathcal{O}(|\alpha||u|^{6}+|\alpha|^{2}|u|^{5}+|\alpha|^{3}|u|^{4})\,.$

If we now set

$G_{k}(u):=[P_{k}||{H}^{5}]=[\widehat{H}^{5}_{k}||{H}^{5}]\,,$

then we indeed get

$\dot{u}_{k}=\gamma_{k}u_{k}-\beta_{k}u_{k}|u_{k}|^{2}-\alpha^{2}G_{k}(u)+\mathcal{O}(|\alpha||u|^{6}+|\alpha|^{2}|u|^{5}+|\alpha|^{3}|u|^{4})\,.$

This completes the proof. ∎

##  3 Anomalous synchronization on a 4-node ring

Consider the four node ring network with a coupling function
$h(z,w)=z\bar{w}.$ leading to

$\displaystyle\dot{z}_{1}$ $\displaystyle=\gamma_{1}z_{1}-\beta
z_{1}|z_{1}|^{2}+\alpha(z_{1}\overline{z}_{2}+z_{1}\overline{z}_{4})$ (S155)
$\displaystyle\dot{z}_{2}$ $\displaystyle=\gamma_{2}z_{2}-\beta
z_{2}|z_{2}|^{2}+\alpha(z_{2}\overline{z}_{3}+z_{2}\overline{z}_{1})$
$\displaystyle\dot{z}_{3}$ $\displaystyle=\gamma_{3}z_{3}-\beta
z_{3}|z_{3}|^{2}+\alpha(z_{3}\overline{z}_{4}+z_{3}\overline{z}_{2})$
$\displaystyle\dot{z}_{4}$ $\displaystyle=\gamma_{4}z_{4}-\beta
z_{4}|z_{4}|^{2}+\alpha(z_{4}\overline{z}_{1}+z_{4}\overline{z}_{3})\,$

We set the parameters $\beta_{k}=-1$, $\gamma_{k}=\lambda+i\omega_{k}$,
$\lambda=1$ and $\omega_{1}=1+\delta$, $\omega_{2}=1$ and $\omega_{3}=5$ and
$\omega_{4}=6$ for performing the simulation of Eq. (S155). We then vary the
mismatch $\delta$ and coupling $\alpha$. Notice that by a naive inspection of
the original equations we obtain

$\dot{\theta}_{k}=\omega_{k}+\alpha\sum_{\ell=1}^{4}A_{k\ell}\sin\theta_{\ell}$

thus, instead of a diffusive interaction we would obtain a interaction akin to
forcing [37]. For each set of $\alpha$ and $\delta$ values, we simulate the
network for 50000s with 0.01 time step. We remove first 10000s from as
transient and compute the unwrapped phases. As we are interested in the phase
synchronization, we introduce a new variable for the phase differences,
$\phi=\theta_{1}-\theta_{2}$, and a naive calculation leads to

$\dot{\phi}=\delta+\alpha[\sin\theta_{2}+\sin\theta_{4}-\sin\theta_{1}-\sin\theta_{3}]+O(\alpha^{2}).$

Because $\omega_{2},\omega_{3}\ll\omega_{1}$, we can average over the fast
oscillations and neglect contributions from these phases then
$\theta_{2}=\theta_{1}+\phi$ we would obtain a interaction term as
$\alpha[\sin\theta_{1}\cos\phi+\cos\theta_{1}\sin\phi]+O(\alpha^{2}),$
however, since $\theta_{1}$ is a fast variable for $\phi$ we also average over
$\theta_{1}$. Thus the whole interaction term linear in $\alpha$ vanishes.

We calculate the mean synchronization error as

$\displaystyle E=\frac{1}{T}\sum_{t=1}^{T}|\phi(t)|.$ (S156)

The synchronization error for varying $\delta=-0.2$ to $0.2$ with $0.01$ step
size and $\alpha=0.0$ to $0.5$ with step size $0.025$, we observed a
synchronization tongue scales with $\alpha\propto\sqrt{\delta}$ (Supplementary
Fig. 7).

Figure 7: Arnold tongue represents the phase synchronization between nodes 1
and 2. The tongue is illustrated for varying mismatch $\delta$ and coupling
strength $\alpha$. As the absolute value of $|\delta|\ll 1$, synchronization
occurs for small coupling strength. However, if the mismatch is high then also
stronger coupling is needed to emerge synchronization. The borderline between
synchrony and asynchrony scales with $\alpha\propto\sqrt{\delta}$.

###  3.1 Explanation of anomalous synchronization via normal form theory

Recall that we may bring this ODE into the form

$\displaystyle\dot{u}_{1}$ $\displaystyle=\gamma_{1}u_{1}-\beta
u_{1}|u_{1}|^{2}$ (S157)
$\displaystyle-\alpha^{2}u_{1}\left(\frac{\overline{u}_{2}\overline{u}_{4}}{\overline{\gamma}_{2}}+\frac{\overline{u}_{4}\overline{u}_{2}}{\overline{\gamma}_{4}}+\frac{\overline{u}_{2}\overline{u}_{2}}{\overline{\gamma}_{2}}+\frac{\overline{u}_{4}\overline{u}_{4}}{\overline{\gamma}_{4}}+\frac{\overline{u}_{2}u_{3}}{\overline{\gamma}_{2}}+\frac{\overline{u}_{2}u_{1}}{\overline{\gamma}_{2}}+\frac{\overline{u}_{4}u_{3}}{\overline{\gamma}_{4}}+\frac{\overline{u}_{4}u_{1}}{\overline{\gamma}_{4}}\right)+\text{h.o.t.}$
$\displaystyle\dot{u}_{2}$ $\displaystyle=\gamma_{2}u_{2}-\beta
u_{2}|u_{2}|^{2}$
$\displaystyle-\alpha^{2}u_{2}\left(\frac{\overline{u}_{1}\overline{u}_{3}}{\overline{\gamma}_{1}}+\frac{\overline{u}_{3}\overline{u}_{1}}{\overline{\gamma}_{3}}+\frac{\overline{u}_{1}\overline{u}_{1}}{\overline{\gamma}_{1}}+\frac{\overline{u}_{3}\overline{u}_{3}}{\overline{\gamma}_{3}}+\frac{\overline{u}_{3}u_{4}}{\overline{\gamma}_{3}}+\frac{\overline{u}_{3}u_{2}}{\overline{\gamma}_{3}}+\frac{\overline{u}_{1}u_{4}}{\overline{\gamma}_{1}}+\frac{\overline{u}_{1}u_{2}}{\overline{\gamma}_{1}}\right)+\text{h.o.t.}$
$\displaystyle\dot{u}_{3}$ $\displaystyle=\gamma_{3}u_{3}-\beta
u_{3}|u_{3}|^{2}$
$\displaystyle-\alpha^{2}u_{3}\left(\frac{\overline{u}_{2}\overline{u}_{4}}{\overline{\gamma}_{2}}+\frac{\overline{u}_{4}\overline{u}_{2}}{\overline{\gamma}_{4}}+\frac{\overline{u}_{2}\overline{u}_{2}}{\overline{\gamma}_{2}}+\frac{\overline{u}_{4}\overline{u}_{4}}{\overline{\gamma}_{4}}+\frac{\overline{u}_{2}u_{3}}{\overline{\gamma}_{2}}+\frac{\overline{u}_{2}u_{1}}{\overline{\gamma}_{2}}+\frac{\overline{u}_{4}u_{3}}{\overline{\gamma}_{4}}+\frac{\overline{u}_{4}u_{1}}{\overline{\gamma}_{4}}\right)+\text{h.o.t.}$
$\displaystyle\dot{u}_{4}$ $\displaystyle=\gamma_{4}u_{4}-\beta
u_{4}|u_{4}|^{2}$
$\displaystyle-\alpha^{2}u_{4}\left(\frac{\overline{u}_{1}\overline{u}_{3}}{\overline{\gamma}_{1}}+\frac{\overline{u}_{3}\overline{u}_{1}}{\overline{\gamma}_{3}}+\frac{\overline{u}_{1}\overline{u}_{1}}{\overline{\gamma}_{1}}+\frac{\overline{u}_{3}\overline{u}_{3}}{\overline{\gamma}_{3}}+\frac{\overline{u}_{3}u_{4}}{\overline{\gamma}_{3}}+\frac{\overline{u}_{3}u_{2}}{\overline{\gamma}_{3}}+\frac{\overline{u}_{1}u_{4}}{\overline{\gamma}_{1}}+\frac{\overline{u}_{1}u_{2}}{\overline{\gamma}_{1}}\right)+\text{h.o.t.}$

We will assume that $\gamma_{1}\approx\gamma_{2}$. In particular, we consider
the possibility that $\gamma_{1}=\gamma_{2}$. In that case we have
$\gamma_{1}+\overline{\gamma_{2}}=2\operatorname{Re}(\gamma_{1})$, which may
be arbitrarily small. Therefore, we may not assume that
$\gamma_{1}+\overline{\gamma_{2}}\not=0$. Other than this, there are no
relevant restrictions. I.e.,
$\gamma_{1},\dots,\gamma_{4},\gamma_{2}-\overline{\gamma_{3}}\dots\gamma_{3}-\overline{\gamma_{4}}$
are all sufficiently large. It follows that we may bring equation (S157) into
the form

$\displaystyle\dot{v}_{1}$ $\displaystyle=\gamma_{1}v_{1}-\beta
v_{1}|v_{1}|^{2}-\epsilon\frac{v_{1}\overline{v}_{2}v_{1}}{\overline{\gamma}_{2}}+\mathcal{O}(|\epsilon,v|^{5})$
(S158) $\displaystyle\dot{v}_{2}$ $\displaystyle=\gamma_{2}v_{2}-\beta
v_{2}|v_{2}|^{2}-\epsilon\frac{v_{2}\overline{v}_{1}v_{2}}{\overline{\gamma}_{1}}+\mathcal{O}(|\epsilon,v|^{5})$
$\displaystyle\dot{v}_{3}$ $\displaystyle=\gamma_{3}v_{3}-\beta
v_{3}|v_{3}|^{2}-\epsilon\frac{v_{3}\overline{v}_{2}v_{1}}{\overline{\gamma}_{2}}+\mathcal{O}(|\epsilon,v|^{5})$
$\displaystyle\dot{v}_{4}$ $\displaystyle=\gamma_{4}v_{4}-\beta
v_{4}|v_{4}|^{2}-\epsilon\frac{v_{4}\overline{v}_{1}v_{2}}{\overline{\gamma}_{1}}+\mathcal{O}(|\epsilon,v|^{5})\,,$

where $\epsilon=\alpha^{2}$. Note that the network topology has changed
drastically. Moreover, performing a phase reduction and introducing the phase
difference $\psi=\theta_{1}-\theta_{2}$ we obtain

$\dot{\psi}=\delta-c\alpha^{2}\sin\psi$

where $c$ is a constant depending on $\gamma_{1}$. By analyzing the fixed
points of this equation we obtain the synchronization tongue behavior where
the critical coupling $\alpha_{c}$ for synchronization scales as
$\sqrt{\delta}$.

##  4 Phase reduction for $h=(z^{2}+z)\bar{w}$ and resonance
$\omega_{1}-\omega_{2,4}+\omega_{3}=0$

For simplicity we fix $\beta_{k}=-1$ and obtain
$\gamma_{k}=r_{0}^{2}+i\omega_{k}.$ We also introduce
$\Delta_{pq}=\omega_{p}-\omega_{q}$. Note that $r(t)=r_{0}+h.o.t$ along with
$\bar{z}\dot{z}=ir_{0}^{2}\dot{\theta}+h.o.t.$ and

$\displaystyle\frac{1}{\gamma_{p}+\bar{\gamma}_{q}}$ $\displaystyle=$
$\displaystyle\frac{2r_{0}^{2}-i\Delta_{pq}}{4r_{0}^{4}+\Delta_{pq}^{2}}$
(S159)

Replacing these observations into Eq. (6) of the main manuscript and
performing the reduction we obtain the functions in Eq. (7-8) of the main
manuscript as

$\displaystyle\rho_{pq}(\phi)$ $\displaystyle=$
$\displaystyle-\frac{\Delta_{pq}}{4r_{0}^{4}+\Delta_{pq}^{2}}\cos\phi+\frac{2r_{0}^{2}}{4r_{0}^{4}+\Delta_{pq}^{2}}\sin\phi$
$\displaystyle\sigma_{pqr}(\phi)$ $\displaystyle=$
$\displaystyle-\chi_{pqr}\sin\phi+\upsilon_{pqr}\cos\phi$

where

$\chi_{pqr}=r_{0}^{2}\left(\frac{4}{4r_{0}^{4}+\Delta_{pq}^{2}}+\frac{4}{4r_{0}^{4}+\Delta_{pr}^{2}}+\frac{1}{r_{0}^{4}+\omega_{q}^{2}}+\frac{1}{r_{0}^{4}+\omega_{r}^{2}}\right)$

and

$\upsilon_{pqr}=-\frac{2\Delta_{pq}}{4r_{0}^{4}+\Delta_{pq}^{2}}-\frac{2\Delta_{pr}}{4r_{0}^{4}+\Delta_{pr}^{2}}+\frac{\omega_{q}}{r_{0}^{4}+\omega_{q}^{2}}+\frac{\omega_{r}}{r_{0}^{4}+\omega_{r}^{2}}$

##  5 Slow phase dynamics for $h=(z^{2}+z)\bar{w}$ and resonance
$\omega_{1}-\omega_{2,4}+\omega_{3}=0$

We consider networks of $n$ coupled oscillators

$\dot{z}_{k}=f_{k}(z_{k})+\alpha\sum_{\ell=1}^{n}A_{k\ell}h_{k}(z_{k},z_{\ell})$
(S161)

where $z_{k}\in\mathbb{C}$ is the state of the $k$th oscillator,
$f_{k}:\mathbb{C}\rightarrow\mathbb{C}$ is its isolated vector field,
$h_{k}:\mathbb{C}\times\mathbb{C}\rightarrow\mathbb{C}$ is the pairwise
coupling function, $\bm{A}=(A_{ij})_{i,j=1}^{n}$ is the adjacency matrix
describing the network structure, and $\alpha>0$ is the coupling strength. We
then generate a multivariate time series for a four-node ring network, as
illustrated in Figure 8 (a), with nonlinear pairwise coupling function

$h(z,w)=(z+z^{2})\bar{w}.$ (S162)

Figure 8: Emergent hypernetworks with triplet interactions. (a) The original
ring network. (b) Hypernetwork recovered from original dynamics data. The
state of a given node is influenced by triplet interactions. We show that as
along as some sparsity is imposed when the model is obtained from data, only
hypernetworks can be found. Our theory predicts the emergence of the
hypernetwork determined by the original coupling function $h$, the network,
and the resonance relations among the isolated frequencies.

We fix $\lambda=0.15$, $\omega_{1}=1.01$, $\omega_{2}=2.5$, $\omega_{3}=1.5$,
$\omega_{4}=2.49$, and $\alpha=0.18$. Numerical integration of complex
differential equations is used to solve the differential equations for
10000-sec with 0.01-sec time-step. We discard the first 5000-sec points as
transient, and we obtain a multivariate time series
$\\{z_{1}(t),z_{2}(t),z_{3}(t),z_{4}(t)\\}_{t=1}^{5000}$.

Next, we aim at obtaining a model from the multivariate time series of $z$.
Because $\alpha$ is small and the isolated orbit is exponentially stable, the
amplitude of each time series is slightly affected by
$|z_{k}|\approx\sqrt{\lambda}+O(\alpha)$ as illustrated in Figure. 9 (a), and
the dynamics is captured by the phases $\theta_{k}(t)$ of $z_{k}(t)$.
Therefore, we perform a polar decomposition
$z_{k}(t)=r_{k}(t)e^{i\theta_{k}(t)}$ to get the unwrapped phase of each time
series and obtain governing equations of the model from its phase dynamics.

Each phase $\theta_{k}$ has a frequency close to $\omega_{k}$, as illustrated
in Figure 9 (b). This means that the growth of the phases is almost linear
with coupling terms as perturbations. In fact, the coupling terms generically
contain fast variables such as phases $\theta$’s and slow variables involving
the resonant combinations of phases such as $\theta_{1}-\theta_{2}+\theta_{3}$
that change slowly in time. Therefore, we subtract the linear growth of the
phases to analyse the effects of the coupling. To this end, we introduce

$\vartheta_{k}(t)=\theta_{k}(t)-\Omega_{k}t,$ (S163)

where $\Omega_{k}$ is obtained from data under the resonance condition
$\Omega_{1}-\Omega_{i}+\Omega_{3}=0$, with $i=2,4$. In the new phases the
coupling has the same magnitude as the frequency mismatch
$\omega_{1}-\omega_{i}+\omega_{3}$, with $i=2,4$. Finally, we obtain a model
for $\vartheta_{k}$. We assume the model

$\dot{\vartheta}_{k}=\varepsilon_{k}+H_{k}(\vartheta_{1},\vartheta_{2},\vartheta_{3},\vartheta_{4})$

where
$H_{k}=\sum[c^{k}_{p}\sin\vartheta_{p}+d_{p}^{k}\cos\vartheta_{p}]+\sum[c^{k}_{p,q}\sin(\vartheta_{p}-\vartheta_{q})+d^{k}_{p,q}\cos(\vartheta_{p}-\vartheta_{q}]+\sum[c^{k}_{pq}\sin(\vartheta_{p}+\vartheta_{q}-\vartheta_{k})+d^{k}_{pq}\cos(\vartheta_{p}+\vartheta_{q}-\vartheta_{k})]$.
Note that this includes pairwise and triplet interactions. We solve for the
coefficients to obtain the least square approximation and we impose sparsity
by eliminating coefficients below a threshold $\tau=10^{-4}$. The technique is
discussed along with the package to perform the recovery as discussed in the
main text. The model recovery yields

$\displaystyle\dot{\vartheta}_{1,3}$ $\displaystyle=$
$\displaystyle\varepsilon_{1,3}+r_{1,3}(\vartheta_{1},\vartheta_{2},\vartheta_{3})+s_{1,3}(\vartheta_{1},\vartheta_{4},\vartheta_{3})$
(S164) $\displaystyle\dot{\vartheta}_{2,4}$ $\displaystyle=$
$\displaystyle\varepsilon_{2,4}+r_{2,4}(\vartheta_{1},\vartheta_{2,4},\vartheta_{3})$

where $s$ and $r$ correspond to triplets in $H_{k}$ with nonzero coefficients.

At first sight, the model recovery with triplets is remarkable because the
original equations have only pairwise interactions. Nonetheless, a
hypernetwork describes the data Figure 8 (b). We show that when $\lambda\ll 1$
and $\alpha\ll 1$ recovering a hypernetwork from data is not a coincidence. As
long as the coupling $h$ is nonlinear, by measuring the original variables of
Eq. (S161) and performing a sparse model recovery only hypernetworks can be
found as they are normal forms of the original equations.

Figure 9: Time series and emergent hypernetwork prediction. The simulation was
performed on a ring network (see Supplementary Fig. 8 (a)). (a) Amplitudes
(solid lines) are slightly affected by the coupling and remaining close to a
circle with radius $r=\sqrt{\lambda}$ (dashed circle). (b) Unwrapped phases
$\theta_{i}$ growth. (c) Time series of the slow phase $\phi_{1}$ from data
(solid) and the prediction of the emergent hypernetwork (dashed) capturing
higher-order interactions (see Supplementary Fig. 8 (b)).

###  5.1 Emergent hypernetwork predicts data behaviour

To illustrate prediction capabilities of emergent hypernetworks, we introduce
the slow phases

$\displaystyle\phi_{1}$ $\displaystyle=$
$\displaystyle\theta_{1}-\theta_{2}+\theta_{3},$ (S165)
$\displaystyle\phi_{2}$ $\displaystyle=$
$\displaystyle\theta_{1}-\theta_{4}+\theta_{3},$

where the coupling strength is comparable to the frequency mismatch
$\omega_{1}-\omega_{2}+\omega_{3}$. Our normal form theory predicts the
emergent hypernetwork phase dynamics described in Eq. (15). We obtain the
vector fields of the coupled slow phases $\phi_{1}$ and $\phi_{2}$
analytically as described in the section above.

Next, we simulate the vector fields obtained from first principles using an
adaptative Runge-Kutta method of 4th order. We treat the initial condition as
unknown and perform a optimization to obtain the the initial condition that
provides the minimum least square error between the data of the slow phase and
the simulations and theory. In Figure 9 c), we compare our predictions and
slow phases estimated from data. The theoretical prediction is in excellent
agreement with the data with an error in the prediction of less than $5\%$ per
cycle of the slow phase.

##  6 Model recovery of a 3-path with coupling $h=(z^{2}+z)\bar{w}$ and
resonance $\omega_{1}-\omega_{2}+\omega_{3}=0$

Now we consider the model for 3-nodes on a chain (Supplementary Fig. 10) that
reads as We consider the network ODE

$\displaystyle\dot{z}_{1}$ $\displaystyle=\gamma_{1}z_{1}-\beta
z_{1}|z_{1}|^{2}+\alpha(z_{1}\overline{z_{2}}+z_{1}^{2}\overline{z_{2}})$
(S166) $\displaystyle\dot{z}_{2}$ $\displaystyle=\gamma_{2}z_{2}-\beta
z_{2}|z_{2}|^{2}+\alpha([z_{2}\overline{z_{1}}+z_{2}^{2}\overline{z_{1}}]+[z_{2}\overline{z_{3}}+z_{2}^{2}\overline{z_{3}}])$
$\displaystyle\dot{z}_{3}$ $\displaystyle=\gamma_{3}z_{3}-\beta
z_{3}|z_{3}|^{2}+\alpha(z_{3}\overline{z_{2}}+z_{3}^{2}\overline{z_{2}})\,,$

Figure 10: Emergent hypernetworks with triplet interaction. Inset a) shows the
original chain network. Each isolated node dynamics is close to a Hopf
bifurcation. The pairwise coupling function $h$ is nonlinear given with
$h=(z^{2}+\bar{z})w$. Inset b) shows the hypernetwork learnt from the phase
dynamics data of the original dynamics. The state of a given node is
influenced by the triple interaction of its own state in combination of the
incoming links. Our theory also predicts the emergent of such hypernetwork.
The hypernetwork emerges as a combination of the original coupling function
$h$, the network, and the resonance relations of the isolated frequencies.

To integrate Eq. (S166) of the main manuscript for 3-node chain with
$\omega_{1}=1.01$, $\omega_{2}=2.5$ and $\omega_{3}=1.5$, we employed a
wrapper of ODEPACK routine. Numerical integration for $\alpha=0.18$ and
$\delta=0.01$ was performed for 10000s with 0.01 time step. We discard the
first 5000s points as transient. Using the simulated phases $\theta_{i}$, we
introduce new phases in Eq. (3) of the main manuscript where $\Omega_{1}=1.0$,
$\Omega_{2}=2.5$, $\Omega_{3}=1.5$. Applying the sequential thresholded least-
squares method on these new phases $\vartheta_{i}$ with thresholding parameter
$\lambda=10^{-4}$ we obtain

$\displaystyle\dot{\vartheta}_{1}$
$\displaystyle=\,\,\,\,\,\,\,\,0.01-0.001\cos(\vartheta_{1}-\vartheta_{2}+\vartheta_{3})$
(S167) $\displaystyle\dot{\vartheta}_{2}$
$\displaystyle=-0.001+0.005\cos(\vartheta_{1}-\vartheta_{2}+\vartheta_{3})$
(S168) $\displaystyle\dot{\vartheta}_{3}$
$\displaystyle=-0.001-0.001\cos(\vartheta_{1}-\vartheta_{2}+\vartheta_{3})$
(S169)

Because the norms of the functions $r_{k}$ and $s_{k}$ are small, we introduce
the slow phases

$\displaystyle\phi$ $\displaystyle=$
$\displaystyle\theta_{1}-\theta_{2}+\theta_{3}$

We then also perform a reconstruction for the slow phases $\phi$ using the
same method and obtain

$\displaystyle\dot{\phi}=0.010+0.001\sin(\phi)-0.006\cos(\phi)$ (S170)

We show the model prediction and data for the slow phase in Figure 11.

Figure 11:  New slow phase variable $\phi$ (blue curve) were computed from the
data collected from the simulations of Eq. (1) on a path. Then Eq. S198
(orange dashed curve) reconstructed from data of $\phi$ using sequential
thresholded least squares method.

###  6.1 Emergent network explanation

The normal-form for this system is given by

$\displaystyle\dot{u}_{1}$ $\displaystyle=\gamma_{1}u_{1}-\beta
u_{1}|u_{1}|^{2}-\alpha^{2}\left(\frac{1}{\gamma_{1}+\overline{\gamma}_{2}}\right)u_{1}^{2}\overline{u}_{2}u_{3}$
(S171) $\displaystyle\dot{u}_{2}$ $\displaystyle=\gamma_{2}u_{2}-\beta
u_{2}|u_{2}|^{2}-\alpha^{2}\left(\frac{2}{\gamma_{2}+\overline{\gamma}_{3}}+\frac{2}{\gamma_{2}+\overline{\gamma}_{1}}+\frac{1}{\overline{\gamma}_{3}}+\frac{1}{\overline{\gamma}_{1}}\right)u_{2}^{2}\overline{u}_{1}\overline{u}_{3}$
$\displaystyle\dot{u}_{3}$ $\displaystyle=\gamma_{3}u_{3}-\beta
u_{3}|u_{3}|^{2}-\alpha^{2}\left(\frac{1}{\gamma_{3}+\overline{\gamma}_{2}}\right)u_{3}^{2}\overline{u}_{2}u_{1}\,.$

The interaction now becomes forth order in $u$. A phase reduction leads to the
triplet interaction recovered numerically.

##  7 6 nodes network examples with $h(z,w)=z\bar{w}$

We consider the network presented in Figure 12a) with the coupling
$h(z,w)=z\bar{w}$ leading to

Figure 12: Emergent hypernetworks with triplet interaction. Inset a) shows the