A parametrisation method for high-order phase reduction in coupled oscillator networks

Many systems in science and engineering consist of coupled periodic processes. Examples vary from the motion of the planets, to the synchronous flashing of fireflies [5], and from the activity of neurons in the brain [11], to power grids and electronic circuits. The functioning and malfunctioning of these coupled systems is often determined by a form of collective behaviour of its constituents, perhaps most notably their synchronisation [1, 26]. For example, synchronisation of neurons plays a critical role in cognitive processes [23, 24].

In this paper, we consider the situation where the coupling between the periodic processes is weak, a case that is amenable to rigorous mathematical analysis. Specifically, we assume that the evolution of the processes can be modelled by a system of differential equations of the form

$\displaystyle\dot{x}_{j}=F_{j}(x_{j})+\varepsilon\,G_{j}(x_{1},\ldots,x_{m})\ \mbox{for}\ x_{j}\in\mathbb{R}^{M_{j}}\ \mbox{and}\ j=1,\ldots,m\,.$

(1.1)

The vector fields $F_{j}:\mathbb{R}^{M_{j}}\to\mathbb{R}^{M_{j}}$ in (1.1) determine the dynamics of the uncoupled oscillators: we assume that each $F_{j}$ possesses a hyperbolic $T_{j}$-periodic orbit $X_{j}(t)$. In the uncoupled limit—when $\varepsilon=0$—equations (1.1) thus admit a normally hyperbolic periodic or quasi-periodic invariant torus $\mathbb{T}_{0}\subset\mathbb{R}^{M}$ (where $M:=M_{1}+\ldots+M_{m}$), consisting of the product of these periodic orbits. The functions $G_{j}$ in (1.1) model the interaction between the oscillators, for example through a (hyper-)network. The interaction strength $0\leq\varepsilon\ll 1$ is assumed small, so that the unperturbed torus $\mathbb{T}_{0}$ persists as an invariant manifold $\mathbb{T}_{\varepsilon}$ for (1.1), depending smoothly on $\varepsilon$, as is guaranteed by Fénichel’s theorem [8, 30].

The process of finding the equations of motion that govern the dynamics on the persisting torus $\mathbb{T}_{\varepsilon}$ is usually referred to as phase reduction [20, 21, 25]. Phase reduction has proved a powerful tool in the study of the synchronisation of coupled oscillators, especially because it often realises a considerable reduction of the dimension—and hence complexity—of the system. Various methods of phase reduction have been introduced over the past decades, the most well-known appearing perhaps in the work on chemical oscillations of Kuramoto [17]. We refer to [25] for an extensive overview of established phase reduction techniques, and refrain from providing an overview of these methods here.

Most existing phase reduction methods provide a first-order approximation of the dynamics on the persisting invariant torus in terms of the small coupling parameter. However, there are various instances where such a first-order approximation is insufficient, see [4, 16, 18, 27, 22], in particular when the first-order reduced dynamics is structurally unstable. For instance, it was observed in [16] that “remote synchronisation” [3] cannot be analysed with first-order methods. More accurate “high-order phase reduction” techniques (that go beyond the first-order approximation) have only been introduced very recently [4, 10, 18]. They have already been applied successfully, for example to predict remote synchronisation [16]. However, to the best of our knowledge, mathematically rigorous high-order phase reduction methods have only been derived in the special case that the unperturbed oscillators are either Stuart-Landau oscillators [10, 18] or deformations thereof [2, 4]. In that setting, phase reduction can be performed by computing an expansion of the phase-amplitude relation that defines the invariant torus. However, this procedure does not generalise to arbitrary systems of the form (1.1).

This paper presents a novel method for high-order phase reduction, that applies to general coupled oscillator systems of the form (1.1). Our method works by computing an expansion (in the small parameter $\varepsilon$) of an embedding

$$e:(\mathbb{R}/2\pi\mathbb{Z})^{m}\to\mathbb{R}^{M}\,$$

of the persisting invariant torus $\mathbb{T_{\varepsilon}}$. In addition, it computes an expansion of the dynamics on $\mathbb{T}_{\varepsilon}$ in local coordinates, in the form of a so-called “reduced phase vector field”

