High-order phase reduction for coupled oscillators

We first briefly recall the basic ideas of the phase reduction in the first order in the coupling parameter ɛ (see references [16, 18] for more details). One starts by considering an oscillatory dynamical system

$$\begin{equation*}\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t}=\mathbf{f}\left(\mathbf{y}\right)\end{equation*}$$

possessing a stable limit cycle Y(t) = Y(t + T) in its state space. As is well-known, on the limit cycle and in its basin of attraction one can define the phase φ = Φ(y) which grows uniformly in time

$$\begin{equation*}\dot {\varphi }=\omega =\frac{\partial {\Phi}}{\partial \mathbf{y}}\mathbf{f}\left(\mathbf{y}\right)\end{equation*}$$

with basic frequency ω = 2π/T. Notice, that on the limit cycle the system's state is uniquely defined by the phase: y = Y(φ). Outside of the limit cycle, this is not true: one also has to know the deviation from the cycle.

In the context of oscillatory networks, one considers many coupled oscillators that we label by index k:

$$\begin{equation}\frac{\mathrm{d}{\mathbf{y}}_{k}}{\mathrm{d}t}={\mathbf{f}}_{k}\left({\mathbf{y}}_{k}\right)+\varepsilon {\mathbf{G}}_{k}\left({\mathbf{y}}_{1},{\mathbf{y}}_{2},\dots \right).\end{equation} \tag{ 1 }$$

Here ɛ is the small parameter governing the strength of the coupling. The equation for the phases is obtained by exploiting their definition φ_k = Φ_k (y_k ):

$$\begin{equation}{\dot {\varphi }}_{k}={\omega }_{k}+\varepsilon \frac{\partial {{\Phi}}_{k}}{\partial {\mathbf{y}}_{k}}{\mathbf{G}}_{k}\left({\mathbf{y}}_{1},{\mathbf{y}}_{2},\dots \right).\end{equation} \tag{ 2 }$$

This equation is exact, but it contains the full state space trajectories {y_j }. However, for a small perturbation these trajectories are close to the limit cycle, up to deviations of order ∼ɛ. Thus, substituting the zero-order approximation y_j ≈ Y_j into the rhs of (2), we obtain equations, where only the states on the limit cycles appear, and these states are unambiguously determined by the phases:

$$\begin{equation*}{\dot {\varphi }}_{k}={\omega }_{k}+\varepsilon \frac{\partial {{\Phi}}_{k}}{\partial {\mathbf{Y}}_{k}}{\mathbf{G}}_{k}\left({\mathbf{Y}}_{1},{\mathbf{Y}}_{2},\dots \right)={\omega }_{k}+\varepsilon {G}_{k}\left({\varphi }_{k},{\varphi }_{1},{\varphi }_{2},\dots \right).\end{equation*}$$

As is clear from the presented consideration, to extend the reduction beyond the first-order approximation, one has to calculate the deviations from the limit cycle of orders ∼ɛ, ɛ², ... and to express these deviations in terms of the phases. Below, we do not provide a general approach valid for arbitrary systems, but restrict ourselves to the simplest case where the phase can be introduced analytically.

It is instructive to introduce the Stuart–Landau oscillators in dimensional variables, to understand the meaning of dimensionless parameters in the problem. We write the equation for the complex amplitude of oscillations Z for a particular oscillator as

$$\begin{equation}\frac{\mathrm{d}Z}{\mathrm{d}\tau }=\left(\mu +\text{i}\nu \right)Z-\left(\beta +\text{i}\gamma \right)\vert Z{\vert }^{2}Z+{\epsilon}G\left({Z}_{1},{Z}_{2},\dots \right),\end{equation} \tag{ 3 }$$

where G is the coupling term depending on the states of other elements of the network. Notice that for brevity and without loss of generality we omit the index for the considered oscillator and label other units by the index k = 1, 2, ....

In equation (3), all the parameters are dimensional, and it is convenient to perform a transformation to dimensionless variables and parameters. Here one should also take into account, that generally the parameters might be different for different oscillators. Using only local parameters, one can introduce the dimensionless amplitude $Z=\sqrt{\mu /\beta }A$, so that all uncoupled oscillators will have amplitude one. Another variable that we can scale is the time, and here one has to choose some common scaling for all oscillators. It appears convenient to assume that the growth rate of linear oscillations μ is the same for all oscillators (the results can be straightforwardly generalized also to the case of different μ), and use it to introduce dimensionless time as t = μτ. Then we obtain

$$\begin{equation}\dot {A}=\left(1+\text{i}\omega \right)A-\vert A{\vert }^{2}A-\text{i}\alpha A\left(\vert A{\vert }^{2}-1\right)+\varepsilon G\left({A}_{1},{A}_{2},\dots \right),\end{equation} \tag{ 4 }$$

where $\dot {A}$ now means derivative with respect to t. Here ω = ν/μ − γ/β is the dimensionless frequency of the limit cycle oscillation and the limit cycle has unit amplitude |A⁽⁰⁾| = 1. Parameter α = γ/β measures non-isochronicity, as we will see below, it determines the phase definition function Φ(A). The dimensionless coupling parameter ɛ depends on the scaling of the function G. If G contains first powers of the amplitudes Z, then ɛ = /μ, where we assume that μ/β ≈ μ_k /β_k , i.e. all the coupled oscillators Z have similar amplitudes. As we see, all the dimensionless parameters entering the problem, namely ω, α, ɛ, are the original parameters of the system normalized by the linear growth rate of oscillations μ.

For further analysis it is instructive to write the Stuart–Landau equations in polar coordinates, with the amplitude R and the angle θ, so that A = R e^iθ and

$$\begin{equation}\dot {R}=R-{R}^{3}+\varepsilon \enspace \mathrm{R}\mathrm{e}\left[{\text{e}}^{-\text{i}\theta }G\right],\end{equation} \tag{ 5 }$$

$$\begin{equation}\dot {\theta }=\omega -\alpha \left({R}^{2}-1\right)+\varepsilon {R}^{-1}\enspace \mathrm{I}\mathrm{m}\left[{\text{e}}^{-\text{i}\theta }G\right].\end{equation} \tag{ 6 }$$

It is straightforward to check that, in the absence of coupling, the quantity φ = θ − α  ln  R grows uniformly in time and, hence, is the true phase. Therefore, we rewrite equations (5) and (6) in terms of R, φ:

$$\begin{equation}\dot {R}=R-{R}^{3}+\varepsilon \enspace \mathrm{R}\mathrm{e}\left[{\text{e}}^{-\text{i}\left(\varphi +\alpha \enspace \mathrm{ln}\enspace R\right)}\right]G\left({R}_{1},{\varphi }_{1},{R}_{2},{\varphi }_{2},\dots \right),\end{equation} \tag{ 7 }$$

$$\begin{align}\hfill \dot {\varphi }& =\omega +\frac{\varepsilon }{R}\left[\right.\mathrm{I}\mathrm{m}\left[{\text{e}}^{-\text{i}\left(\varphi +\alpha \enspace \mathrm{ln}\enspace R\right)}\enspace G\left({R}_{1},{\varphi }_{1},{R}_{2},{\varphi }_{2},\dots \right)\right]\hfill \\ \hfill & \quad -\alpha \mathrm{R}\mathrm{e}\left[{\text{e}}^{-\text{i}\left(\varphi +\alpha \enspace \mathrm{ln}\enspace R\right)}\enspace G\left({R}_{1},{\varphi }_{1},{R}_{2},{\varphi }_{2},\dots \right)\right]\left.\right].\hfill \end{align} \tag{ 8 }$$

Here we have explicitly used that the coupling terms depend on the amplitudes and the phases of other oscillators.

