Optimizing stability of mutual synchronization between a pair of limit-cycle oscillators with weak cross coupling

TL;DR

This paper develops a method to optimize the linear stability of synchronized states in weakly coupled limit-cycle oscillators by designing the coupling matrix, enhancing synchronization robustness.

Contribution

It introduces a phase reduction-based approach to derive the optimal coupling matrix for maximizing stability under specific constraints.

Findings

Optimal coupling matrices improve synchronization stability.

The method applies to various types of limit-cycle oscillators.

Enhanced stability demonstrated through numerical examples.

Abstract

We consider optimization of linear stability of synchronized states between a pair of weakly coupled limit-cycle oscillators with cross coupling, where different components of state variables of the oscillators are allowed to interact. On the basis of the phase reduction theory, the coupling matrix between different components of the oscillator states that maximizes the linear stability of the synchronized state under given constraints on overall coupling intensity and on stationary phase difference is derived. The improvement in the linear stability is illustrated by using several types of limit-cycle oscillators as examples.

Figures

Equations

$${\dot{\mathbf{X}}}_1(t)$$

$${\dot{\mathbf{X}}}_2(t)$$

$${\mathbf{F}}_{1,2}(\mathbf{X})=\mathbf{F}(\mathbf{X})+ϵ{\mathbf{f}}_{1,2}(\mathbf{X}),$$

$${\dot{\theta }}_1(t)$$

$${\dot{\theta }}_2(t)$$

$${\omega }_{1,2}$$

$$=\omega +ϵ⟨\mathbf{Z}(\psi )\cdot {\mathbf{f}}_{1,2}({\mathbf{X}}_0(\psi )){⟩}_{\psi }$$

$$\Gamma (\varphi )$$

$$=⟨\mathbf{Z}(\varphi +\psi )\cdot K\{{\mathbf{X}}_0(\psi )-{\mathbf{X}}_0(\varphi +\psi )\}{⟩}_{\psi }.$$

$$⟨A(\psi ){⟩}_{\psi }=\frac{1}{2\pi }{\int }_0^{2\pi }A(\psi )d\psi ,$$

$$\Gamma (\varphi )$$

$$=\mathrm{Tr}(KW(\varphi {)}^{\mathrm{T}}),$$

$$W(\varphi )=⟨\mathbf{Z}(\varphi +\psi )\otimes \{{\mathbf{X}}_0(\psi )-{\mathbf{X}}_0(\varphi +\psi )\}{⟩}_{\psi }$$

$$\dot{\varphi }=ϵ\{\Delta \omega +{\Gamma }_a(\varphi )\},{\Gamma }_a(\varphi )=\Gamma (\varphi )-\Gamma (-\varphi ).$$

$$\Delta \omega +{\Gamma }_a({\varphi }^{\ast })=0$$

$${{\Gamma }_a^{\prime }({\varphi }^{\ast })=\frac{d}{d\varphi }{\Gamma }_a(\varphi )|}_{\varphi ={\varphi }^{\ast }}<0,$$

$$V(\varphi )=W(\varphi )-W(-\varphi ),$$

$${\Gamma }_a(\varphi )=\mathrm{Tr}\{KV(\varphi {)}^{\mathrm{T}}\},\frac{d}{d\varphi }{\Gamma }_a(\varphi )=\mathrm{Tr}\{K{V}^{\prime }(\varphi {)}^{\mathrm{T}}\},$$

$$V(\varphi )=⟨\{\mathbf{Z}(\varphi +\psi )-\mathbf{Z}(-\varphi +\psi )\}\otimes {\mathbf{X}}_0(\psi ){⟩}_{\psi },$$

$$V^{\prime }(\varphi )$$

$$=⟨\{{\mathbf{Z}}^{\prime }(\varphi +\psi )+{\mathbf{Z}}^{\prime }(-\varphi +\psi )\}\otimes {\mathbf{X}}_0(\psi ){⟩}_{\psi },$$

$$W^{\prime }(\varphi )$$

$$-ϵ{\Gamma }_a^{\prime }(0),$$

$$\Vert K{\Vert }^2=P(P>0).$$

$$-{\Gamma }_a^{\prime }(0)=-{\frac{d}{d\varphi }{\Gamma }_a(\varphi )|}_{\varphi =0},$$

$$S(K,\lambda )$$

$$\displaystyle=-{\rm Tr\ }\Big{(}KV^{\prime}(0)^{\rm T}\Big{)}+\lambda\Big{(}\|K\|^{2}-P\Big{)},$$

$$\frac{\partial }{\partial K}S(K,\lambda )=$$

$$K=\frac{1}{2\lambda }{V}^{\prime }(0).$$

$$\lambda =\pm \frac{1}{2\sqrt{P}}\Vert V^{\prime }(0)\Vert .$$

$$K_{\mathrm{opt}}=-\sqrt{P}\frac{V^{\prime }(0)}{\Vert V^{\prime }(0)\Vert }.$$

$${\Gamma }_a(\varphi )=\mathrm{Tr}\{K_{\mathrm{opt}}V(\varphi {)}^{\mathrm{T}}\}=-\frac{\sqrt{P}}{\Vert V^{\prime }(0)\Vert }\mathrm{Tr}\{V^{\prime }(0)V(\varphi {)}^{\mathrm{T}}\}$$

$$-{\Gamma }_a^{\prime }(0)$$

$$-{\Gamma }_a^{\prime }({\varphi }^{\ast })=-{\frac{d}{d\varphi }{\Gamma }_a(\varphi )|}_{\varphi ={\varphi }^{\ast }},$$

Full text

I Introduction