$${\bf f}:(\mathbb{R}/2\pi\mathbb{Z})^{m}\to\mathbb{R}^{m}$$

on the standard torus $(\mathbb{R}/2\pi\mathbb{Z})^{m}$. We find these $e$ and ${\bf f}$ by solving a so-called “conjugacy equation”. Our method is thus inspired by the work of De la Llave et al. [13], who popularised the idea of finding invariant manifolds by solving conjugacy equations. In fact, this idea was used in [6] to design a quadratically convergent iterative scheme for finding normally hyperbolic invariant tori. However, in [6] these tori are required to carry Diophantine quasi-periodic motion, not only before but also after the perturbation.

The phase reduction method presented in this paper is more similar in nature to the parametrisation method developed in [19]. There the idea of parametrisation is used to calculate expansions of slow manifolds and their flows in geometric singular perturbation problems [30]. Just like the method in [19], the phase reduction method presented here yields asymptotic expansions to finite order, but it poses no restrictions on the nature of the dynamics on the invariant torus.

We now sketch the idea behind our method. Let us write ${\bf F}_{0}$ for the vector field on $\mathbb{R}^{M}=\mathbb{R}^{M_{1}}\times\ldots\times\mathbb{R}^{M_{m}}$ that governs the dynamics of the uncoupled oscillators in (1.1), that is,

$\displaystyle{\bf F}_{0}(x_{1},\ldots,x_{m}):=(F_{1}(x_{1}),\ldots,F_{m}(x_{m}))\,.$

(1.2)

Our starting point is an embedding of the invariant torus $\mathbb{T}_{0}$ for this ${\bf F}_{0}$. Recall our assumption that every $F_{j}$ possesses a hyperbolic periodic orbit $X_{j}(t)$ of minimal period $T_{j}>0$. We denote the frequency of this orbit by $\omega_{j}:=\frac{2\pi}{T_{j}}$. An obvious embedding of $\mathbb{T}_{0}$ is the map $e_{0}:(\mathbb{R}/2\pi\mathbb{Z})^{m}\to\mathbb{R}^{M}$ defined by

$\displaystyle e_{0}(\phi)=e_{0}(\phi_{1},\ldots,\phi_{m}):=\left(X_{1}\left(\omega_{1}^{-1}\phi_{1}\right),\ldots,X_{m}\left(\omega_{m}^{-1}\phi_{m}\right)\right)\,.$

(1.3)

In fact, this $e_{0}$ sends the periodic or quasi-periodic solutions of the ODEs

$$\dot{\phi}=\omega:=(\omega_{1},\ldots,\omega_{m})$$

on $(\mathbb{R}/2\pi\mathbb{Z})^{m}$ to integral curves of ${\bf F}_{0}$. In other words—see also Lemma 2.1 below—it satisfies the conjugacy equation

$$e_{0}^{\prime}\cdot\omega={\bf F}_{0}\circ e_{0}\,.$$

The idea is now that we search for an asymptotic approximation of an embedding of the persisting torus $\mathbb{T}_{\varepsilon}$ by solving a similar conjugacy equation. We do this by making a series expansion ansatz for such an embedding, of the form

$$e=e_{0}+\varepsilon e_{1}+\varepsilon^{2}e_{2}+\ldots:(\mathbb{R}/2\pi\mathbb{Z})^{m}\to\mathbb{R}^{M}\,,$$

as well as for a reduced phase vector field

$${\bf f}=\omega+\varepsilon{\bf f}_{1}+\varepsilon^{2}{\bf f}_{2}+\ldots:(\mathbb{R}/2\pi\mathbb{Z})^{m}\to\mathbb{R}^{m}\,.$$

Indeed, writing ${\bf F}={\bf F_{0}}+\varepsilon{\bf F}_{1}:\mathbb{R}^{M}\to\mathbb{R}^{M}$, with ${\bf F}_{0}$ as above, and

$${\bf F}_{1}(x):=(G_{1}(x),\ldots,G_{m}(x))\,$$

denoting the coupled part of (1.1), we have that $e$ maps integral curves of ${\bf f}$ to solutions of (1.1), exactly when the conjugacy equation

$$e^{\prime}\cdot{\bf f}={\bf F}\circ e$$