Now, we outline the exploited perturbation procedure, while in the next subsection we will elaborate on it for a particular example. The idea is to look for a dynamics, in which the amplitudes are enslaved by the phases. Namely, we assume that the amplitude R is a function of the phases:

$$\begin{equation}R=1+\varepsilon {r}^{\left(1\right)}\left(\varphi ,{\varphi }_{1},{\varphi }_{2},\dots \right)+{\varepsilon }^{2}{r}^{\left(2\right)}\left(\varphi ,{\varphi }_{1},{\varphi }_{2},\dots \right)+\cdots \enspace ,\end{equation} \tag{ 9 }$$

and that the dynamics of the phases is also represented as a power series in ɛ:

$$\begin{equation}\dot {\varphi }=\omega +\varepsilon {{\Psi}}^{\left(1\right)}\left(\varphi ,{\varphi }_{1},{\varphi }_{2},\dots \right)+{\varepsilon }^{2}{{\Psi}}^{\left(2\right)}\left(\varphi ,{\varphi }_{1},{\varphi }_{2},\dots \right)+\cdots \enspace .\end{equation} \tag{ 10 }$$

We substitute these expressions, together with the similar expressions for R_k , φ_k , into equations (7) and (8). Equating the terms in each power of ɛ, we obtain the unknown functions r⁽¹⁾, r⁽²⁾, Ψ⁽¹⁾, Ψ⁽²⁾, ....

We illustrate here the first few steps. In the equation for the phase, the first-order approximation, as discussed above, corresponds to taking into account only the leading term in equation (9), this yields

$$\begin{equation*}{{\Psi}}^{\left(1\right)}=\mathrm{I}\mathrm{m}\left[{\text{e}}^{-\text{i}\varphi }\enspace G\left(1,{\varphi }_{1},1,{\varphi }_{2},\dots \right)\right]-\alpha \enspace \mathrm{R}\mathrm{e}\left[{\text{e}}^{-\text{i}\varphi }\enspace G\left(1,{\varphi }_{1},1,{\varphi }_{2},\dots \right)\right].\end{equation*}$$

The equation for the amplitude in the first order is obtained by substituting equation (9) in equation (7):

$$\begin{equation*}\frac{\mathrm{d}{r}^{\left(1\right)}}{\mathrm{d}t}=-2{r}^{\left(1\right)}+\mathrm{R}\mathrm{e}\left[{\text{e}}^{-\text{i}\varphi }\enspace G\left(1,{\varphi }_{1},1,{\varphi }_{2},\dots \right)\right].\end{equation*}$$

We express the time derivative of r⁽¹⁾ via partial derivatives with respect to the phases, and insert the zero-order expressions for these derivatives:

$$\begin{align}\hfill \frac{\mathrm{d}{r}^{\left(1\right)}}{\mathrm{d}t}& =\frac{\partial {r}^{\left(1\right)}}{\partial \varphi }\dot {\varphi }+\frac{\partial {r}^{\left(1\right)}}{\partial {\varphi }_{1}}{\dot {\varphi }}_{1}+\frac{\partial {r}^{\left(1\right)}}{\partial {\varphi }_{2}}{\dot {\varphi }}_{2}+\cdots \hfill \\ \hfill & \approx \frac{\partial {r}^{\left(1\right)}}{\partial \varphi }\omega +\frac{\partial {r}^{\left(1\right)}}{\partial {\varphi }_{1}}{\omega }_{1}+\frac{\partial {r}^{\left(1\right)}}{\partial {\varphi }_{2}}{\omega }_{2}+\cdots \enspace .\hfill \end{align} \tag{ 11 }$$

Thus, the problem of determining first-order correction to the amplitude reduces to the partial differential equation

$$\begin{align}\hfill & 2{r}^{\left(1\right)}+\frac{\partial {r}^{\left(1\right)}}{\partial \varphi }\omega +\frac{\partial {r}^{\left(1\right)}}{\partial {\varphi }_{1}}{\omega }_{1}+\frac{\partial {r}^{\left(1\right)}}{\partial {\varphi }_{2}}{\omega }_{2}+\cdots =\mathrm{R}\mathrm{e}\left[{\text{e}}^{-\text{i}\varphi }\enspace G\left(1,{\varphi }_{1},1,{\varphi }_{2},\dots \right)\right]\hfill \\ \hfill & =\sum _{m,{m}_{1},{m}_{2},\dots }{g}_{m{m}_{1}{m}_{2}\dots }{\text{e}}^{\text{i}\left(m\varphi +{m}_{1}{\varphi }_{1}+{m}_{2}{\varphi }_{2}+\cdots \enspace \right)}.\hfill \end{align} \tag{ 12 }$$

Here, we used that the function G is 2π-periodic with respect to the phases and, hence, the r.h.s. can be written as a multiple Fourier series. Similarly, we expand the function r⁽¹⁾ as:

$$\begin{equation*}{r}^{\left(1\right)}=\sum _{m,{m}_{1},{m}_{2},\dots }{\rho }_{m{m}_{1}{m}_{2}\dots }{\text{e}}^{\text{i}\left(m\varphi +{m}_{1}{\varphi }_{1}+{m}_{2}{\varphi }_{2}+\cdots \enspace \right)},\end{equation*}$$

which finally yields an expression for the Fourier coefficients of r⁽¹⁾

$$\begin{equation}{\rho }_{m{m}_{1}{m}_{2}\dots }=\frac{{g}_{m{m}_{1}{m}_{2}\dots }}{2+\text{i}m\omega +\text{i}{m}_{1}{\omega }_{1}+\text{i}{m}_{2}{\omega }_{2}+\cdots }.\end{equation} \tag{ 13 }$$

These are the basic steps in the perturbation expansion. The expressions for ${r}^{\left(1\right)},{r}_{k}^{\left(1\right)}$, being substituted in equation (8) provide phase equations in order ∼ɛ². Equations for r⁽²⁾ are partial differential equations of type (12) with r.h.s. containing also Ψ⁽¹⁾, r⁽¹⁾, etc.

In this section, we exemplify the perturbative procedure for the derivation of phase dynamics equations in a higher order of the perturbation parameter ɛ with a particular configuration of three coupled Stuart–Landau oscillators.

2.4.1. Configuration of the network

We consider an array of three elements with the coupling structure 1 ↔ 2 ↔ 3 and write the equations in the form (4):

$$\begin{equation}\begin{aligned}\hfill {\dot {A}}_{1}& =\left(1+\mathrm{i}{\omega }_{1}\right){A}_{1}-\vert {A}_{1}{\vert }^{2}{A}_{1}-\mathrm{i}\alpha {A}_{1}\left(\vert {A}_{1}{\vert }^{2}-1\right)+\varepsilon {c}_{2,1}\enspace {\text{e}}^{\mathrm{i}{\beta }_{2,1}}{A}_{2},\hfill \\ \hfill {\dot {A}}_{2}& =\left(1+\mathrm{i}{\omega }_{2}\right){A}_{2}-\vert {A}_{2}{\vert }^{2}{A}_{2}-\mathrm{i}\alpha {A}_{2}\left(\vert {A}_{2}{\vert }^{2}-1\right)+\varepsilon \left({c}_{1,2}\enspace {\text{e}}^{\mathrm{i}{\beta }_{1,2}}{A}_{1}+{c}_{3,2}\enspace {\text{e}}^{\mathrm{i}{\beta }_{3,2}}{A}_{3}\right),\hfill \\ \hfill {\dot {A}}_{3}& =\left(1+\mathrm{i}{\omega }_{3}\right){A}_{3}-\vert {A}_{3}{\vert }^{2}{A}_{3}-\mathrm{i}\alpha {A}_{3}\left(\vert {A}_{3}{\vert }^{2}-1\right)+\varepsilon {c}_{2,3}\enspace {\text{e}}^{\mathrm{i}{\beta }_{2,3}}{A}_{2}.\hfill \end{aligned}\end{equation} \tag{ 14 }$$