Synchronization of nonlinear oscillators is widely observed and often plays important functional roles in a variety of real-world systems ref:pikovsky01 ; ref:strogatz03 ; strogatz15 ; ref:winfree80 ; Kuramoto ; ref:hoppensteadt97 ; ref:ermentrout10 ; Kiss ; Kiss2 . Exploration of efficient methods for realizing stable synchronization between coupled oscillators or between oscillators and driving signals is both fundamentally and practically important. Improvement in the efficiency of collective synchronization in networks of coupled oscillators has been extensively studied in the literature ref:tanaka08 ; ref:yanagita10 ; ref:yanagita12 ; ref:yanagita14 ; ref:skardal14 ; ref:skardal16 ; ref:nishikawa06a ; ref:nishikawa06b ; ref:nishikawa10 ; Ravoori for both Kuramoto-type phase models and chaotic oscillators, where optimization of coupling networks connecting the oscillators has been the main target.

In the analysis of synchronization dynamics between weakly coupled nonlinear oscillators undergoing limit-cycle oscillations, the phase reduction theory has played a dominant role ref:winfree80 ; Kuramoto ; ref:hoppensteadt97 ; ref:ermentrout10 ; ref:brown04 ; ref:nakao16 ; ref:ashwin16 . It allows us to simplify the dynamics of a pair of limit-cycle oscillators with weak coupling to a simple scalar equation for their phase difference. The phase reduction theory, originally developed for finite-dimensional smooth limit-cycle oscillators, has recently been generalized to non-conventional limit-cycling systems such as collectively oscillating populations of coupled oscillators Kawamura1 , systems with time delay Kotani ; Novicenko1 ; Novicenko2 , reaction-diffusion systems Nakao , oscillatory fluid convection Kawamura2 , and hybrid dynamical systems Shirasaka . Recently, methods for optimizing periodic external driving signals for efficient injection locking and controlling of a single nonlinear oscillator (or a population of uncoupled oscillators) have also been proposed on the basis of the phase reduction theory ref:moehlis06 ; ref:harada10 ; ref:dasanayake11 ; ref:zlotnik12 ; ref:zlotnik13 ; ref:pikovsky15 ; ref:tanaka14a ; ref:tanaka14b ; ref:tanaka15 ; ref:zlotnik2016 ; ref:hasegawa14a ; ref:hasegawa14b . In this study, we consider a pair of coupled limit-cycle oscillators and try to optimize the linear stability of the synchronized state using the phase reduction theory.

In the analysis of mutual synchronization of coupled oscillators, linear diffusive coupling between the oscillators is a common setup. However, in most cases, only the same vector component of the state variables can interact between the oscillators and different vector components of the oscillator states are usually not allowed to interact. In this study, we analyze a pair of oscillators with weak cross coupling, where different vector components of the oscillator states are allowed to interact, that is, differences in each vector component of the oscillator states can be feed-backed to every other component with a linear gain specified by a coupling matrix, and optimize the coupling matrix so that the linear stability of the mutually synchronized state is maximized.

We use the phase reduction theory to simplify the dynamics of a pair of weakly coupled limit-cycle oscillators to a scalar equation for the phase difference, and use the method of Lagrange multipliers to derive the optimal coupling matrix for the cases with and without frequency mismatch between the oscillators. Using three examples of simple limit-cycle oscillators, we illustrate that the linear stability of the synchronized state is actually improved and also that the stationary phase difference can be controlled by appropriately choosing the coupling matrix.

This paper is organized as follows: in Sec. II, we introduce the coupled-oscillator model and derive the equation for the phase difference by using the phase reduction theory. In Sec. III, we formulate the optimization problem for improving linear stability of the phase-locked states. In Sec. IV, the theoretical results are illustrated by several examples of limit-cycle oscillators. Sec. V gives summary and discussion.

II Model

In this section, we introduce a pair of nearly identical limit-cycle oscillators with weak cross coupling, reduce the dynamical equations to coupled phase equations by using the phase reduction theory ref:winfree80 ; Kuramoto ; ref:hoppensteadt97 ; ref:brown04 ; ref:ermentrout10 ; ref:nakao16 ; ref:ashwin16 , and derive the equation for the phase difference.

II.1 A pair of cross-coupled oscillators

We consider a pair of weakly and symmetrically coupled, nearly identical limit-cycle oscillators described by

[TABLE]

where ${\bm{X}}_{1}$ and ${\bm{X}}_{2}$ are the $m$-dimensional state vectors of the oscillators $1$ and $2$, respectively, ${\bm{F}}_{1}$ and ${\bm{F}}_{2}$ are $m$-dimensional vector-valued functions representing the dynamics of the oscillators, $K$ is a $m\times m$ matrix of coupling intensities between the components of the state variables, and $\epsilon$ is a small positive parameter $(0<\epsilon\ll 1)$ indicating that the interaction is sufficiently small.

Here, although the oscillators are “diffusively” coupled, we assume that the matrix $K$ is generally not diagonal and can possess non-diagonal elements. That is, differences in each vector component of the oscillator states are returned to other components as feedback signals with appropriate gains. Therefore, different components of the state variables of the oscillators can mutually interact. This gives the possibility to improve the stability of the synchronized state by adjusting the non-diagonal elements of the coupling matrix, exceeding the stability that is achievable only with the diagonal coupling. We assume linear diffusive coupling in the following, but the argument can be straightforwardly generalized to nonlinear coupling; see Sec. V.

We assume that the properties of the oscillators are nearly identical and their difference is $O(\epsilon)$. That is, the functions ${\bm{F}}_{1,2}$ can be split into a common part ${\bm{F}}$ and deviations ${\bm{f}}_{1,2}$ as

[TABLE]