holds. If this is the case, then $\mathbb{T}_{\varepsilon}=e((\mathbb{R}/2\pi\mathbb{Z})^{m})$ is the persisting invariant torus, whereas the vector field ${\bf f}$ on $(\mathbb{R}/2\pi\mathbb{Z})^{m}$ represents the dynamic on $\mathbb{T}_{\varepsilon}$ in local coordinates, that is, it determines the reduced phase dynamics.

We will see that the conjugacy equation for $(e,{\bf f})$ translates into a sequence of iterative equations for $(e_{1},{\bf f}_{1}),(e_{2},{\bf f}_{2}),\ldots$. We will show how to solve these iterative equations, which then allows us to compute the expansions for $e$ and ${\bf f}$ to any desired order in the small parameter. Because the embedding of the torus $\mathbb{T}_{\varepsilon}$ is not unique, neither are the solutions $(e_{j},{\bf f}_{j})$ to these iterative equations. We characterize the extent to which one is free to choose these solutions, and we show how this freedom can be exploited to obtain ${\bf f}_{j}$ that are in normal form. This means that “nonresonant” terms have been removed from the reduced phase equations to high order.

A crucial requirement for the solvability of the iterative equations is that the torus $\mathbb{T}_{0}$ is reducible. Reducibility is a property of the unperturbed dynamics normal to $\mathbb{T}_{0}$. We shall define it at the hand of an embedding of the so-called fast fibre bundle of $\mathbb{T}_{0}$. We call such an embedding a fast fibre map. The fast fibre map is an important ingredient of our method. An invariant torus for an uncoupled oscillator system is always reducible. We show in Section 5 how, in this case, the fast fibre map can be obtained from the Floquet decompositions of the fundamental matrix solutions of the periodic orbits $X_{j}(t)$. We remark that by using fast fibre maps, we are able to avoid the use of isochrons [12] to characterise the dynamics normal to $\mathbb{T}_{0}$. Our parametrisation method is therefore not restricted to the case where the periodic orbits $X_{j}(t)$ are stable limit cycles—it suffices if they are hyperbolic. We also stress that our method is not restricted to weakly coupled oscillator systems: it applies whenever the unperturbed embedded torus $\mathbb{T}_{0}$ is quasi-periodic, normally hyperbolic and reducible.

This paper is organised as follows. In section 2 we discuss the conjugacy problem for $(e,{\bf f})$ in more detail, and derive the iterative equations for $(e_{j},{\bf f}_{j})$. In section 3 we introduce fast fibre maps and use them to define when an embedded (quasi-)periodic torus is reducible. In section 4 we explain how the fast fibre map can be used to solve the iterative equations for $(e_{j},{\bf f}_{j})$. We give formulas for the solutions, and discuss their properties. Section 5 shows how to compute the fast fibre map for a coupled oscillator system, treating the Stuart-Landau oscillator as an example. We finish with an application/illustration of our method in section 6, in which we prove that remote synchronisation occurs in a chain of weakly coupled Stuart-Landau oscillators.

We start this section with a proof of our earlier claim about the embedding $e_{0}$. In the formulation of Lemma 2.1 below, we use the notation

$\displaystyle\partial_{\omega}e_{0}:=e_{0}^{\prime}\cdot\omega=\left.\frac{d}{ds}\right|_{s=0}\!\!\!\!\!\!\!e_{0}(\,\cdot+s\omega)\,$

(2.1)

for the (directional) derivative of $e_{0}$ in the direction of the vector $\omega\in\mathbb{R}^{m}$. Like $e_{0}$ itself, $\partial_{\omega}e_{0}$ is a smooth map from $(\mathbb{R}/2\pi\mathbb{Z})^{m}$ to $\mathbb{R}^{M}$.