Notice that the oscillators have different frequencies. Introducing the amplitudes and the angle variables according to A_k = R_k  exp[iθ_k ] and the phases φ_k = θ_k − α  ln  R_k , we obtain a system of equations in the form (7) and (8):

$$\begin{equation}\begin{aligned}\hfill {\dot {R}}_{1}& ={R}_{1}-{R}_{1}^{3}+\varepsilon {c}_{2,1}{R}_{2}\enspace \mathrm{cos}\left({\theta }_{2}-{\theta }_{1}+{\beta }_{2,1}\right),\hfill \\ \hfill {\dot {R}}_{2}& ={R}_{2}-{R}_{2}^{3}+\varepsilon {c}_{1,2}{R}_{1}\enspace \mathrm{cos}\left({\theta }_{1}-{\theta }_{2}+{\beta }_{1,2}\right)\hfill \\ \hfill & \quad +\varepsilon {c}_{3,2}{R}_{3}\enspace \mathrm{cos}\left({\theta }_{3}-{\theta }_{2}+{\beta }_{3,2}\right),\hfill \\ \hfill {\dot {R}}_{3}& ={R}_{3}-{R}_{3}^{3}+\varepsilon {c}_{2,3}{R}_{2}\enspace \mathrm{cos}\left({\theta }_{2}-{\theta }_{3}+{\beta }_{2,3}\right),\hfill \\ \hfill \dot {{\varphi }_{1}}& ={\omega }_{1}+\varepsilon {c}_{2,1}\frac{{R}_{2}}{{R}_{1}}\left[\mathrm{sin}\left({\theta }_{2}-{\theta }_{1}+{\beta }_{2,1}\right)-\alpha \enspace \mathrm{cos}\left({\theta }_{2}-{\theta }_{1}+{\beta }_{2,1}\right)\right],\hfill \\ \hfill \dot {{\varphi }_{2}}& ={\omega }_{2}+\varepsilon {c}_{1,2}\frac{{R}_{1}}{{R}_{2}}\left[\mathrm{sin}\left({\theta }_{1}-{\theta }_{2}+{\beta }_{1,2}\right)-\alpha \enspace \mathrm{cos}\left({\theta }_{1}-{\theta }_{2}+{\beta }_{1,2}\right)\right]\hfill \\ \hfill & \quad +\varepsilon {c}_{3,2}\frac{{R}_{3}}{{R}_{2}}\left[\mathrm{sin}\left({\theta }_{3}-{\theta }_{2}+{\beta }_{3,2}\right)-\alpha \enspace \mathrm{cos}\left({\theta }_{3}-{\theta }_{2}+{\beta }_{3,2}\right)\right],\hfill \\ \hfill \dot {{\varphi }_{3}}& ={\omega }_{3}+\varepsilon {c}_{2,3}\frac{{R}_{2}}{{R}_{3}}\left[\mathrm{sin}\left({\theta }_{2}-{\theta }_{3}+{\beta }_{2,3}\right)-\alpha \enspace \mathrm{cos}\left({\theta }_{2}-{\theta }_{3}+{\beta }_{2,3}\right)\right].\hfill \end{aligned}\end{equation} \tag{ 15 }$$

2.4.2. Small parameter expansion

Now we expand the amplitude deviations in powers of ɛ (below k = 1, 2, 3):

$$\begin{equation}{R}_{k}=1+\varepsilon {r}_{k}^{\left(1\right)}+{\varepsilon }^{2}{r}_{k}^{\left(2\right)}+{\varepsilon }^{3}{r}_{k}^{\left(3\right)}+\cdots \enspace .\end{equation} \tag{ 16 }$$

According to this, the angles θ_k can be represented as

$$\begin{equation}{\theta }_{k}={\varphi }_{k}+\alpha \left[\varepsilon {r}_{k}^{\left(1\right)}+{\varepsilon }^{2}{r}_{k}^{\left(2\right)}-0.5{\varepsilon }^{2}{\left({r}_{k}^{\left(1\right)}\right)}^{2}\right]+\cdots \enspace .\end{equation} \tag{ 17 }$$

Also the ratios of the amplitudes, entering (15), can be expressed as

$$\begin{equation}\frac{{R}_{m}}{{R}_{k}}=1+\varepsilon \left[{r}_{m}^{\left(1\right)}-{r}_{k}^{\left(1\right)}\right]+{\varepsilon }^{2}\left[{r}_{m}^{\left(2\right)}-{r}_{k}^{\left(2\right)}-{r}_{m}^{\left(1\right)}{r}_{k}^{\left(1\right)}+{\left({r}_{k}^{\left(1\right)}\right)}^{2}\right]+\cdots \enspace .\end{equation} \tag{ 18 }$$

Substituting these expansions in equation (15), we obtain the following expressions for the dynamics of the phases, up to the order ɛ³:

$$\begin{equation}\begin{aligned}\hfill {\dot {\varphi }}_{1}& ={\omega }_{1}+\varepsilon {c}_{2,1}\left[\mathrm{sin}\left({\varphi }_{2}-{\varphi }_{1}+{\beta }_{2,1}\right)-\alpha \enspace \mathrm{cos}\left({\varphi }_{2}-{\varphi }_{1}+{\beta }_{2,1}\right)\right]\hfill \\ \hfill & \quad +{\varepsilon }^{2}{c}_{2,1}\left(1+{\alpha }^{2}\right)\mathrm{sin}\left({\varphi }_{2}-{\varphi }_{1}+{\beta }_{2,1}\right)\left({r}_{2}^{\left(1\right)}-{r}_{1}^{\left(1\right)}\right)\hfill \\ \hfill & \quad +{\varepsilon }^{3}{c}_{2,1}\left(1+{\alpha }^{2}\right)\left[\left({r}_{2}^{\left(2\right)}-{r}_{1}^{\left(2\right)}-{r}_{2}^{\left(1\right)}{r}_{1}^{\left(1\right)}+{\left({r}_{1}^{\left(1\right)}\right)}^{2}\right)\mathrm{sin}\left({\varphi }_{2}-{\varphi }_{1}+{\beta }_{2,1}\right)\right.\hfill \\ \hfill & \quad \left.+\alpha \left(0.5{\left({r}_{1}^{\left(1\right)}\right)}^{2}+0.5{\left({r}_{2}^{\left(1\right)}\right)}^{2}-{r}_{2}^{\left(1\right)}{r}_{1}^{\left(1\right)}\right)\mathrm{cos}\left({\varphi }_{2}-{\varphi }_{1}+{\beta }_{2,1}\right)\right]+\cdots \enspace ,\hfill \\ \hfill {\dot {\varphi }}_{2}& ={\omega }_{2}+\varepsilon {c}_{1,2}\left[\mathrm{sin}\left({\varphi }_{1}-{\varphi }_{2}+{\beta }_{1,2}\right)-\alpha \enspace \mathrm{cos}\left({\varphi }_{1}-{\varphi }_{2}+{\beta }_{1,2}\right)\right]\hfill \\ \hfill & \quad +\varepsilon {c}_{3,2}\left[\mathrm{sin}\left({\varphi }_{3}-{\varphi }_{2}+{\beta }_{3,2}\right)-\alpha \enspace \mathrm{cos}\left({\varphi }_{3}-{\varphi }_{2}+{\beta }_{3,2}\right)\right]\hfill \\ \hfill & \quad +{\varepsilon }^{2}{c}_{1,2}\left(1+{\alpha }^{2}\right)\mathrm{sin}\left({\varphi }_{1}-{\varphi }_{2}+{\beta }_{1,2}\right)\left({r}_{1}^{\left(1\right)}-{r}_{2}^{\left(1\right)}\right)\hfill \\ \hfill & \quad +{\varepsilon }^{2}{c}_{3,2}\left(1+{\alpha }^{2}\right)\mathrm{sin}\left({\varphi }_{3}-{\varphi }_{2}+{\beta }_{3,2}\right)\left({r}_{3}^{\left(1\right)}-{r}_{2}^{\left(1\right)}\right)\hfill \\ \hfill & \quad +{\varepsilon }^{3}{c}_{1,2}\left(1+{\alpha }^{2}\right)\left[\left({r}_{1}^{\left(2\right)}-{r}_{2}^{\left(2\right)}-{r}_{1}^{\left(1\right)}{r}_{2}^{\left(1\right)}+{\left({r}_{2}^{\left(1\right)}\right)}^{2}\right)\mathrm{sin}\left({\varphi }_{1}-{\varphi }_{2}+{\beta }_{1,2}\right)\right.\hfill \\ \hfill & \quad \left.+\alpha \left(0.5{\left({r}_{2}^{\left(1\right)}\right)}^{2}+0.5{\left({r}_{1}^{\left(1\right)}\right)}^{2}-{r}_{1}^{\left(1\right)}{r}_{2}^{\left(1\right)}\right)\mathrm{cos}\left({\varphi }_{1}-{\varphi }_{2}+{\beta }_{1,2}\right)\right]\hfill \\ \hfill & \quad +{\varepsilon }^{3}{c}_{3,2}\left(1+{\alpha }^{2}\right)\left[\left({r}_{3}^{\left(2\right)}-{r}_{2}^{\left(2\right)}-{r}_{3}^{\left(1\right)}{r}_{2}^{\left(1\right)}+{\left({r}_{2}^{\left(1\right)}\right)}^{2}\right)\mathrm{sin}\left({\varphi }_{3}-{\varphi }_{2}+{\beta }_{3,2}\right)\right.\hfill \\ \hfill & \quad \left.+\alpha \left(0.5{\left({r}_{2}^{\left(1\right)}\right)}^{2}+0.5{\left({r}_{3}^{\left(1\right)}\right)}^{2}-{r}_{3}^{\left(1\right)}{r}_{2}^{\left(1\right)}\right)\mathrm{cos}\left({\varphi }_{3}-{\varphi }_{2}+{\beta }_{3,2}\right)\right]+\cdots \enspace ,\hfill \\ \hfill {\dot {\varphi }}_{3}& ={\omega }_{3}+\varepsilon {c}_{2,3}\left[\mathrm{sin}\left({\varphi }_{2}-{\varphi }_{3}+{\beta }_{2,3}\right)-\alpha \enspace \mathrm{cos}\left({\varphi }_{2}-{\varphi }_{3}+{\beta }_{2,3}\right)\right]\hfill \\ \hfill & \quad +{\varepsilon }^{2}{c}_{2,3}\left(1+{\alpha }^{2}\right)\mathrm{sin}\left({\varphi }_{2}-{\varphi }_{3}+{\beta }_{2,3}\right)\left({r}_{2}^{\left(1\right)}-{r}_{3}^{\left(1\right)}\right)\hfill \\ \hfill & \quad +{\varepsilon }^{3}{c}_{2,3}\left(1+{\alpha }^{2}\right)\left[\left({r}_{2}^{\left(2\right)}-{r}_{3}^{\left(2\right)}-{r}_{2}^{\left(1\right)}{r}_{3}^{\left(1\right)}+{\left({r}_{3}^{\left(1\right)}\right)}^{2}\right)\mathrm{sin}\left({\varphi }_{2}-{\varphi }_{3}+{\beta }_{2,3}\right)\right.\hfill \\ \hfill & \quad \left.+\alpha \left(0.5{\left({r}_{3}^{\left(1\right)}\right)}^{2}+0.5{\left({r}_{2}^{\left(1\right)}\right)}^{2}-{r}_{2}^{\left(1\right)}{r}_{3}^{\left(1\right)}\right)\mathrm{cos}\left({\varphi }_{2}-{\varphi }_{3}+{\beta }_{2,3}\right)\right]+\cdots \enspace .\hfill \end{aligned}\end{equation} \tag{ 19 }$$

Next, we have to evaluate corrections ${r}_{k}^{\left(1\right)},{r}_{k}^{\left(2\right)},\dots$. This is accomplished by substituting expressions (16) in the equations for the amplitudes in (15). Here the time derivatives are calculated according to the chain rule, as the corrections ${r}_{k}^{\left(1\right)},{r}_{k}^{\left(2\right)},\dots$ are assumed to be functions of the phases φ_k . For the sake of brevity, we present only the formulas for r₁:

$$\begin{equation}\begin{aligned}\hfill & {\omega }_{1}\frac{\partial {r}_{1}^{\left(1\right)}}{\partial {\varphi }_{1}}+{\omega }_{2}\frac{\partial {r}_{1}^{\left(1\right)}}{\partial {\varphi }_{2}}+{\omega }_{3}\frac{\partial {r}_{1}^{\left(1\right)}}{\partial {\varphi }_{3}}+2{r}_{1}^{\left(1\right)}={c}_{2,1}\enspace \mathrm{cos}\left({\varphi }_{2}-{\varphi }_{1}+{\beta }_{2,1}\right),\hfill \\ \hfill & {\omega }_{1}\frac{\partial {r}_{1}^{\left(2\right)}}{\partial {\varphi }_{1}}+{\omega }_{2}\frac{\partial {r}_{1}^{\left(2\right)}}{\partial {\varphi }_{2}}+{\omega }_{3}\frac{\partial {r}_{1}^{\left(2\right)}}{\partial {\varphi }_{3}}+2{r}_{1}^{\left(2\right)}=-3{\left({r}_{1}^{\left(1\right)}\right)}^{2}\hfill \\ \hfill & -\alpha {c}_{2,1}\left({r}_{2}^{\left(1\right)}-{r}_{1}^{\left(1\right)}\right)\mathrm{sin}\left({\varphi }_{2}-{\varphi }_{1}+{\beta }_{2,1}\right)+{c}_{2,1}{r}_{2}^{\left(1\right)}\enspace \mathrm{cos}\left({\varphi }_{2}-{\varphi }_{1}+{\beta }_{2,1}\right)\hfill \\ \hfill & +{c}_{2,1}\left[\mathrm{sin}\left({\varphi }_{2}-{\varphi }_{1}+{\beta }_{2,1}\right)-\alpha \enspace \mathrm{cos}\left({\varphi }_{2}-{\varphi }_{1}+{\beta }_{2,1}\right)\right]\frac{\partial {r}_{1}^{\left(1\right)}}{\partial {\varphi }_{1}}\hfill \\ \hfill & +\left[{c}_{1,2}\left[\mathrm{sin}\left({\varphi }_{1}-{\varphi }_{2}+{\beta }_{1,2}\right)-\alpha \enspace \mathrm{cos}\left({\varphi }_{1}-{\varphi }_{2}+{\beta }_{1,2}\right)\right]\right.\hfill \\ \hfill & \left.+{c}_{3,2}\left[\mathrm{sin}\left({\varphi }_{3}-{\varphi }_{2}+{\beta }_{3,2}\right)-\alpha \enspace \mathrm{cos}\left({\varphi }_{3}-{\varphi }_{2}+{\beta }_{3,2}\right)\right]\right]\frac{\partial {r}_{1}^{\left(1\right)}}{\partial {\varphi }_{2}}\hfill \\ \hfill & +{c}_{2,3}\left[\mathrm{sin}\left({\varphi }_{2}-{\varphi }_{3}+{\beta }_{2,3}\right)-\alpha \enspace \mathrm{cos}\left({\varphi }_{2}-{\varphi }_{3}+{\beta }_{2,3}\right)\right]\frac{\partial {r}_{1}^{\left(1\right)}}{\partial {\varphi }_{3}}.\hfill \end{aligned}\end{equation} \tag{ 20 }$$