where ${\bm{F}}$, ${\bm{f}}_{1}$, and ${\bm{f}}_{2}$ are assumed to be $O(1)$. We also assume that the common part of the oscillator dynamics, $\dot{\bm{X}}(t)={\bm{F}}({\bm{X}})$, possesses a stable limit-cycle solution ${\bm{X}}_{0}(t)={\bm{X}}_{0}(t+T)$ of period $T$ and frequency $\omega=2\pi/T$, and that the dynamics of the oscillator is only slightly deformed and persists even if small perturbations from the deviations ${\bm{f}}_{1,2}$ and mutual coupling are introduced. These assumptions are necessary for the phase reduction that we rely on in the present study.

II.2 Phase reduction

Under the above assumptions, we can simplify the dynamics of the coupled oscillators to coupled phase equations by applying the phase reduction theory ref:winfree80 ; Kuramoto ; ref:hoppensteadt97 ; ref:brown04 ; ref:ermentrout10 ; ref:nakao16 ; ref:ashwin16 . That is, we introduce a phase $\theta$ ($0\leq\theta<2\pi$) of the oscillator state near the limit-cycle solution ${\bm{X}}_{0}(t)$ that increases with a constant frequency $\omega$ in the absence of perturbations, and represent the oscillator state on the limit cycle as a function of the phase $\theta(t)$ as ${\bm{X}}_{0}(\theta(t))$.

In the present case, we introduce phase variables $\theta_{1,2}$ of the two oscillators, represent the oscillator states near the limit-cycle orbit as ${\bm{X}}_{1,2}(t)={\bm{X}}_{0}(\theta_{1,2}(t))+O(\epsilon)$ as functions of $\theta_{1,2}(t)$ at $t$, and approximately describe their dynamics by using only $\theta_{1,2}$. By following the standard phase reduction and averaging procedures, we can derive a pair of coupled phase equations, which is correct up to $O(\epsilon)$, as

[TABLE]

The frequencies $\omega_{1,2}$ of the oscillators are given by

[TABLE]

and the phase coupling function $\Gamma(\phi)$ is given by

[TABLE]

Here, we introduced an abbreviation for the average over phase from [math] to $2\pi$,

[TABLE]

where $A(\psi)$ is a $2\pi$-periodic function of $\psi$. In the following, without loss of generality, we assume that $\omega_{1}\geq\omega_{2}$, and denote the frequency difference between the oscillators as $\epsilon\Delta\omega=\omega_{1}-\omega_{2}\geq 0$, where $\Delta\omega$ is $O(1)$.

The function ${\bm{Z}}(\theta)$ in Eqs. (7) and (9) is a phase sensitivity function of the limit cycle ${\bm{X}}_{0}(\theta)$ of the common part, $\dot{\bm{X}}(t)={\bm{F}}({\bm{X}})$. It is given by a $2\pi$-periodic solution to the adjoint equation $\partial{\bm{Z}}(\theta)/\partial\theta=-J(\theta)^{\rm T}{\bm{Z}}(\theta)$, where $J(\theta)$ is a Jacobi matrix of the vector field ${\bm{F}}({\bm{X}})$ at ${\bm{X}}={\bm{X}}_{0}(\theta)$ and T denotes the matrix transpose, and is normalized as ${\bm{Z}}(\theta)\cdot{\bm{F}}({\bm{X}}_{0}(\theta))=\omega$. By using the adjoint method by Ermentrout ref:ermentrout10 ; ref:brown04 ; ref:nakao16 , i.e., by backwardly evolving the adjoint equation with occasional renormalization, ${\bm{Z}}(\theta)$ can be calculated numerically.

For convenience, we rewrite the phase coupling function as

[TABLE]

where

[TABLE]

is a correlation matrix between the vector components of the phase sensitivity function and the state difference between the oscillators. Here, the symbol $\otimes$ represents a tensor product and ${\rm Tr}$ denotes the trace of a matrix. See Appendix A for the definition and related matrix formulas. Because ${\bm{X}}_{0}(\theta)$ and ${\bm{Z}}(\theta)$ are $2\pi$-periodic functions, $\Gamma(\phi)$ and $W(\phi)$ are also $2\pi$-periodic.

II.3 Stability of the synchronized state

From Eq. (5), the phase difference $\phi=\theta_{1}-\theta_{2}$ (restricted to $-\pi\leq\phi\leq\pi$ hereafter) between the two oscillators obeys

[TABLE]

Here, $\Gamma_{a}(\phi)$ is the antisymmetric part of the phase coupling function $\Gamma(\phi)$; it is also $2\pi$-periodic and satisfies $\Gamma_{a}(0)=\Gamma_{a}(\pm\pi)=0$. Therefore, if $\Delta\omega$ satisfies $-{\rm max}_{\phi}\Gamma_{a}(\phi)<\Delta\omega<-{\rm min}_{\phi}\Gamma_{a}(\phi)$, Eq. (14) has at least one stable fixed point at the phase differences satisfying $\Delta\omega+\Gamma_{a}(\phi)=0$. We denote one of such fixed points as $\phi^{*}$.

From Eq. (14), the linear stability of $\phi^{*}$ is given by $\epsilon\Gamma_{a}^{\prime}(\phi^{*})$, where $\Gamma_{a}^{\prime}(\phi^{*})$ is the slope of $\Gamma_{a}(\phi)$ at $\phi=\phi^{*}$. Thus, if a fixed point $\phi=\phi^{*}$ satisfies

[TABLE]

and

[TABLE]

the phase difference $\phi$ can take a stationary phase difference, and the two oscillators can mutually synchronize, or phase-locked to each other, with a stable phase difference $\phi=\phi^{*}$ (within the phase-reduction approximation).

By defining a new matrix

[TABLE]

the antisymmetric part of the phase coupling function and its slope can be expressed as

[TABLE]