At this point we temporarily abandon the setting of coupled oscillators and consider a general ODE $\dot{x}={\bf F}_{0}(x)$ defined by a smooth vector field ${\bf F}_{0}:\mathbb{R}^{M}\to\mathbb{R}^{M}$. That is, we do not assume that this ODE decouples into mutually independent ODEs. However, we will assume throughout this paper that ${\bf F}_{0}$ possesses a normally hyperbolic periodic or quasi-periodic invariant torus $\mathbb{T}_{0}$ which admits an embedding $e_{0}:(\mathbb{R}/2\pi\mathbb{Z})^{m}\to\mathbb{R}^{M}$ that semi-conjugates the constant vector field $\omega$ on $(\mathbb{R}/2\pi\mathbb{Z})^{m}$ to ${\bf F}_{0}$. In other words, we assume that $e_{0}$ and ${\bf F}_{0}$ satisfy

$\displaystyle\partial_{\omega}e_{0}={\bf F}_{0}\circ e_{0}\,,$

(2.2)

just as in Lemma 2.1. We return to coupled oscillator systems in section 5.

We now study any smooth perturbation of ${\bf F}_{0}$ of the form

$${\bf F}={\bf F}(x)={\bf F}_{0}(x)+\varepsilon\,{\bf F}_{1}(x)+\varepsilon^{2}\,{\bf F}_{2}(x)+\ldots:\ \mathbb{R}^{M}\to\mathbb{R}^{M}\,.$$

Fénichel’s theorem [8, 30] guarantees that, for $0\leq\varepsilon\ll 1$, the perturbed ODE $\dot{x}={\bf F}(x)$ admits an invariant torus $\mathbb{T}_{\varepsilon}$ close to $\mathbb{T}_{0}$, that depends smoothly on $\varepsilon$. Our strategy to find $\mathbb{T}_{\varepsilon}$ will be to search for an embedding $e:(\mathbb{R}/2\pi\mathbb{Z})^{m}\to\mathbb{R}^{M}$ close to $e_{0}$, and a reduced vector field ${\bf f}:(\mathbb{R}/2\pi\mathbb{Z})^{m}\to\mathbb{R}^{m}$ close to $\omega$ satisfying the conjugacy equation

$\displaystyle\mathfrak{C}(e,{\bf f}):=e^{\prime}\cdot{\bf f}-{\bf F}\circ e=0\,.$

(2.3)

Any solution $(e,{\bf f})$ to (2.3) indeed yields an embedded ${\bf F}$-invariant torus $\mathbb{T}_{\varepsilon}:=e((\mathbb{R}/2\pi\mathbb{Z})^{m})\subset\mathbb{R}^{M}$, as we see from (2.3) that at any point $x=e(\phi)\in\mathbb{T}_{\varepsilon}$ the vector ${\bf F}(x)$ lies in the image of the derivative $e^{\prime}(\phi)$, and is thus tangent to $\mathbb{T}_{\varepsilon}$. Moreover, $e$ semi-conjugates ${\bf f}$ to ${\bf F}$, that is, ${\bf f}$ is the restriction of ${\bf F}$ to $\mathbb{T}_{\varepsilon}$ represented in (or “pulled back to”) the local coordinate chart $(\mathbb{R}/2\pi\mathbb{Z})^{m}$.

As explained in the introduction, we try to find solutions to (2.3) by making a series expansion ansatz

$$e=e_{0}+\varepsilon e_{1}+\varepsilon^{2}e_{2}+\ldots\ \mbox{and}\ {\bf f}=\omega+\varepsilon{\bf f}_{1}+\varepsilon^{2}{\bf f}_{2}+\ldots$$

for $e_{1},e_{2},\ldots:(\mathbb{R}/2\pi\mathbb{Z})^{m}\to\mathbb{R}^{M}$ and ${\bf f}_{1},{\bf f}_{2},\ldots:(\mathbb{R}/2\pi\mathbb{Z})^{m}\to\mathbb{R}^{m}$. Substitution of this ansatz in (2.3), and Taylor expansion to $\varepsilon$, yields the following list of recursive equations for the $e_{j}$ and ${\bf f}_{j}$:

$\displaystyle\begin{array}[]{ccc}(\partial_{\omega}-{\bf F}_{0}^{\prime}\circ e_{0})\cdot e_{1}+e_{0}^{\prime}\cdot{\bf f}_{1}=&{\bf F}_{1}\circ e_{0}&=:{\bf G}_{1}\\ (\partial_{\omega}-{\bf F}_{0}^{\prime}\circ e_{0})\cdot e_{2}+e_{0}^{\prime}\cdot{\bf f}_{2}=&{\bf F}_{2}\circ e_{0}+({\bf F}_{1}^{\prime}\circ e_{0})\cdot e_{1}\\ &+\frac{1}{2}({\bf F}_{0}^{\prime\prime}\circ e_{0})(e_{1},e_{1})-e_{1}^{\prime}\cdot{\bf f}_{1}&=:{\bf G}_{2}\\ \vdots&\vdots&\vdots\\ (\partial_{\omega}-{\bf F}_{0}^{\prime}\circ e_{0})\cdot e_{j}+e_{0}^{\prime}\cdot{\bf f}_{j}=&\ldots&=:{\bf G}_{j}\\ \vdots&\vdots&\vdots\end{array}$

(2.10)

Here, each ${\bf G}_{j}:(\mathbb{R}/2\pi\mathbb{Z})^{m}\to\mathbb{R}^{M}$ is an “inhomogeneous term” that can iteratively be determined and depends on ${\bf F}_{1},\ldots,{\bf F}_{j},{\bf f}_{1},\ldots,{\bf f}_{j-1}$ and $e_{1},\ldots,e_{j-1}$. Concretely, ${\bf G}_{j}$ is given by

$\displaystyle\!\!\!\!{\bf G}_{j}\!:=\!\!\left.\frac{1}{j!}\frac{d^{j}}{d\varepsilon^{j}}\right|_{\varepsilon=0}\!\!\!\!\!\!\begin{array}[]{c}\!\!\!({\bf F}_{0}+\varepsilon{\bf F}_{1}+\ldots+\varepsilon^{j}{\bf F}_{j})(e_{0}+\varepsilon e_{1}+\ldots+\varepsilon^{j-1}e_{j-1})\\ -(e_{0}+\varepsilon e_{1}\ldots+\varepsilon^{j-1}e_{j-1})^{\prime}\cdot(\omega+\varepsilon{\bf f}_{1}+\ldots+\varepsilon^{j-1}{\bf f}_{j-1})\end{array}\!\!\!\!\!\!.$

(2.13)

Explicit formulas for ${\bf G}_{1}$ and ${\bf G}_{2}$ are given in (2.10). Note that equations (2.10) are all of the form

$\displaystyle\mathfrak{c}(e_{j},{\bf f}_{j})={\bf G}_{j}\ \mbox{for}\ j=1,2,\ldots\,,$

(2.14)

in which

$\displaystyle\mathfrak{c}(e_{j},{\bf f}_{j}):=(\partial_{\omega}-{\bf F}_{0}^{\prime}\circ e_{0})\cdot e_{j}+e_{0}^{\prime}\cdot{\bf f}_{j}\,$

(2.15)

is the linearisation of the operator $\mathfrak{C}$ defined in (2.3) at the point $(e,{\bf f})=(e_{0},\omega)$, where $\varepsilon=0$. This linearisation $\mathfrak{c}$ is not invertible, but we will see that $\mathfrak{c}$ is surjective under the assumption that $\mathbb{T}_{0}$ is reducible. This implies that equations (2.10) can iteratively be solved.

As was indicated in Remarks 2 and 3, the solutions to the iterative equations $\mathfrak{c}(e_{j},{\bf f}_{j})={\bf G}_{j}$ are not unique. However, we show in section 4 that solutions can be found if we assume that the unperturbed torus $\mathbb{T}_{0}$ is reducible. We define this concept by means of a parametrisation of the linearised dynamics of ${\bf F}_{0}$ normal to $\mathbb{T}_{0}$. But we start with the observation that the linearised dynamics tangent to $\mathbb{T}_{0}$ is trivial. Recall that if $e_{0}:(\mathbb{R}/2\pi\mathbb{Z})^{m}\to\mathbb{R}^{M}$ is an embedding of $\mathbb{T}_{0}\subset\mathbb{R}^{M}$, then the tangent map ${\bf T}e_{0}:(\mathbb{R}/2\pi\mathbb{Z})^{m}\times\mathbb{R}^{m}\to\mathbb{R}^{M}\times\mathbb{R}^{M}$ defined by