The partial differential equations for ${r}_{k}^{\left(1\right)},{r}_{k}^{\left(2\right)},\dots$ are straightforwardly solved in the Fourier representation, because the variations of the amplitudes are 2π-periodic functions of the phases, as outlined in discussion of equation (13) above. As a result, we obtain in the 1st order in ɛ:

$$\begin{equation}\begin{aligned}\hfill {r}_{1}^{\left(1\right)}& =\frac{2{c}_{2,1}}{4+{\left({\omega }_{2}-{\omega }_{1}\right)}^{2}}\enspace \mathrm{cos}\left({\varphi }_{2}-{\varphi }_{1}+{\beta }_{2,1}\right)+\frac{\left({\omega }_{2}-{\omega }_{1}\right){c}_{2,1}}{4+{\left({\omega }_{2}-{\omega }_{1}\right)}^{2}}\enspace \mathrm{sin}\left({\varphi }_{2}-{\varphi }_{1}+{\beta }_{2,1}\right),\hfill \\ \hfill {r}_{2}^{\left(1\right)}& =\frac{2{c}_{1,2}}{4+{\left({\omega }_{1}-{\omega }_{2}\right)}^{2}}\enspace \mathrm{cos}\left({\varphi }_{1}-{\varphi }_{2}+{\beta }_{1,2}\right)+\frac{\left({\omega }_{1}-{\omega }_{2}\right){c}_{1,2}}{4+{\left({\omega }_{1}-{\omega }_{2}\right)}^{2}}\enspace \mathrm{sin}\left({\varphi }_{1}-{\varphi }_{2}+{\beta }_{1,2}\right)\hfill \\ \hfill & \quad +\frac{2{c}_{3,2}}{4+{\left({\omega }_{3}-{\omega }_{2}\right)}^{2}}\enspace \mathrm{cos}\left({\varphi }_{3}-{\varphi }_{2}+{\beta }_{3,2}\right)+\frac{\left({\omega }_{3}-{\omega }_{2}\right){c}_{3,2}}{4+{\left({\omega }_{3}-{\omega }_{2}\right)}^{2}}\enspace \mathrm{sin}\left({\varphi }_{3}-{\varphi }_{2}+{\beta }_{3,2}\right),\hfill \\ \hfill {r}_{3}^{\left(1\right)}& =\frac{2{c}_{2,3}}{4+{\left({\omega }_{2}-{\omega }_{3}\right)}^{2}}\enspace \mathrm{cos}\left({\varphi }_{2}-{\varphi }_{3}+{\beta }_{2,3}\right)+\frac{\left({\omega }_{2}-{\omega }_{3}\right){c}_{2,3}}{4+{\left({\omega }_{2}-{\omega }_{3}\right)}^{2}}\enspace \mathrm{sin}\left({\varphi }_{2}-{\varphi }_{3}+{\beta }_{2,3}\right).\hfill \end{aligned}\end{equation} \tag{ 21 }$$

Substitution of these expression in equation (19) completes the second-order phase reduction and yields closed equations:

$$\begin{equation}\begin{aligned}\hfill {\dot {\varphi }}_{1}& ={\omega }_{1}+\varepsilon {c}_{2,1}\left[\mathrm{sin}\left({\varphi }_{2}-{\varphi }_{1}+{\beta }_{2,1}\right)-\alpha \enspace \mathrm{cos}\left({\varphi }_{2}-{\varphi }_{1}+{\beta }_{2,1}\right)\right]\hfill \\ \hfill & \quad +{\varepsilon }^{2}\left[{a}_{1;0}^{\left(2\right)}+{a}_{1;-2,2,0}^{\left(2\right)}\enspace \mathrm{cos}\left(2{\varphi }_{2}-2{\varphi }_{1}\right)+{b}_{1;-2,2,0}^{\left(2\right)}\enspace \mathrm{sin}\left(2{\varphi }_{2}-2{\varphi }_{1}\right)\right.\hfill \\ \hfill & \quad \left.+{a}_{1;-1,2,-1}^{\left(2\right)}\enspace \mathrm{cos}\left(2{\varphi }_{2}-{\varphi }_{1}-{\varphi }_{3}\right)+{b}_{1;-1,2,-1}^{\left(2\right)}\enspace \mathrm{sin}\left(2{\varphi }_{2}-{\varphi }_{1}-{\varphi }_{3}\right)\right.\hfill \\ \hfill & \quad \left.+{a}_{1;-1,0,1}^{\left(2\right)}\enspace \mathrm{cos}\left({\varphi }_{3}-{\varphi }_{1}\right)+{b}_{1;-1,0,1}^{\left(2\right)}\enspace \mathrm{sin}\left({\varphi }_{3}-{\varphi }_{1}\right)\right],\hfill \end{aligned}\end{equation} \tag{ 22 }$$

$$\begin{equation}\begin{aligned}\hfill {\dot {\varphi }}_{2}& ={\omega }_{2}+\varepsilon {c}_{1,2}\left[\mathrm{sin}\left({\varphi }_{1}-{\varphi }_{2}+{\beta }_{1,2}\right)-\alpha \enspace \mathrm{cos}\left({\varphi }_{1}-{\varphi }_{2}+{\beta }_{1,2}\right)\right]\hfill \\ \hfill & \quad +\varepsilon {c}_{3,2}\left[\mathrm{sin}\left({\varphi }_{3}-{\varphi }_{2}+{\beta }_{3,2}\right)-\alpha \enspace \mathrm{cos}\left({\varphi }_{3}-{\varphi }_{2}+{\beta }_{3,2}\right)\right]\hfill \\ \hfill & \quad +{\varepsilon }^{2}\left[{a}_{2;0}^{\left(2\right)}+{a}_{2;2,-2,0}^{\left(2\right)}\enspace \mathrm{cos}\left(2{\varphi }_{1}-2{\varphi }_{2}\right)+{b}_{2;2,-2,0}^{\left(2\right)}\enspace \mathrm{sin}\left(2{\varphi }_{1}-2{\varphi }_{2}\right)\right.\hfill \\ \hfill & \quad \left.+{a}_{2;0,-2,2}^{\left(2\right)}\enspace \mathrm{cos}\left(2{\varphi }_{3}-2{\varphi }_{2}\right)+{b}_{2;0,-2,2}^{\left(2\right)}\enspace \mathrm{sin}\left(2{\varphi }_{3}-2{\varphi }_{2}\right)\right.\hfill \\ \hfill & \quad \left.+{a}_{2;-1,2,-1}^{\left(2\right)}\enspace \mathrm{cos}\left(2{\varphi }_{2}-{\varphi }_{1}-{\varphi }_{3}\right)+{b}_{2;-1,2,-1}^{\left(2\right)}\enspace \mathrm{sin}\left(2{\varphi }_{2}-{\varphi }_{1}-{\varphi }_{3}\right)\right.\hfill \\ \hfill & \quad \left.+{a}_{2;1,0,-1}^{\left(2\right)}\enspace \mathrm{cos}\left({\varphi }_{1}-{\varphi }_{3}\right)+{b}_{2;1,0,-1}^{\left(2\right)}\enspace \mathrm{sin}\left({\varphi }_{1}-{\varphi }_{3}\right)\right],\hfill \end{aligned}\end{equation} \tag{ 23 }$$