where $V^{\prime}(\phi)$ represents derivative of $V(\phi)$ with respect to $\phi$. Using $2\pi$-periodicity of ${\bm{Z}}(\theta)$ and ${\bm{X}}_{0}(\theta)$, the matrices $V(\phi)$ and $V^{\prime}(\phi)$ can be expressed as

[TABLE]

and

[TABLE]

where the derivative of $W(\phi)$ with respect to $\phi$ is given by

[TABLE]

See Appendix B for the calculations. We use these expressions in the next section.

III Optimizing the coupling matrix

In this section, we derive the optimal coupling matrix $K_{\rm opt}$ for stable synchronization (phase locking) of the two oscillators.

III.1 Optimality condition and constraint on the coupling matrix

Our aim is to maximize the linear stability of the synchronized (phase-locked) state, characterized by $-\epsilon\Gamma_{a}^{\prime}(\phi^{*})$, by adjusting the coupling matrix $K$. Other types of optimality conditions for synchronization have also been considered in the literature for nonlinear oscillators driven by periodic signals, such as maximization of frequency difference between the oscillator and signal for fixed coupling intensity ref:harada10 ; ref:tanaka14a and minimization of the phase diffusion constant under the effect of noise ref:pikovsky15 , in addition to the maximization of linear stability ref:dasanayake11 ; ref:zlotnik12 ; ref:zlotnik13 that we generalize to coupled oscillators in the present study footnote .

We first consider the simple case where the two oscillators are identical and share the same frequency, and optimize the stability of the in-phase synchronized state with zero phase difference, $\phi^{*}=0$. We then consider the general case with a frequency mismatch $\Delta\omega\geq 0$ and optimize the stability of the synchronized state with a given stationary phase difference $\phi^{*}$, which is not necessarily [math]. In both cases, as a constraint on the overall connection intensity between the oscillators, we fix the Frobenius norm (see Appendix A) of the coupling matrix $K$ as $\|K\|^{2}=P$, where $P>0$ is a given constant. In the latter case, the stationary phase difference $\phi^{*}$ is also constrained.

III.2 Optimization for identical oscillators without frequency mismatch

We first consider the simple case where the oscillators are identical, ${\bm{F}}_{1}={\bm{F}}_{2}$, and their frequencies are equal to each other, $\omega_{1}=\omega_{2}=\omega$ and $\Delta\omega=0$. In this case, the in-phase and anti-phase synchronized states $\phi^{*}=0$ and $\phi^{*}=\pi$ are always stationary solutions to Eq. (14) because $\Gamma_{a}(0)=\Gamma_{a}(\pm\pi)=0$.

We thus try to find the coupling matrix $K$ that gives the maximum of linear stability of $\phi^{*}=0$,

[TABLE]

subject to the constraint on the Frobenius norm of $K$,

[TABLE]

Because $\epsilon>0$, we divide this quantity by $\epsilon$ and simply try to maximize

[TABLE]

which we also call “linear stability” for simplicity in the following. We introduce an action,

[TABLE]

where $\lambda$ is a Lagrange multiplier. The first term of $S$ represents the stability of the fixed point and the second term represents the constraint.

By differentiating $S$ by $K$ and $\lambda$, we obtain

[TABLE]

and the constraint $\|K\|^{2}=P$. Therefore, the optimal $K$ should satisfy

[TABLE]

Plugging this $K$ into Eq. (24) yields

[TABLE]

It turns out that the negative sign should be chosen (see Eq. (33)), so that the optimal coupling matrix is given by

[TABLE]

The antisymmetric part of the phase coupling function with this $K_{\rm opt}$ is given, from Eq. (18), by

[TABLE]

and the optimal linear stability of the in-phase fixed point $\phi=0$ is given by

[TABLE]

III.3 Optimization for nonidentical oscillators with frequency mismatch

We next consider the general case with nonidentical oscillators with a frequency mismatch $\Delta\omega\geq 0$. We constrain the Frobenius norm of $K$ as $\|K\|^{2}=P$ as before, and also require that the given $\phi^{*}$ satisfies Eq. (36), i.e., $\Delta\omega+\Gamma_{a}(\phi^{*})=0$, so that $\phi^{*}$ is actually the stationary phase difference of the oscillators.

We thus seek for the optimal coupling matrix $K_{\rm opt}$ that maximizes

[TABLE]

now for a given stationary phase difference $\phi^{*}$, subject to

[TABLE]

and

[TABLE]

Here, we exclude the cases with $\phi^{*}=0$ and $\phi^{*}=\pm\pi$, because these states can never be realized when $\Delta\omega>0$ as we argue later (the case with $\Delta\omega=0$ and $\phi^{*}=0$ was already considered in the previous subsection, and $\Delta\omega=0$ and $\phi^{*}=\pm\pi$ can be analyzed similarly).

Using Lagrange multipliers $\lambda$ and $\mu$, we introduce an action (in the rest of this subsection, shorthand notations $V_{*}=V(\phi_{*})$ and $V^{\prime}_{*}=V^{\prime}(\phi_{*})$ are used),

[TABLE]

Differentiating $S$ by $K$, $\lambda$, and $\mu$, we obtain

[TABLE]

and the two constraints, Eqs. (36) and (24). Thus, the optimal $K$ should satisfy

[TABLE]

and plugging this into Eq. (36) yields

[TABLE]

Solving this equation for $\mu$, we obtain

[TABLE]

and therefore

[TABLE]

where $\lambda$ has yet to be determined from the constraint on the Frobenius norm.

Plugging this $K_{\rm opt}$ into $\|K\|^{2}=P$ and using ${\rm Tr\ }(V^{\prime}V^{\rm T})={\rm Tr\ }(VV^{\prime\rm T})$, we obtain

[TABLE]

which gives

[TABLE]

It turns out that the minus sign should be chosen to maximize the linear stability (see below).