$\displaystyle{\bf T}e_{0}(\phi,u)=(e_{0}(\phi),e_{0}^{\prime}(\phi)\cdot u)$

(3.1)

is an embedding as well. Its image is the tangent bundle ${\bf T}\mathbb{T}_{0}\subset\mathbb{R}^{M}\times\mathbb{R}^{M}$.

We finish this section with an alternative characterisation of property $\it ii)$ in Definition 3.2.

We now return to solving the iterative equations (2.10), assuming from here on out that $\mathbb{T}_{0}$ is an embedded (quasi-)periodic reducible and normally hyperbolic invariant torus for ${\bf F}_{0}$. The main result of this section can be summarised (at this point still somewhat imprecisely) as follows.

To prove the theorem, recall that (because $\mathbb{T}_{0}$ is reducible) we have at our disposal a fast fibre map ${\bf N}e_{0}$ for $\mathbb{T}_{0}$, defined by a family of injective matrices $N=N(\phi)$ that satisfies $\mathbb{R}^{M}={\rm im}\,e_{0}^{\prime}(\phi)\oplus{\rm im}\,N(\phi)$ for every $\phi\in(\mathbb{R}/2\pi\mathbb{Z})^{m}$. This enables us to make the ansatz

$\displaystyle\underbrace{e_{j}(\phi)}_{\footnotesize\begin{array}[]{c}\in\\ \mathbb{R}^{M}\end{array}}\ \ =\!\!\underbrace{e_{0}^{\prime}(\phi)}_{\footnotesize\begin{array}[]{c}\in\\ \mathcal{L}(\mathbb{R}^{m},\mathbb{R}^{M})\end{array}}\!\!\!\!\!\cdot\ \underbrace{{\bf g}_{j}(\phi)}_{\footnotesize\begin{array}[]{c}\in\\ \mathbb{R}^{m}\end{array}}\ \ +\!\!\!\!\!\underbrace{N(\phi)}_{\footnotesize\begin{array}[]{c}\in\\ \mathcal{L}(\mathbb{R}^{M-m},\mathbb{R}^{M})\end{array}}\!\!\!\!\!\!\!\cdot\ \underbrace{{\bf h}_{j}(\phi)}_{\footnotesize\begin{array}[]{c}\in\\ \mathbb{R}^{M-m}\end{array}}\,,$

(4.11)

for (unknown) smooth functions ${\bf g}_{j}:(\mathbb{R}/2\pi\mathbb{Z})^{m}\to\mathbb{R}^{m}$ and ${\bf h}_{j}:(\mathbb{R}/2\pi\mathbb{Z})^{m}\to\mathbb{R}^{M-m}$. This ansatz decomposes $e_{j}$ into components in the direction of the tangent bundle ${\bf T}e_{0}$ and the fast fibre bundle ${\bf N}e_{0}$.

Applying $\pi$ and $1-\pi$ to (4.12) produces, respectively,

	$\displaystyle e_{0}^{\prime}\cdot(\partial_{\omega}{\bf g}_{j}+{\bf f}_{j})=\pi\cdot{\bf G}_{j}\,,$
	$\displaystyle N\cdot(\partial_{\omega}-L)({\bf h}_{j})=(1-\pi)\cdot{\bf G}_{j}\,.$

Because $e_{0}^{\prime}(\phi)$ and $N(\phi)$ are injective, these equations are equivalent to

$\displaystyle\begin{array}[]{rll}\partial_{\omega}{\bf g}_{j}+{\bf f}_{j}=&(e_{0}^{\prime})^{+}\cdot\pi\cdot{\bf G}_{j}&=:U_{j}\,,\\ (\partial_{\omega}-L)({\bf h}_{j})=&N^{+}\cdot(1-\pi)\cdot{\bf G}_{j}&=:V_{j}\,.\end{array}$

(4.15)

Here, $A^{+}:=(A^{T}A)^{-1}A^{T}$ denotes the Moore-Penrose pseudo-inverse, which is well-defined for an injective linear map $A$. Clearly, $(e_{0}^{\prime})^{+}$ and $N^{+}$ depend smoothly on $\phi\in(\mathbb{R}/2\pi\mathbb{Z})^{m}$. We give these equations a special name.