$$\begin{equation}\begin{aligned}\hfill {\dot {\varphi }}_{3}& ={\omega }_{3}+\varepsilon {c}_{2,3}\left[\mathrm{sin}\left({\varphi }_{2}-{\varphi }_{3}+{\beta }_{2,3}\right)-\alpha \enspace \mathrm{cos}\left({\varphi }_{2}-{\varphi }_{3}+{\beta }_{2,3}\right)\right]\hfill \\ \hfill & \quad +{\varepsilon }^{2}\left[{a}_{3;0}^{\left(2\right)}+{a}_{3;0,2,-2}^{\left(2\right)}\enspace \mathrm{cos}\left(2{\varphi }_{2}-2{\varphi }_{3}\right)+{b}_{3;0,2,-2}^{\left(2\right)}\enspace \mathrm{sin}\left(2{\varphi }_{2}-2{\varphi }_{3}\right)\right.\hfill \\ \hfill & \quad \left.+{a}_{3;-1,2,-1}^{\left(2\right)}\enspace \mathrm{cos}\left(2{\varphi }_{2}-{\varphi }_{3}-{\varphi }_{1}\right)+{b}_{3;-1,2,-1}^{\left(2\right)}\enspace \mathrm{sin}\left(2{\varphi }_{2}-{\varphi }_{3}-{\varphi }_{1}\right)\right.\hfill \\ \hfill & \quad \left.+{a}_{3;1,0,-1}^{\left(2\right)}\enspace \mathrm{cos}\left({\varphi }_{1}-{\varphi }_{3}\right)+{b}_{3;1,0,-1}^{\left(2\right)}\enspace \mathrm{sin}\left({\varphi }_{1}-{\varphi }_{3}\right)\right].\hfill \end{aligned}\end{equation} \tag{ 24 }$$

The coefficients of the second-order coupling terms, denoted in equations (22)–(24) by ${a}_{k;\boldsymbol{l}}^{\left(2\right)},{b}_{k;\boldsymbol{l}}^{\left(2\right)}$, are listed in tables A1–A3. Here the three-component vector l = (l₁, l₂, l₃) is used to signify the term with the combination of the phases l₁ φ₁ + l₂ φ₂ + l₃ φ₃. Furthermore, in table B1 we list the terms (without coupling coefficients) appearing in orders ɛ³, ɛ⁴.

We now shortly discuss the physical meaning of the terms that appear in the leading and the higher orders in ɛ.

• (a)  
  In the first order in ɛ, we obtain the standard phase-reduced description of the coupled oscillators, as described, e.g., in [15, 16]. If the focus is on first-order in ɛ effects, no further treatment is needed.
• (b)  
  In the second-order approximation, there are no correction terms to the first-order couplings. These corrections appear in the third order (and, presumably, in all odd orders).
• (c)  
  There are terms, which can be roughly described as 'squares' of the basic coupling terms; for the dynamics of the 1st oscillator φ₁ these are constant terms and terms containing the second harmonics of the phase difference, e.g. ∼sin(2φ₂ − 2φ₁). These high-order terms do not arise from the interaction within the whole network; they appear already in a system of two coupled oscillators.
• (d)  
  Terms containing combinations of all three phases, e.g., ∼sin(2φ₂ − φ₁ − φ₃), mean that an effective hypernetwork with non-pairwise coupling appears already in the second-order reduction (cf. studies of synchronization on hypernetworks of phase oscillators [41, 42]).
• (e)  
  Terms containing phase differences for not directly coupled oscillators, (e.g., the term ∼sin(φ₃ − φ₁) on the r.h.s. of the equation for ${\dot {\varphi }}_{1}$) mean that connections in terms of the phase dynamics do not coincide with the 'structural' connections in the original formulation. This phenomenon is well-known in the context of inverse problems of phase dynamics reconstruction from data [39]. In this reference, the second-order terms in the phase coupling reconstruction of small networks, which do not exist as 'structural' couplings, have been estimated numerically. The theory above provides an analytic justification of these findings.
• (f)  
  While the first-order coupling terms are frequency-independent, the second-order terms depend explicitly on frequency differences. In a more general setup (cf. a system of coupled van der Pol equations treated below) we expect that coupling will depend on the frequencies themselves.

The first step in numerical analysis is the determination of phases φ_i and their derivatives ${\dot {\varphi }}_{i}$ for all elements of the analyzed network, i = 1, ..., N. For this purpose, we extend the technique suggested in [43], where it was described how to obtain phases from an arbitrary trajectory.

Consider a particular oscillator within a network described by equation (1). Omitting the index of this unit for simplicity of presentation, we write

$$\begin{equation}\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t}=\mathbf{f}\left(\mathbf{y}\right)+\varepsilon \mathbf{G},\end{equation} \tag{ 25 }$$

where the last term describes the coupling to all other units. As a preparatory step we compute the autonomous period T. It is, for ɛ = 0 we take an arbitrary point on the limit cycle Y, and assign to it φ = 0; the return time to this point is exactly T. Now, for any other point on the limit cycle we can compute the time τ(Y) required to reach the zero point where φ = 0 or, equivalently, φ = 2π. Since the true phase grows linearly in time and gains 2π with one revolution around the cycle, its value can be obtained as $\varphi \left(\mathbf{Y}\right)=2\pi \frac{T-\tau \left(\mathbf{Y}\right)}{T}$.

The next step is to compute φ(t) and $\dot {\varphi }\left(t\right)$ for an arbitrary trajectory, on which at some time t the oscillator has the state u = y(t), and the time derivative of this state is $\mathbf{v}=\dot {\mathbf{y}}\left(t\right)$. To this end, we introduce an autonomous copy of the investigated unit:

$$\begin{equation}\frac{\mathrm{d}\mathbf{w}}{\mathrm{d}t}=\mathbf{f}\left(\mathbf{w}\right),\end{equation} \tag{ 26 }$$

and let this auxiliary system evolve from initial conditions w(0) = u, for the time interval nT, where the integer n shall be large enough to ensure that the trajectory attracts to the limit cycle. (Practically, we stop the evolution when ${\Vert}\mathbf{w}\left(\left(n-1\right)T\right)-\mathbf{w}\left(nT\right){\Vert}$ is smaller than a given tolerance.) Since the time interval of the evolution is a multiple of the period, the initial point w(0) = u and the end point $\mathbf{w}\left(nT\right)=\bar{\mathbf{w}}$ have same value of the phase. The point $\bar{\mathbf{w}}$ is on the limit cycle, and, hence, its phase φ can be easily computed as described above, and $\varphi \left(\mathbf{u}\right)=\varphi \left(\mathbf{w}\left(0\right)\right)=\varphi \left(\bar{\mathbf{w}}\right)=2\pi \frac{T-\tau \left(\bar{\mathbf{w}}\right)}{T}$ is exactly the desired phase of the oscillator at the state u.