Note here that $\|V^{\prime}_{*}\|\|V_{*}\|\geq{\rm Tr\ }(V^{\prime}_{*}V_{*}^{\rm T})$ holds by the Schwartz inequality (see Appendix A), so the condition

[TABLE]

is necessary for $\lambda$ and hence $K_{\rm opt}$ to exist. Note also that

[TABLE]

should hold strictly for the existence of $K_{\rm opt}$ in Eq. (45), that is, $\|V^{\prime}_{*}\|\|V_{*}\|$ should not be equal to ${\rm Tr\ }(V^{\prime}_{*}V_{*}^{\rm T})$, because then $\lambda=0$ and $K_{\rm opt}$ does not exist. Therefore, the optimization problem cannot be solved in the case that $V^{\prime}_{*}$ and $V_{*}$ are parallel to each other.

The antisymmetric part of the phase coupling function for $K_{\rm opt}$ is given by

[TABLE]

and the maximal possible linear stability is given by

[TABLE]

Because $\|V_{*}\|^{2}>0$ and $\|V^{\prime}_{*}\|^{2}\|V_{*}\|^{2}-[{\rm Tr\ }(V^{\prime}_{*}V_{*}^{\rm T})]^{2}\geq 0$, the first term is positive only when $\lambda<0$. Therefore, the minus sign should be chosen for $\lambda$ in Eq. (47) to realize the maximal stability, and the optimal coupling matrix is given by Eq. (45) with the negative $\lambda$.

Note that even if we choose the minus sign for $\lambda$, the above quantity can still be negative if the second term on the right-hand side is negative, i.e., ${\rm Tr\ }(V^{\prime}_{*}V_{*}^{\rm T})<0$. If so, the fixed point with phase difference $\phi^{*}$ is unstable and cannot be realized. Thus, in this case, as can be shown by comparing the two terms on the right-hand side of Eq. (55), $P$ should additionally satisfy

[TABLE]

for $\phi^{*}$ to be linearly stable.

Depending on the conditions, the present optimization problem may or may not possess an appropriate solution. For example, when $\Delta\omega\neq 0$, it is impossible to realize completely in-phase ($\phi^{*}=0$) or anti-phase ($\phi^{*}=\pm\pi$) synchronization, because $\Gamma_{a}(\phi)$ satisfies $\Gamma_{a}(0)=\Gamma_{a}(\pm\pi)=0$, so $\Delta\omega+\Gamma_{a}(\phi)=0$ can never be satisfied at $\phi=0$ or $\phi=\pm\pi$. Also, when $\Delta\omega\neq 0$, it is generally difficult (very large $P$ is required) to realize the synchronized state with a stationary phase difference $\phi^{*}$ very close to [math] or $\pi$. This will be illustrated in the next section. The equation $\Delta\omega+\Gamma(\phi)=0$ may also have multiple solutions, so not only the fixed point with the given phase difference but also spurious fixed points with other phase differences may arise.

IV Examples

In this section, we illustrate the improvement in the linear stability of coupled oscillators by optimizing $K$ with a few types of limit-cycle oscillators as examples.

IV.1 Stuart-Landau oscillator

As the first example, we consider the Stuart-Landau (SL) oscillator, a normal form of the supercritical Hopf bifurcation Kuramoto . All necessary quantities can be analytically calculated for this model. The SL oscillator has a two-dimensional state variable, ${\bm{X}}=(x,y)^{\rm T}$, whose dynamics is specified by a vector field

[TABLE]

where $\alpha$ and $\beta$ are parameters. It possesses a single stable limit-cycle orbit of frequency $\omega=\alpha-\beta$ given by

[TABLE]

with $\theta(t)=\omega t$ (mod $2\pi$). The phase sensitivity function of this limit cycle can be explicitly calculated as Kuramoto ; ref:nakao16

[TABLE]

We consider a pair of symmetrically coupled SL oscillators with identical properties obeying Eq. (2), which is explicitly described by

[TABLE]

where ${\bm{X}}_{1}=(x_{1},y_{1})^{\rm T}$ and ${\bm{X}}_{2}=(x_{2},y_{2})^{\rm T}$ are the state variables of the oscillators. The frequency of the oscillators in the absence of mutual coupling is given by $\omega=\alpha-\beta$.

From Eqs. (66) and (71), the matrix $V(\phi)$ and its derivative $V^{\prime}(\phi)$ can be calculated as

[TABLE]

and

[TABLE]

so they are parallel to each other.

Because the oscillators are identical, $\omega_{1}=\omega_{2}$ and $\Delta\omega=0$, the in-phase synchronized state $\phi^{*}=0$ always exists. From Eq. (31), the stability of this state can be maximized by choosing $K$ as

[TABLE]

and the maximum possible linear stability is given, from Eq. (33), by

[TABLE]

For comparison, suppose that the coupling matrix is a multiple of the identity matrix with the same Frobenius norm as $K_{\rm opt}$, i.e.,

[TABLE]

which we call the “identity coupling” hereafter. The linear stability of $\phi^{*}=0$ with this $K_{I}$ is given by

[TABLE]

so the linear stability improves by a factor of $\sqrt{\beta^{2}+1}$ by using the optimal coupling matrix $K_{\rm opt}$ from the case with $K_{I}$.

Figure 1 shows the antisymmetric part $\Gamma_{a}(\phi)$ of the phase coupling function calculated for coupling matrices $K=K_{\rm opt}$ and $K=K_{I}$, with parameters $\alpha=3$ and $\beta=2$ and overall coupling intensity $P=0.1$. It can be seen that the linear stability $-\Gamma_{a}^{\prime}(0)$ of the in-phase fixed point $\phi^{*}=0$ is higher in the case with $K_{\rm opt}$ than in the case with $K_{I}$. Figure 2 compares the time courses of the phase difference $\phi(t)$ from the initial condition $\phi(0)=0.5$ for $K_{\rm opt}$ and $K_{I}$ obtained by direct numerical simulations of the coupled SL oscillators and by numerical integration of the reduced phase equation. Figure 3 shows the synchronization dynamics obtained numerically, where the $x$ components of the two oscillators are shown for $K_{\rm opt}$ and $K_{I}$. It can be seen that the oscillators synchronize faster in the case with $K_{\rm opt}$, reflecting higher linear stability of $\phi^{*}=0$, than in the case with $K_{I}$.

It is interesting to note that when $\beta=0$, that is, when the instantaneous frequency of the SL oscillator does not depend on its amplitude, $K_{I}$ is already optimal and no improvement can be made by introducing cross coupling between different variables of the oscillators.

For nonidentical SL oscillators with a small parameter mismatch $\Delta\omega>0$, the antisymmetric part of the phase coupling function and its derivative take the form $\Gamma_{a}(\phi)=-C\sin\phi$ and $\Gamma_{a}^{\prime}(\phi)=-C\cos\phi$ from Eqs. (18), (92) and (95), where $C$ is a constant determined by $K$ and $\beta$ ($K$ should also satisfy $\|K\|^{2}=P$). Once the stationary phase difference $\phi^{*}\neq 0$ is specified, the constant $C$ is determined as $C=\Delta\omega/\sin\phi^{*}$ and the linear stability of $\phi^{*}$ is given by a fixed value $-\Gamma_{a}^{\prime}(\phi^{*})=\Delta\omega/\tan\phi^{*}$. Therefore, we cannot consider further optimization of the stability for the coupled SL oscillators. Indeed, we cannot consider the second case in Sec. III, because $V(\phi)$ and $V^{\prime}(\phi)$, given by Eqs. (92) and (95), are strictly parallel to each other, so the Lagrange multiplier $\lambda$ vanishes and $K_{\rm opt}$ does not exist. This is a peculiar property of the SL oscillator with purely sinusoidal limit cycles and phase sensitivity functions.

IV.2 Brusselator

As the second example, we use the Brusselator model of chemical oscillations Kuramoto . It has a two-dimensional state variable, ${\bm{X}}=(x,y)^{\rm T}$, which obeys

[TABLE]

where $a$ and $b$ are parameters. When $a=1$ and $b=3$, the period of the oscillation is $\omega\approx 0.878$. Figure 4 shows the limit-cycle solution ${\bm{X}}_{0}(\theta)=(x_{0}(\theta),y_{0}(\theta))^{\rm T}$ and the phase sensitivity function ${\bm{Z}}(\theta)=(Z_{x}(\theta),Z_{y}(\theta))^{\rm T}$ for $0\leq\theta<2\pi$ obtained numerically. Other quantities such as $V(\phi)$ and $V^{\prime}(\phi)$ can also be numerically calculated from these ${\bm{X}}_{0}(\theta)$ and ${\bm{Z}}(\theta)$.

We consider a pair of Brusselators with parameters $b=3\pm\delta$, where $\delta$ is a small number representing parameter mismatch, and couple them in the same way as in the previous SL case, Eq. (89). We seek for the optimal $K_{\rm opt}$ that gives the maximum stability of the synchronized state, and compare the results for $K_{\rm opt}$ with those for the identity coupling, i.e., $K_{I}$ given by Eq. (102). The overall intensity $P$ is fixed at $P=0.1$ in the following.

We first consider the case without parameter mismatch, $\delta=0$. The frequencies of the oscillator are identical in this case, $\Delta\omega=0$. The optimal and identity coupling matrices with $\|K\|^{2}=P=0.1$ are calculated as

[TABLE]

We can see that, in the optimal case, the feedback from the difference in $y$ component to the dynamics of $x$ and $y$ components is stronger than that in the opposite direction. This reflects the waveforms of the oscillation and phase sensitivity function, in particular, that the variation in $y$ is generally larger than that in $x$, as shown in Fig. 4.

Figure 5 shows the antisymmetric part $\Gamma_{a}(\phi)$ of the phase coupling function for $K=K_{\rm opt}$ and $K=K_{I}$. The linear stability of the in-phase synchronize state $\phi^{*}=0$ is approximately $-\Gamma_{a}^{\prime}(0)=0.621$ for $K_{\rm opt}$ and $-\Gamma_{a}^{\prime}(0)=0.448$ for $K_{I}$. Figure 6 shows the evolution of phase differences for $K_{\rm opt}$ and $K_{I}$ obtained by direct numerical simulations of the coupled Brusselators and by numerical integration of the reduced phase model. The parameter $\epsilon$ is fixed at $\epsilon=0.05$ in the numerical simulations. Figure 7 shows the synchronization processes of the Brusselators, where time courses of the differences in $x$ components between the oscillators, i.e., $x_{1}-x_{2}$, are plotted for $K_{\rm opt}$ and $K_{I}$. For comparison, an exponentially decaying curve with the decay rate $\epsilon\Gamma_{a}^{\prime}(0)\ (<0)$ is also shown in each figure. It can been seen that in-phase synchronization is established faster when $K_{\rm opt}$ is used, and the exponential decay rate of the state difference matches with the linear stability $\Gamma^{\prime}_{a}(0)$ of $\phi^{*}=0$.

We next consider the case with parameter mismatch, $\delta=0.01$. The frequencies of the oscillators are $\omega_{1}\approx 0.8797$ ($b=2.99$) and $\omega_{2}\approx 0.8762$ ($b=3.01$). We assume $\epsilon=0.02$ in the following calculations, so the frequency mismatch parameter is $\Delta\omega\approx 0.175$. Using the results obtained in the previous section, we calculate the optimal coupling matrix $K_{\rm opt}$ for a given phase difference $\phi^{*}$ in $(-\pi,\pi)$.