In order to obtain the phase derivative we also have to follow the autonomous evolution (26) towards the limit cycle of the initial condition u + vdt. (Practically, it can be performed simultaneously with the evolution of the point u.) Since dt is (infinitely) small, this can be done by tracing the linear evolution of vdt to $\bar{\mathbf{v}}\mathrm{d}t$ within time interval nT. The law of this linear evolution is given by the Jacobian of the original equation (26). Thus, two states u and u + vdt of the coupled system (25) map to two points $\bar{\mathbf{w}}$ and $\bar{\mathbf{w}}+\bar{\mathbf{v}}\mathrm{d}t$ on the limit cycle of the autonomous system (26). These points are characterized by the phases φ and φ + dφ, respectively. On the other hand, let us consider the evolution of the autonomous system, from the point $\bar{\mathbf{w}}$ to $\bar{\mathbf{w}}+\bar{\mathbf{v}}\mathrm{d}t$. This evolution occurs along the limit cycle of the system (26) within time interval $\bar{\mathrm{d}t}$. (Notice that generally $\bar{\mathrm{d}t}\ne \mathrm{d}t$). The evolution is governed by the flow on the cycle, i.e. $\bar{\mathbf{w}}+\bar{\mathbf{v}}\mathrm{d}t=\bar{\mathbf{w}}+\mathbf{f}\left(\bar{\mathbf{w}}\right)\bar{\mathrm{d}t}$, what yields $\bar{\mathrm{d}t}=\mathrm{d}t\left(\bar{\mathbf{v}}\cdot \mathbf{f}\left(\bar{\mathbf{w}}\right)\right)/{\Vert}\mathbf{f}\left(\bar{\mathbf{w}}\right){{\Vert}}^{2}$. Accordingly, phase growth is determined by the natural frequency ω = 2π/T, i.e. $\mathrm{d}\varphi =\omega \bar{\mathrm{d}t}$, which finally yields

$$\begin{equation*}\frac{\mathrm{d}\varphi }{\mathrm{d}t}=\omega \frac{\bar{\mathbf{v}}\cdot \mathbf{f}\left(\bar{\mathbf{w}}\right)}{{\Vert}\mathbf{f}\left(\bar{\mathbf{w}}\right){{\Vert}}^{2}}.\end{equation*}$$

Thus, with the presented algorithm we can compute phases and their derivatives as time series with an arbitrary time step and of sufficient length. Certainly, this can be done for all elements of the network.

Now, we discuss how the phase dynamics equations of a network can be constructed numerically from given phases and their derivatives. To explain the procedure, it is convenient to denote the r.h.s. of equation (10) as Q_k (φ₁, φ₂, ..., φ_N ), where the subscript k = 1, 2, ..., N labels the oscillator. Q_k are commonly referred to as the coupling functions. Since each Q_k is 2π-periodic with respect to all its arguments, we can write it as a multiple Fourier series

$$\begin{equation}{\dot {\varphi }}_{k}={a}_{k;\boldsymbol{0}}+\sum _{\boldsymbol{l}\ne \mathbf{0}}\left[{a}_{k;\boldsymbol{l}}\enspace \mathrm{cos}\left(\boldsymbol{\varphi }\cdot \boldsymbol{l}\right)+{b}_{k;\boldsymbol{l}}\enspace \mathrm{sin}\left(\boldsymbol{\varphi }\cdot \boldsymbol{l}\right)\right],\end{equation} \tag{ 27 }$$

where l = (l₁, l₂, ..., l_N ) and φ = (φ₁, φ₂, ..., φ_N ) are N-dimensional vectors of integer indices and phases, respectively. $\boldsymbol{\varphi }\cdot \boldsymbol{l}={\sum }_{j=1}^{N}{l}_{j}{\varphi }_{j}$ denotes the scalar product between the vectors of phases and the mode indices. Notice the relation between the Fourier coefficients a_{k; l} , b_{k; l} , and coefficients in equations (22)–(24):

$$\begin{equation*}{a}_{k;\boldsymbol{l}}={a}_{k;\boldsymbol{l}}^{\left(0\right)}+\varepsilon {a}_{k;\boldsymbol{l}}^{\left(1\right)}+{\varepsilon }^{2}{a}_{k;\boldsymbol{l}}^{\left(2\right)}\dots , {b}_{k;\boldsymbol{l}}=\varepsilon {b}_{k;\boldsymbol{l}}^{\left(1\right)}+{\varepsilon }^{2}{b}_{k;\boldsymbol{l}}^{\left(2\right)}\dots .\end{equation*}$$

Thus, we have time series of the phases φ_1,n , φ_2,n , ..., φ_N,n and their derivatives ${\dot {\varphi }}_{1,n},{\dot {\varphi }}_{2,n},\dots ,{\dot {\varphi }}_{N,n}$, where n = 1, 2, ..., L is the point index, equation (27) becomes a system of L linear equations for the unknown coupling coefficients a, b. It is natural to approximate Q_k by a finite Fourier series, preserving only M ≪ L terms in the sum in equation (27). Practically, we vary the mode indices in the range |l_j | ⩽ m, which leaves M = 2(m + 1)(2m(2m + 1) + 1) − 1 unknown coefficients. We apply the least square fit method to find the Fourier modes in the coupling from the over-determined linear system. Practically this is accomplished via the singular value decomposition (SVD) as described in [44].

We denote the Fourier coefficients of the truncated series obtained from this procedure by capitals A_{k; l} , B_{k; l} . Thus, as a result of the numerical evaluation we obtain an approximation of the phase dynamics as

$$\begin{equation}{\dot {\varphi }}_{k}={A}_{k;\boldsymbol{0}}+\sum _{\boldsymbol{l}\ne \boldsymbol{0},\vert {l}_{j}\vert {\leqslant}m}\left[{A}_{k;\boldsymbol{l}}\enspace \mathrm{cos}\left(\boldsymbol{\varphi }\cdot \boldsymbol{l}\right)+{B}_{k;\boldsymbol{l}}\enspace \mathrm{sin}\left(\boldsymbol{\varphi }\cdot \boldsymbol{l}\right)\right].\end{equation} \tag{ 28 }$$

The sum in equation (27) goes over the combination of mode indices such that either l or − l is counted.

The success of this numerical approach depends on how how interdependent the series φ_k are. Here we use two different protocols.

Asynchronous case. Suppose that the network does not synchronize. It means that the trajectories of the system (1) are dense on the N-dimensional torus. Then one just has to calculate one long trajectory of the system (1) starting from some initial conditions and compute phases and instantaneous frequencies for each point of the solution, as described above. If the data series is sufficiently long, the computed trajectory covers the surface of the torus spanned by φ, φ₁, φ₂, ..., φ_N . Hence, we obtain enough information to recover the function Q of N variables and to determine the Fourier coefficients.

Synchronous case. If the network synchronizes, the trajectory on the torus collapses to a closed line and equation (27) for the Fourier coefficients become dependent and cannot be solved. However, there is a way to obtain enough data to solve the problem even in this case. For this goal one considers not one long trajectory, but an ensemble of short transients (cf. [45, 46]). Namely, one starts numerical integration with some asynchronous initial conditions and follows the trajectory unless it is on the invariant torus. The corresponding transient time should be larger than the amplitude relaxation time, but smaller than the synchronization time of the phases. Then the procedure is repeated with new initial conditions. Thus, instead of one long record one collects many dynamical states on the torus, unless the sufficient number of points is obtained. Certainly, this protocol can be use for the asynchronous case as well.

We exemplify these two protocols in the next section, but before proceeding with the examples, we have to discuss an important issue. As we mentioned, the least squares optimization requires that data points fill the surface of an N-dimensional torus. Thus, we face the curse of dimensionality: the data requirement grows fast with the network size N and becomes hardly feasible already for N > 3. Some information about the network, e.g., strength of directed links, can, however, be revealed for N > 3 as well. A possible approach is to perform the triplet-based analysis [40, 47–49].

Here, we perform the numerical analysis for the system (14) with the goal to verify main results in section 2.4 as well as the numerical approach. We choose two sets of parameters that correspond to asynchronous and synchronous dynamics, respectively, and employ the corresponding protocols.