Figure 8 shows the necessary conditions for $P$ given by Eq. (48) and Eq. (56) as functions of the phase difference $\phi^{*}$ for $\Delta\omega=0.175$, where the latter applies only when ${\rm Tr\ }(V^{\prime}_{*}V_{*}^{\rm T})<0$. Both conditions are satisfied in the non-shaded regions. We see that $P$ should not be too small and that the regions near $\phi^{*}=0$ and $\phi^{*}=\pm\pi$ are difficult to realize, as argued in the previous section.

Figure 9 shows the elements of the optimal coupling matrix $K_{\rm opt}$ and the corresponding stability of the fixed point as functions of the phase difference $\phi^{*}$. For comparison, the results for $K_{I}$, which gives a stable phase difference $\phi^{*}\approx 0.378$ and negative slope $0.487$, are also indicated in the figure. In this particular example, $K_{I}$ yields reasonably high stability close to the negative slope $0.493$ with the optimal coupling matrix $K_{\rm opt}$ at $\phi^{*}=0.378$ footnote2 . In the blank regions where the data are not shown, any of the necessary conditions is not satisfied. The stability varies with $\phi^{*}$ and, in this case, nearly anti-phase synchronized state yields the highest stability. Elements of the coupling can be positive or negative depending on $\phi^{*}$.

Figure 10 shows the antisymmetric parts $\Gamma_{a}(\phi)$ of the obtained phase coupling functions for given phase differences, $\phi^{*}=-2.5,-2.0,-1.5,-1.0,0.5,1.0,1.5$, and $2.0$. The given phase differences are actually realized with the optimal $\Gamma_{a}(\phi)$ as stable fixed points. From this figure, we can clearly see why stationary phase differences close to [math] or $\pi$ are difficult to realize, in consistent with the conditions shown in Fig. 8. Figure 11 plots the results of direct numerical simulations for coupled Brusselators using the optimal coupling matrices $K_{\rm opt}$, where the convergence of the phase differences to given values is shown.

IV.3 Lorenz model

Finally, as a simple three-dimensional example, we consider the Lorenz model in the limit-cycling regime strogatz15 , whose state variable ${\bm{X}}=(x,y,z)^{\rm T}$ evolves with the vector field

[TABLE]

with $\sigma=10$, $b=8/3$, and $r=350$. The frequency of the limit-cycle oscillation is $\omega=16.18$. Figure 12 shows the evolution of ${\bm{X}}_{0}(\theta)=(x_{0},y_{0},z_{0})^{\rm T}$ for one period of oscillation and the corresponding phase sensitivity function ${\bm{Z}}(\theta)=(Z_{x},Z_{y},Z_{z})^{\rm T}$ obtained by the adjoint method ref:ermentrout10 .

We consider two Lorenz models without frequency mismatch and couple them via the coupling matrix $K$ as in Eq. (2). We compare the results for the optimal coupling matrix $K_{\rm opt}$ with those for $K_{I}$. From the results in the previous section, for $P=0.1$, $K_{\rm opt}$ and $K_{I}$ are estimated as

[TABLE]

and

[TABLE]

The linear stability $-\Gamma_{a}^{\prime}(0)$ is approximately $0.872$ for the optimal coupling and $0.365$ for the identity coupling, respectively. Figure 13 shows the antisymmetric parts $\Gamma_{a}(\phi)$ of the phase coupling functions for $K=K_{\rm opt}$ and $K=K_{I}$, and compares the time courses of the phase differences $\phi$ obtained numerically for $\epsilon=0.5$. We can clearly see that the stability of the in-phase state is higher and correspondingly the phase difference decays to zero faster in the optimal case.

It is notable that $K_{\rm opt}$ has several zero components, indicating no feedback from $z$ component to $x$ or $y$ component nor from $x$ or $y$ component to $z$ component arise even after optimization. This is because $z$ component exhibits qualitatively different dynamics from those of $x$ and $y$ components in the Lorenz model. As can be seen from Fig. 12, the fundamental frequency of $z$ component is exactly twice that of $x$ and $y$ components. Reflecting the symmetry of the Lorenz model (invariance under $x\to-x$, $y\to-y$, $z\to z$), the waveforms of $z_{0}(\theta)$ and $Z_{z}(\theta)$ exhibit the same pulse-like oscillations exactly twice while other quantities, $x_{0}(\theta)$, $y_{0}(\theta)$, $Z_{x}(\theta)$, and $Z_{y}(\theta)$, undergo one period of smooth oscillation that is symmetric to $(x,y)\to(-x,-y)$. Therefore, when averaged over one period, feedback from $z$ to $x$ or $y$ (characterized by $Z_{x}$ or $Z_{y}$ multiplied by the difference in $z$ components) vanishes and does not help improve the stability of the synchronized state for the coupled Lorenz oscillators. Similarly, feedback from $x$ or $y$ to $z$ (characterized by $Z_{z}$ multiplied by the difference in $x$ or $y$) does not contribute to the stability.

V Summary and discussion

We have considered a pair of limit-cycle oscillators with weak cross coupling, where different components of the oscillator states are allowed to interact, and optimized the coupling matrix so that the stability of the synchronized state is improved. For oscillators without frequency mismatch, the optimal coupling matrix yields higher linear stability of the in-phase synchronized state. For oscillators with frequency mismatch, a range of phase-locked state with given stationary difference can be realized by choosing the coupling matrix appropriately. Necessary conditions for realizability of a given phase difference are also derived.

In this paper, we have derived the optimal coupling matrix that yields the highest linear stability of the synchronized state for linear diffusive coupling given by Eq. (2). This result can be straightforwardly extended to coupled oscillators with general coupling functions, described by

[TABLE]

where ${\bm{G}}$ represents general nonlinear coupling between the oscillators $1$ and $2$. In this case, the phase coupling function in the reduced phase equations (5) is given by

[TABLE]

instead of Eq. (9). Thus, by defining the function $W(\phi)$ as

[TABLE]

in place of Eq. (13) and calculating $V(\phi)=W(\phi)-W(-\phi)$ and $V^{\prime}(\phi)=W^{\prime}(\phi)+W^{\prime}(-\phi)$ from this $W(\phi)$, the optimization can be performed in a similar way to the linear diffusive case. For example, the optimal coupling matrix for the case without frequency mismatch is given by Eq. (31) with the above $W(\phi)$.

Also, though we have considered only the simple case where all components of the oscillator states can interact with all other components in this paper, it is straightforward to restrict the pairs of components that can actually interact by constraining certain components of $K$ to zero, in order to incorporate realistic physical situations. It would also be interesting to generalize the theory to incorporate different constraints on $K$, for example, to reduce the number of non-zero components by assuming sparsity constraint on $K$.

Although we have considered only the most fundamental two-oscillator problem in this paper, synchronization of a network of many oscillators have attracted much attention Kiss ; Kiss2 ; ref:tanaka08 ; ref:yanagita10 ; ref:yanagita12 ; ref:yanagita14 ; ref:skardal14 ; ref:skardal16 ; ref:nishikawa06a ; ref:nishikawa06b ; ref:nishikawa10 ; Ravoori ; Dorfler ; LiWong ; Stankovski , and generalization of the present framework to many-oscillator networks would be an interesting future problem. For the simplest globally coupled population of $N$ identical oscillators described by

[TABLE]

it is expected that the optimal coupling matrix for the two-oscillator case would also provide faster convergence to global synchrony than the identity coupling matrix. To illustrate this, we simulated $N=400$ SL oscillators with the same parameter values as in Sec. IV, starting from uniformly random initial conditions on the limit cycle. Figure 14 shows synchronization processes for $K=K_{\rm opt}$ and for $K=K_{I}$ in Eqs. (98) and (102), where evolution of the modulus of the Kuramoto order parameter, estimated by $R=\left|(1/N)\sum_{i=1}^{N}\exp[\sqrt{-1}\arctan(y_{i}/x_{i})]\right|$, is plotted. We can observe that the oscillators exhibit much faster convergence to complete synchrony ($R=1$) with $K=K_{\rm opt}$ than with $K=K_{I}$, as expected. Of course, for more complex oscillator networks with frequency heterogeneity and coupling randomness, the result of optimization for the two-oscillator case would not apply due to many-body effects and further investigation should be necessary.

Finally, synchronization between spatiotemporal rhythms in chemical systems has been studied recently ref:mikhailov06 ; ref:mikhailov13 ; ref:mikhailov14 ; hildebrand ; Fukushima ; Epstein , and generalization of the phase reduction theory to reaction-diffusion equations exhibiting stable spatiotemporal oscillations has also been performed Nakao . The present framework can also be extended to such situations and can be used to derive the optimal coupling schemes between two coupled spatiotemporal oscillations. A study in this direction is reported in our forthcoming article kawamura , where improvement in the stability of synchronized states between reaction-diffusion systems by introducing linear spatial filters into mutual coupling is considered.

Acknowledgements.

S.S. acknowledges financial support from JSPS (Japan) KAKENHI Grant Number JP15J12045. Y.K. acknowledges financial support from JSPS (Japan) KAKENHI Grant Number JP16K17769. H.N. acknowledges financial support from JSPS (Japan) KAKENHI Grant Numbers JP16H01538, JP16K13847, and JP17H03279.

Appendix

V.1 Matrix formulas

The tensor product of $m$-dimensional vectors ${\bm{a}}=(a_{1},...,a_{m})^{\rm T}$ and ${\bm{b}}=(b_{1},...,b_{m})^{\rm T}$ gives a $m\times m$ matrix whose $(i,j)$-component is $[{\bm{a}}\otimes{\bm{b}}]_{ij}=a_{i}b_{j}$ for $i,j=1,...,m$. The inner product of $m\times m$ matrices $A$ and $B$ is defined as

[TABLE]

where $A_{ij}$ and $B_{ij}$ represent $(i,j)$-components of the matrices $A$ and $B$, respectively. The Frobenius norm of a $m\times m$ matrix $A$ is defined as

[TABLE]

and the inner product of the matrices $A$ and $B$ is defined as

[TABLE]

Derivative of the inner product of matrices is given by

[TABLE]

and derivative of the Frobenius norm is given by

[TABLE]

For arbitrary matrices $A$ and $B$, the Schwartz inequality

[TABLE]

holds, which can be shown by plugging $\lambda={\rm Tr\ }(AB^{\rm T})/\|B\|^{2}$ into an inequality $\|A-\lambda B\|^{2}\geq 0$ that holds for arbitrary $\lambda$.

V.2 Calculation of $W(\phi)$ and $V(\phi)$

Using ${\bm{Z}}(\theta)$ and ${\bm{X}}_{0}(\theta)$, $W(\phi)$ is explicitly given as

[TABLE]

From this $W(\phi)$, the function $V(\phi)$ can be calculated as

[TABLE]

where $2\pi$-periodicity of the functions ${\bm{Z}}(\theta)$ and ${\bm{X}}_{0}(\theta)$ was used. Similarly, the derivatives of $W(\phi)$ and $W(-\phi)$ can be calculated as

[TABLE]

and

[TABLE]

where $2\pi$-periodicity was used again. Therefore, the derivative of $V(\phi)$ can be calculated as

[TABLE]