The parameter values are: α = 0.1, while ω₁, ω₂ and ω₃ are varied to change the synchronization behaviour. Uncoupled oscillators with these parameters have a limit cycle with R = 1 and initial conditions were chosen on it so that the relaxation time is significantly reduced. The coupling parameters are all equal c_j,k = 1 for (j, k) ∈ {(1, 2), (2, 1), (2, 3), (3, 2)} and the phase lags are β_1,2 = 0.32, β_2,1 = 0.44, β_2,3 = 0.43 and β_3,2 = 0.18. We consider two sets of frequencies: case I with ${\omega }_{1}=-\sqrt{5}/2$, ${\omega }_{2}=\left(\sqrt{2}-1\right)/10$ and ω₃ = 0.8; here the dynamics is asynchronous quasiperiodic; while for case II with ω₁ = −0.055, ω₂ = 0 and ω₃ = 0.33 the dynamics is synchronous with a long transient. In both cases, we generated the set of points as follows: we started the dynamics of system (14) with random phases and amplitudes equal to one. Then, after an initial transient time Δt = 20 (which has been chosen to ensure that the relaxation of the amplitude to the invariant torus is over, but locking of the phases still does not occur), the values of the phases and their velocities were stored. Altogether we constructed a set of L = 10⁶ data points.

Next, we computed the coefficients of the truncated Fourier series, see equation (28), for different values of ɛ, using the SVD approach as described above. The system of coupled Stuart–Landau oscillators is invariant with respect to a phase shift φ_k → φ_k + ϕ which means that only modes that fulfill the condition

$$\begin{equation}\sum _{j=1}^{3}{l}_{j}=0\end{equation} \tag{ 29 }$$

can exist. To incorporate this into the analysis, we exploited two approaches:

• (a)  
  only modes satisfying (29) are taken into account, all other modes are set to zero; these modes constitute a small subset of all possible modes and therefore we determined the Fourier coefficients up to harmonics |l_j | ⩽ m = 8.
• (b)  
  all modes with |l_j | ⩽ m = 4 are determined.

Here we present the results for the case (a), while the results for (b) are presented in appendix C.

First of all, we compare the theoretical findings presented in equations (22)–(24) with numerical results. To this end, we show in figures 1 and 2 the differences |a_{k; l} − A_{k; l} |, |b_{k; l} − B_{k; l} |, where k = 1, 2, 3, for the theoretically known modes (see equations (22)–(24)), for the asynchronous and synchronous configurations, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We see that for weak coupling the difference is on the level of numerical precision; it becomes of the order of one only for such strong coupling as ɛ = 0.4. Next, we see that the difference for the first-order terms grows proportionally to ɛ³, in correspondence with our theoretical conclusion that there are no second-order correction to the first-order terms. Thus, numerical results exhibit a good correspondence with the theory.

Figures 3 and 4 present an overview of all Fourier coefficients for which we do not have theoretical values. Namely, we show here all coefficients except for those entering equations (22)–(24) and analyzed in figures 1 and 2. The results are in full agreement with the power-series representation of the coupling terms. Indeed, we see four groups of coefficients that scale as ɛ³, ɛ⁴, ɛ⁵ and ɛ⁶, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, we mention that the overall precision of the numerical procedure can be estimated by computing the rest term ξ_k of the Fourier series representation in equation (28). The results presented in figure 5 show that the rest term grows with the coupling strength as ɛ⁹. This is an indication that all the terms of orders from 0 to 8 are within the set of chosen Fourier modes.

Zoom In Zoom Out Reset image size

Up to this point, our analysis was restricted to the case of Stuart–Landau oscillators where expressions for phases and their derivatives are known. In this section we present a purely numerical evaluation of high-order coupling terms for a network of three non-identical van der Pol oscillators:

$$\begin{equation}\begin{aligned}\hfill {\ddot {x}}_{1}& -\mu \left(1-{x}_{1}^{2}\right){\dot {x}}_{1}+{\omega }_{1}^{2}{x}_{1}=\varepsilon {x}_{2}\hfill \\ \hfill {\ddot {x}}_{2}& -\mu \left(1-{x}_{2}^{2}\right){\dot {x}}_{2}+{\omega }_{2}^{2}{x}_{2}=\varepsilon \left({x}_{1}+{x}_{3}\right)\hfill \\ \hfill {\ddot {x}}_{3}& -\mu \left(1-{x}_{3}^{2}\right){\dot {x}}_{3}+{\omega }_{3}^{2}{x}_{3}=\varepsilon {x}_{2}.\hfill \end{aligned}\end{equation} \tag{ 30 }$$

We fixed parameters ω₁ = 1, ω₂ = 1.324 715 957, ${\omega }_{3}={\omega }_{2}^{2}$ (here ω₂ is the spiral mean, the root of cubic equation ${\omega }_{2}^{3}-{\omega }_{2}-1=0$). We fixed μ = 1 and varied the coupling constant ɛ in the range 0.001 ⩽ ɛ ⩽ 0.3, in this range the dynamics of the network is asynchronous. From a trajectory of the system (30) we obtained time series ${\varphi }_{1,2,3},{\dot {\varphi }}_{1,2,3}$ of length L = 10⁶ points each. We used this set to estimate the coefficients of the truncated Fourier series in equation (28) for m = 4. Thus, 729 coupling constants A_{k; l} , B_{k; l} were calculated for each ɛ. Below we restrict our attention only to the strengths of the coupling and ignore the relative phase, so we look on the properties of 364 coupling constants ${H}_{k,\boldsymbol{l}}={\left({A}_{k;\boldsymbol{l}}^{2}+{B}_{k;\boldsymbol{l}}^{2}\right)}^{1/2}$. Together with the free term A_k;0 this constitutes a set of 365 numerically obtained coefficients for each oscillator and for each coupling constant ɛ.

For presentation of the results we perform a preliminary sorting. According to the theory, we expect to obtain terms with leading dependencies ∼ɛ^q , with q = 0, 1, 2, .... Therefore, for each set of indices l and each oscillator, we tried to approximate the coefficients H_{k, l} (ɛ) by the function ∼ɛ^q , and if the fitted value q was close to an integer and the reliability of the fit was high, we attributed the corresponding power to the coupling term. Additionally, for ɛ⁰, ɛ, ɛ² we used the five small values of ɛ = 0.001, 0.002, 0.005, 0.01, 0.02, while for powers ɛ³ and ɛ⁴ we used larger coupling constants ɛ = 0.04, 0.06, 0.08, 0.1, 0.15. We did not look for powers 5 and larger.

Figures 6 and 7 illustrate these findings. It is instructive to look which coupling modes appear in each power of ɛ. For modes up to ∼ɛ² we summarize this in table 1. The modes A_k;0 , describing corrections to the natural frequency are not listed. Notice that in each cell of the table the modes are ordered according to their amplitudes. Below we summarize the properties of the coupling modes of the van der Pol oscillators (Eqs. (30)) :

• (a)  
  Looking at the modes that appear in the first order in ɛ, we notice that the terms with phase differences and the terms with phase sums have nearly the same amplitude. This means that the coupling terms ∼ɛ have nearly the Winfree form: they are products of two functions of the oscillator phase and of the driving oscillator phase, in full agreement with the first-order theory. Recall that it is different for the Stuart–Landau oscillators, where only terms with phase differences appear due to the condition (29).
• (b)  
  Because the variable x of the van der Pol equation possesses odd harmonics of the phase (and the same is true for the phase response curve), terms with the third harmonics appear already in the 1st order in ɛ.
• (c)  
  Inspection of terms ∼ɛ in table 1 reveals some terms that are not expected in the first-order analysis, we show them in italic font. Data in figure 6 show that the amplitudes of these terms are extremely small, comparable to the errors in terms reconstruction. We conclude that these terms are probably spurious and just occasionally possess a scaling ∼ɛ, what is not surprising due to the fact that the total number of terms to be found is large.
• (d)  
  With the size of the time series explored, we could not reliably detect coupling terms appearing in order ɛ⁵. However, as the figures show, the terms ∼4 can be detected with good confidence.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size