Linear response theory for coupled phase oscillators with general coupling functions

TL;DR

This paper introduces a comprehensive linear response theory for coupled phase oscillators that accommodates diverse coupling functions, natural frequency distributions, phase-lags, and time-delays, surpassing previous limited models.

Contribution

The authors develop a general linear response framework applicable to a broad class of coupled oscillators, extending beyond the constraints of prior methods like Ott--Antonsen ansatz.

Findings

Theory accurately predicts oscillator behavior in simulations.

Applicable to systems with arbitrary coupling functions and parameters.

Extends the analytical tools available for studying synchronization.

Abstract

We develop a linear response theory by computing the asymptotic value of the order parameter from the linearized equation of continuity around the nonsynchronized reference state using the Laplace transform in time. The proposed theory is applicable to a wide class of coupled phase oscillator systems and allows for any coupling functions, any natural frequency distributions, any phase-lag parameters, and any values for the time-delay parameter. This generality is in contrast to the limitation of the previous methods of the Ott--Antonsen ansatz and the self-consistent equation for an order parameter, which are restricted to a model family whose coupling function consists of only a single sinusoidal function. The theory is verified by numerical simulations.

Equations

$$\frac{d{\theta }_j}{dt}(t)={\omega }_j+\frac{1}{N}\sum _{k=1}^N\Gamma \left({\theta }_j(t)-{\theta }_k(t-\tau )\right)+H\left({\theta }_j(t),t\right),$$

$$H(\theta ,t)=\Theta (t)\bar{H}(\theta -{\omega }_{\mathrm{ex}}t),$$

$${\int }_{-\infty }^{\infty }d\omega g(\omega )=1.$$

$$\Gamma (\theta )=-\sum _{n>0}{K}_n\sin (n\theta +{\alpha }_n)$$

$$H(\theta ,t)=-\Theta (t)\sum _{n>0}{h}_n\sin [n(\theta -{\omega }_{\mathrm{ex}}t)-{\beta }_n].$$

$$z_m(t)=\frac{1}{N}\sum _{k=1}^N{e}^{im{\theta }_k(t)},(m=1,2,\cdots )$$

$$\chi_{mn}(\omega_{\rm ex})e^{in\omega_{\rm ex}t}=\frac{\partial z_{m}}{\partial h_{n}}\Bigg{|}_{\mathbf{h}\to\mathbf{0}}(t),\quad\mathrm{as}\quad t\to\infty.$$

$$\frac{d{\theta }_j}{dt}(t)={\omega }_j+\sum _{n=-\infty }^{\infty }\left[{\tilde{\Gamma }}_n{z}_{-n}(t-\tau )+{\tilde{H}}_n(t)\right]e^{in{\theta }_j},$$

$$\Gamma (\theta )=\sum _{n=-\infty }^{\infty }{e}^{in\theta }{\tilde{\Gamma }}_n$$

$$H(\theta ,t)=\sum _{n=-\infty }^{\infty }{e}^{in\theta }{\tilde{H}}_n(t).$$

$${\tilde{\Gamma }}_n=\frac{i}{2}{K}_n{e}^{i{\alpha }_n},{\tilde{\Gamma }}_{-n}={\tilde{\Gamma }}_n^{\ast },(n>0),$$

$${\tilde{H}}_n(t)=\frac{i}{2}{h}_n\Theta (t)e^{-i(n{\omega }_{\mathrm{ex}}t+{\beta }_n)},{\tilde{H}}_{-n}(t)={\tilde{H}}_n^{\ast }(t),(n>0),$$

$$\frac{\partial F}{\partial t}+\frac{\partial }{\partial \theta }\left(VF\right)=0,$$

$${\int }_0^{2\pi }d\theta {\int }_{-\infty }^{\infty }d\omega F(\theta ,\omega ,t)=1,$$

$$V(\theta ,\omega ,t)=\omega +\sum _{n=-\infty }^{\infty }{e}^{in\theta }\left[{\tilde{\Gamma }}_n{z}_{-n}(t-\tau )+{\tilde{H}}_n(t)\right].$$

$$z_m(t)={\int }_0^{2\pi }d\theta {\int }_{-\infty }^{\infty }d\omega e^{im\theta }F(\theta ,\omega ,t),(m=1,2,\cdots ).$$

$$F(\theta ,\omega ,t)=F_0(\omega ),\mathrm{for}t<0,$$

$$F_0(\omega )=\frac{g(\omega )}{2\pi }.$$

$$F(\theta ,\omega ,t)=F_0(\omega )+f(\theta ,\omega ,t)\mathrm{for}t\ge 0,$$

$$f(\theta ,\omega ,0)=0.$$

$$z_n(t)={\int }_0^{2\pi }d\theta {\int }_{-\infty }^{\infty }d\omega e^{in\theta }f(\theta ,\omega ,t).$$

$$\frac{\partial f}{\partial t}(\theta ,\omega ,t)+\omega \frac{\partial f}{\partial \theta }(\theta ,\omega ,t)+F_0(\omega )\left[\frac{\partial v}{\partial \theta }(\theta ,t,\tau )+\frac{\partial H}{\partial \theta }(\theta ,t)\right]=0,$$

$$v(\theta ,t,\tau )=\sum _{n=-\infty }^{\infty }{e}^{in\theta }{\tilde{\Gamma }}_n{z}_{-n}(t-\tau ).$$

$$f(\theta ,\omega ,t)=\sum _{n=-\infty }^{\infty }{e}^{in\theta }{\tilde{f}}_n(\omega ,t),$$

$$\frac{\partial {\tilde{f}}_{-n}}{\partial t}(\omega ,t)-in\omega {\tilde{f}}_{-n}(\omega ,t)-in{F}_0(\omega )\left[{\tilde{\Gamma }}_{-n}{z}_n(t-\tau )+{\tilde{H}}_{-n}(t)\right]=0.$$

$$\hat{\phi }(s)={\int }_0^{\infty }dt\phi (t)e^{-st}\mathrm{for}\mathrm{Re}(s)>0,$$

$$(s-in\omega ){\hat{f}}_{-n}(\omega ,s)=in{F}_0(\omega )\left[{\tilde{\Gamma }}_{-n}{e}^{-s\tau }{\hat{z}}_n(s)+{\hat{H}}_{-n}(s)\right]+{\tilde{f}}_{-n}(\omega ,0),$$

$${\hat{z}}_n(s)=2\pi {\int }_{-\infty }^{\infty }d\omega {\hat{f}}_{-n}(\omega ,s).$$

$${\hat{z}}_n(s)=in\left[{\tilde{\Gamma }}_{-n}{e}^{-s\tau }{\hat{z}}_n(s)+{\hat{H}}_{-n}(s)\right]{\int }_{-\infty }^{\infty }d\omega \frac{g(\omega )}{s-in\omega }+2\pi {\int }_{-\infty }^{\infty }d\omega \frac{{\tilde{f}}_{-n}(\omega ,0)}{s-in\omega }.$$

$$I_n(s)=in{\int }_{-\infty }^{\infty }d\omega \frac{g(\omega )}{s-in\omega },$$

Full text

1 Introduction

Synchronization among rhythmic elements can be ubiquitously observed in a wide variety of systems, whose numbers of elements range from small to large [1, 2, 3, 4]. Since synchronization between two pendulum clocks hanging on the same wall was observed by Huygens, it has also been found in many large systems, such as flashing fireflies [1], cardiac cells [5], the circadian rhythm in mammals [2], and neuronal populations [3, 6]. Synchronization is often essential to biological functions, in large systems in particular, and we focus on large systems in this study.

Based on mathematical modeling of a rhythmic element and the phase reduction theory [7, 8, 9], synchronization is described by coupled phase-oscillator models. For example, the Kuramoto model [7], a paradigmatic coupled phase-oscillator model, has two types of states: a nonsynchronized state, in which each oscillator rotates with its proper frequency called the natural frequency, and partially synchronized states, in which part of the oscillators rotate with the same effective frequency. Strengthening the couplings provides a synchronization transition from the nonsynchronized state to the partially synchronized states, and the continuity of the transition is determined by the natural frequency distribution in the Kuramoto model [7, 10, 11, 12, 13]. We can thus reproduce macroscopic dynamics such as the continuous and discontinuous synchronization transitions by controlling the microscopic details that include the coupling strength and the natural frequencies.

However, the single direction from a microscopic model to macroscopic dynamics does not provide a complete picture of the synchronization in real systems because we have often no prior knowledge of the microscopic character, such as the coupling strength or natural frequency distribution. Additionally, it is difficult to access directly the microscopic character. This motivates us to develop a theory that extracts the microscopic character from a macroscopic experiment. One potential candidate is the linear response theory; a clear example of an application is spectroscopy.

One strategy for obtaining the linear response formula in coupled oscillator models is to construct and analyze the self-consistent equation for the order parameter based on knowledge of the stationary states under an external force. This strategy has been developed in models that have a single sinusoidal coupling function[14, 15] and can reach the nonlinear regime beyond the linear response. However, the construction of the self-consistent equation is not easily extended to general systems whose coupling functions consist of many harmonics of sinusoidal functions because there are several stable stationary states for a given set of parameters [16, 17, 18].

Another strategy relies on the Ott–Antonsen ansatz [19, 20]. This ansatz reduces the equation of continuity, which describes the dynamics in the large population limit, to finite-dimensional ordinary differential equations. See ref. [15] for an application to the linear response. However, the benefit of this reduction is limited in the single harmonic case again, and the natural frequency distribution must be rational, e.g., a Lorentzian, so that one can apply the residue theorem.

The above two strategies result in difficulties when we consider a general coupling function, while coupling functions often consist of several harmonics in neuronal networks [21], in electrochemical oscillators [22, 23], and near a Hopf bifurcation [24, 25]. In addition to the generality of coupling functions, a desired linear response theory must provide susceptibility for all rotation frequencies of external forces from the point of view of experimental measurements. Furthermore, the time delay in couplings, which drastically changes its synchrony [26, 27], must be incorporated because it is inevitable in many natural systems. Thus, we propose a linear response theory that is applicable to systems with any coupling functions accompanied by a time delay, any natural frequency distributions, and any rotational frequencies of external forces.

Linear response theory has been developed in statistical mechanics [28, 29], but computing a value of susceptibility is not easy because we must solve the equations of motion for an $N$-body system, which is nonintegrable in general. In the class of globally coupled interactions, this difficulty has been overcome in the large population limit by using the Vlasov equation, which describes the dynamics of a one-body distribution function and is an analogy to the equation of continuity [30, 31]. Our approach has been inspired by this linear response theory.

The construction of this paper is as follows. In Section 2, we describe the model which we analyze in this study. We describe the constructed linear response theory based on the equation of continuity with a general coupling function and a time delay in Section 3. In Section 4, we describe the numerical simulations used to validate the theory. Section 5 is devoted to the summary and discussion.

2 Model

In this section, we introduce our model, relevant quantities, and the large population limit.

2.1 Settings

We consider the system described by the ordinary differential equations

[TABLE]

where $\theta_{j}(t)\in[0,2\pi)$ is the phase of the $j$th phase-oscillator, $\omega_{j}$ is its time-independent natural frequency, $\Gamma(\theta_{j}(t)-\theta_{k}(t-\tau))$ is the coupling function from oscillator $k$ to $j$, $\tau$ refers to the time delay of the couplings, and $H(\theta,t)$ is the external force. We assume that function $H(\theta,t)$ has the form

[TABLE]

where $\omega_{\rm ex}$ is the external frequency and $\Theta(t)$ is the Heviside step function. The external force is turned off for $t<0$ and is turned on at $t=0$. To investigate the linear response, we assume that the external force is sufficiently small, $|H(\theta,t)|\ll 1$.

The natural frequencies follow the natural frequency distribution $g(\omega)$, which satisfies the normalization condition of

[TABLE]

We note that, for a nonzero time delay of $\tau>0$, system (1) is not invariant under the rotating reference frame with frequency $\omega_{0}$ and a shift in $\bar{g}(\omega)=g(\omega-\omega_{0})$, whereas this invariance holds for $\tau=0$.

Referring to the works by Kuramoto [33] and Daido [34], we expand the coupling function $\Gamma$ and the external force $H$ into the Fourier series as

[TABLE]

and, from the form of the external force (2),

[TABLE]

Real constants $K_{n}$, $\alpha_{n}$, and $h_{n}$ represent the coupling strength, the phase-lag, and the amplitude of the external force for the Fourier mode $n$, respectively. The real parameter $\beta_{n}$ determines the direction of the external force in the rotating reference frame with frequency $\omega_{\rm ex}$.

2.2 Order parameters and susceptibilities

To study the synchrony of the oscillators, we introduce the order parameters

[TABLE]

which are called the Daido order parameters [34] and detect the clusters whose phases are congruent modulo $2\pi/m$. The first-order parameter, $z_{1}$, is equivalent to the Kuramoto order parameter [33]. Owing to the rotating external force $H(\theta,t)$, the order parameters also rotate with external frequency $\omega_{\rm ex}$. Accordingly, the susceptibility tensor $\boldsymbol{\chi}(\omega_{\rm ex})=(\chi_{mn}(\omega_{\rm ex}))$ depending on the external frequency is asymptotically defined by

[TABLE]

In other words, $\chi_{mn}(\omega_{\rm ex})$ is defined as the rate of change in the $m$th Daido order parameter induced by the $n$th Fourier mode of the external force in the rotating reference frame with frequency $\omega_{\rm ex}$. We aim to obtain the susceptibility $\chi_{mn}(\omega_{\rm ex})$ by taking the large population limit in a coupled oscillator system that has a general coupling function.

2.3 Large population limit

Using the Daido order parameters (6), we have a simple expression for the equations of motion (1) as

[TABLE]

where $\widetilde{\Gamma}_{n}$ and $\widetilde{H}_{n}(t)$ are the Fourier components of the coupling function $\Gamma(\theta)$ and the external force $H(\theta,t)$, which are defined by

[TABLE]

and

[TABLE]

The present forms of $\Gamma$ (4) and $H$ (5) give

[TABLE]

and

[TABLE]

where $\widetilde{\Gamma}_{n}^{\ast}$ is the complex conjugate of $\widetilde{\Gamma}_{n}$, for instance, and we define $\widetilde{\Gamma}_{0}=\widetilde{H}_{0}=0$.

In the large population limit of $N\to\infty$, the dynamics of the system (8) are described by the equation of continuity

[TABLE]

where $F(\theta,\omega,t)$ is the probability density function that satisfies the normalization condition

[TABLE]

and velocity $V$ is given by

[TABLE]

The Daido order parameters are expressed as

[TABLE]

3 Linear response theory

We derive the susceptibility in the nonsynchronized state by solving the linearized equation of continuity around the nonsynchronized state. The solution is obtained by performing the Fourier transform in the phase variable and the Laplace transform in the time variable.

3.1 Linearization around the nonsynchronized state

As the reference state for $t<0$, we consider the nonsynchronized state

[TABLE]

which is the trivial stationary solution to the equation of continuity (13). The normalization condition (14) gives

[TABLE]

The small external force $H(\theta,t)$ that starts at $t=0$ perturbs the state from $F_{0}$ to

[TABLE]

where $|f(\theta,\omega,t)|\ll 1$ is assumed according to $|H(\theta,t)|\ll 1$. The continuity of $F(\theta,\omega,t)$ at $t=0$ implies that the initial condition of $f(\theta,\omega,t)$ is

[TABLE]

The Daido order parameters receive no contribution from the nonsynchronized reference state, $F_{0}(\omega)$, and they have the same order with perturbation $f$ as

[TABLE]

This ordering produces the linearized equation of continuity as

[TABLE]

where we defined

[TABLE]

The linearized equation (22) is valid under a weak external force.

3.2 Fourier–Laplace transforms

To solve the linearized equation (22), we employ the Fourier and Laplace analyses. Substituting the Fourier series expansion of $f$,

[TABLE]

and of $H$ (10) into the linearized equation (22), we find the Fourier-expanded expression for each $n$ as

[TABLE]

Further, we perform the Laplace transform which is defined for $\varphi(t)$ to be

[TABLE]

where the condition ${\rm Re}(s)>0$ is assumed to ensure the convergence of the integral. The Laplace transform modifies equation (25) as

[TABLE]

where we used $z_{n}(t)=0$ for $t<0$. We note that the last term disappears owing to the initial condition (20); however this term is kept tentatively to discuss the stability of the reference state, $F_{0}(\omega)$.

Both sides of Eq. (27) contain $\widehat{f}_{-n}$, and $\widehat{f}_{-n}$ is associated with the Laplace transform of the order parameter, $\widehat{z}_{n}(s)$, as

[TABLE]

Using this relation, we solve Eq. (27) with respect to $\widehat{z}_{n}(s)$. By multiplying (27) by $2\pi(s-in\omega)^{-1}$ and integrating over $\omega$, we have

[TABLE]

By defining the functions

[TABLE]

and

[TABLE]

the Laplace transform of the order parameter $z_{n}(t)$ is expressed as

[TABLE]

The temporal evolution of $z_{n}(t)$ is obtained by performing the inverse Laplace transform.

The inverse Laplace transform of function $\widehat{\varphi}(s)$ is defined by

[TABLE]

where the integral is performed along the Bromwich contour which lies on the right-hand side of any singularities of $\widehat{\varphi}(s)$ on the complex $s$ plane. If $\widehat{\varphi}(s)$ has a simple pole at $s=s_{0}\in\mathbb{C}$, the inverse Laplace transform produces a mode proportional to $\exp(s_{0}t)$ upon picking up this pole in the residue theorem. As a result, a pole of ${\rm Re}(s_{0})<0$ corresponds to a stable mode, and a pole of ${\rm Re}(s_{0})>0$ corresponds to an unstable mode.

We note that integral $I_{n}(s)$, and function $G_{n}(s)$ are defined on the right-half plane of ${\rm Re}(s)>0$ of the complex $s$ plane following the definition of the Laplace transform (26). To apply the discussion above to $\widehat{z}_{n}(s)$, we need to perform the analytic continuation of $I_{n}(s)$. See A for this continuation.

The stability of the nonsynchronized state $F_{0}(\omega)$ is examined by turning off the external force, $H(\theta,t)=0$, and is determined by the roots of $G_{n}(s)$. A demonstration of the stability analysis is provided in B with the use of the Nyquist diagrams. Assuming the stability of the nonsynchronized state $F_{0}(\omega)$ through the rest of this study, we compute the susceptibility from (32), whose last term vanishes because of the initial condition (20).

3.3 Susceptibilities

The susceptibility $\chi_{mn}(\omega_{\rm ex})$ is defined in the limit of $t\to\infty$ as in (7), and the asymptotic temporal evolution of $z_{n}(t)$ is obtained by considering the singularities of (32). The stability assumption of $F_{0}$ implies that all the roots of $G_{n}(s)$ are on the left-half plane ${\rm Re}(s)<0$. Consequently, the asymptotic behavior of the order parameters is determined by the residues of the poles on the imaginary axis, which come from the external force as, from (12),

[TABLE]

The residue theorem gives

[TABLE]

and

[TABLE]

where, for $n>0$,

[TABLE]

and

[TABLE]

PV represents the Cauchy principal value. ¿From the definition (7), the susceptibility tensor $\boldsymbol{\chi}(\omega_{\rm ex})=(\chi_{mn}(\omega_{\rm ex}))$ is diagonal as

[TABLE]

and, for $n>0$,

[TABLE]

and

[TABLE]

The susceptibility (40) is the main result of this study. This main result (40) is an extension of the susceptibility [14] obtained in the Kuramoto model [7, 10], whose coupling function is $\Gamma(\theta)=-\sin\theta$ without the time-delay. We provide three remarks derived from this main result (40).

The first remark is on the divergence of the susceptibility. The divergence appears if the complex denominator of the susceptibility (40) vanishes. For instance, the divergence is observed at $K_{n,{\rm c}}=2/[\pi g(0)]$ if the model has three types of symmetry: a zero phase lag $\alpha_{n}=0$, a zero time delay $\tau=0$, and a symmetric natural frequency $g(\omega)$, where the satisfaction of the last symmetry gives $J(0)=0$. In contrast, the divergence does not appear for any value of $K_{n}$ unless we tune the external frequency $\omega_{\rm ex}$ so that the imaginary part of the denominator vanishes, if one of the three types of symmetry is broken. The absence of divergence is demonstrated in C in the Kuramoto model with an asymmetric natural frequency. In the next section, we discuss our examination of the theory with the tuned $\omega_{\rm ex}$ to observe the divergence because the divergence is a representative phenomenon at the critical point, and it provides a deep numerical examination: the divergence expands the discrepancy of susceptibility between the theoretically predicted value and the numerically observed value that is perturbed by breaking the assumption of the large population limit or the zero external force limit.

The second remark is on the critical exponent $\gamma$, which is defined as

[TABLE]

around the critical coupling strength $K_{\rm c}$ [32]. We choose $\omega_{\rm ex}$ to observe the divergence. The denominator of the susceptibility (40) depends on the coupling constant linearly, and the critical exponent $\gamma$ should be unity: $\gamma=1$. This universality of $\gamma$ gives a sharp contrast with the dependence of the critical exponent $\beta$ on the coupling function: the single harmonic case has $\beta=1/2$ [33] while $\beta=1$ for the general case [35, 36]. The critical exponent $\beta$ is defined here, with absence of the external force, by

[TABLE]

above and around the critical point.

The last remark is regarding the roles of the phase lag and the time delay. The phase lag $\alpha_{n}$ and the time delay $\tau$ are included only in the exponential factor of the denominator and play a similar role in the susceptibility (40). In fact, the phase lag can be eliminated from the time-delayed Kuramoto model if the external frequency, $\omega_{\rm ex}$, is fixed [26]. However, the time delay couples to $\omega_{\rm ex}$, and, from experimental observation of the susceptibility, we can identify which factor is included in the system by varying $\omega_{\rm ex}$.

4 Numerical tests

To verify our theoretical results, we perform numerical simulations, where we also shed light on the difference between the roles of the phase lag $\alpha_{n}$ and the time delay $\tau$. We first examine the linear response in the Sakaguchi–Kuramoto model, which contains the phase lag but not the time delay in a single sinusoidal coupling function. The single harmonics permits us to apply the Ott–Antonsen reduction, which eliminates finite-size fluctuations and is useful for examining the theory. The second model is the time-delayed Daido–Kuramoto model, which has multi-harmonics in the coupling functions. This model highlights the usefulness of the proposed theory for general systems. In the above two models, we use the Lorentz distribution

[TABLE]

as the natural frequency distribution, although the proposed theory can be applied to other distributions. Note that the time delay breaks the invariance of the system with respect to the Galilei transform of $\theta$ and hence the rotation of the reference frame and $\Omega$ cannot be removed by shifting $\omega$ if the coupling function has a nonzero time delay, $\tau$. The principal value, $J(\omega_{\rm ex})$ (38), can be explicitly computed for the Lorentzian (44), as shown in D.

We adopt the 4th-order Runge–Kutta method to integrate the original equations (1) and the reduced equations, where the time length is $T=1000$ with a discrete time size $dt=0.05$. Values of the order parameters are evaluated by averaging over the time interval of $[T/2,T]$.

4.1 The Sakaguchi–Kuramoto model

The Sakaguchi–Kuramoto model is recovered by setting the coupling function, $\Gamma(\theta)$, to be

[TABLE]

without the time-delay, where $0<\alpha<\pi$. The first susceptibility is read as

[TABLE]

We choose the frequency of the external force, $\omega_{\rm ex}$, so that the system has the divergence of $\chi_{1}$ at the critical point. The divergence appears when the denominator of (46) vanishes. The imaginary part provides the condition for $\omega_{\rm ex}$ as

[TABLE]

and the real part determines the critical point, $K_{\rm c}$, as

[TABLE]

For a fixed phase-lag parameter $\alpha=1$, susceptibility $\chi_{1}$ obtained theoretically, (46), is shown in Figure 1 with numerical results in the $N$-body system and in the reduced system derived by the Ott–Antonsen ansatz. The reduced system, corresponding to the limit $N\to\infty$, is in good agreement with the theory for a sufficiently small external force, $h=10^{-4}$, whereas the agreement is not perfect for $h=10^{-2}$ which is not sufficiently small for imitating the limit of $h\to 0$ taken in the definition of susceptibility, (7). This observation concludes that the discrepancy between the theory and numerics results from the nonlinearity of the response with respect to $h$ rather than the finite-size fluctuation in the $N$-body system. We have two remarks in order. First, to observe the linear response clearly, $h$ in the $N$-body system must be larger than $O(1/\sqrt{N})$, which is the expected level of the finite-size fluctuation of the order parameter. Thus, we must use a larger $N$ to use a smaller $h$. Second, the nonlinearity inhibits the divergence at the critical point because the divergence results from the linear response. The nonlinearity, therefore, enhances the discrepancy around the critical point as observed for $h=10^{-2}$ in Figure 1.

4.2 The time-delayed Daido–Kuramoto model

The time-delayed Daido–Kuramoto model has the coupling function of

[TABLE]

with $M>1$ and the nonzero time-delay parameter $\tau$ in (1). The susceptibilities are explicitly written as

[TABLE]

This susceptibility for $n=1$ is obtained by replacing $\alpha$ with $\omega_{\rm ex}\tau$ in the susceptibility (46) of the Sakaguchi–Kuramoto model. We adopt $M=2$ in the numerical simulations for simplicity, although the theory is applicable to any arbitrary $M$. The coupling constant of the second Fourier mode, $K_{2}$, will be sufficiently small positive number so that the nonsynchronized state becomes unstable in the first Fourier mode. We set the time-delay parameter as $\tau=2$ here and choose the value of $\omega_{\rm ex}$ to satisfy the divergence condition.

It is worth commenting on the continuity of the synchronization transition. First, the theory is applicable in the nonsynchronized state even if the transition is discontinuous. One advantage of the continuous transition is that numerical examinations can be robustly performed, while the discontinuous transition, owing to the bistability, may accompany a jump of responses around the critical point for a finite external force. The parameter set, $\gamma=0.1,\Omega=3$, and $\tau=2$, gives the continuous transition in the delayed Kuramoto model [26], but it is not obvious if the time-delayed Daido–Kuramoto model also exhibits the continuous transition. For instance, the multi-harmonic coupling function with $M=2$ results in the discontinuous transition if $K_{2}>0$ and $\tau=0$ [37]. Our numerical computation, however, implies that the transition is continuous with the given parameter set.

The Ott–Antonsen reduction is not applicable to the Daido–Kuramoto model, and we show the results for only $N$-body simulations in Figure 2. The numerical results do not perfectly agree with the theoretical prediction, but they approach the theoretical curves as the external force becomes small. This tendency of discrepancies is very similar to the ones observed in the Sakaguchi–Kuramoto model, in which the validity of the theory has been confirmed with the aid of the reduced system. Thus, we conclude that the theory is also valid in the time-delayed Daido–Kuramoto model.

4.3 Distinction between the phase lag and the time delay

Finally, we illustrate that phase lag $\alpha$ plays the same role as time delay $\tau$ for a fixed $\omega_{\rm ex}$ but a different role for a varying $\omega_{\rm ex}$. For simplicity, we concentrate on the time-delayed Sakaguchi–Kuramoto model and assume that the coupling strength $K_{1}$ and the natural frequency distribution, $g(\omega)$, are known. This assumption implies that, in a rewritten form of (40),

[TABLE]

the first term of the most right-hand side is known, but the second term must be pointwisely observed in an experiment with inevitable errors. We imitate the errors by adding noise $\xi^{(l)}(\omega_{\rm ex})$ to the theoretical value of $\chi_{1}(\omega_{\rm ex})$ at each observing point $\omega_{\rm ex}$: $\chi_{1}^{(l)}(\omega_{\rm ex})=\chi_{1}(\omega_{\rm ex})+\sigma\xi^{(l)}(\omega_{\rm ex})$. The random value, $\xi^{(l)}(\omega_{\rm ex})$, is independently drawn from the standard normal distribution, and superscript $l$ denotes the index of realization, in other words, virtual observation. In Figure 3, we show the means and standard deviations of the right-hand side of (51) averaged over $100$ realizations of $\chi_{1}^{(l)}(\omega_{\rm ex})$ for $\sigma=10^{-3}$ and $10^{-2}$. For $(\tau,\alpha)=(2,0)$, the real and imaginary parts of the right-hand side of (51) are wavy with period $\pi$ according to the choice of $\tau$, whereas no wave is found for $(\tau,\alpha)=(0,2)$, as we expect from (51). The large standard deviations for large $|\omega_{\rm ex}-\Omega|$ are due to the lack of population of oscillators, that is the smallness of $g(\omega)$. This demonstration suggests that a sufficiently precise observation of the susceptibility can identify the existence and value of the time delay. We must remember that knowledge of the coupling strength and the natural frequency is assumed in the above demonstration; however, it is usually difficult to access these in advance. Identification of the coupling strength and the natural frequency distribution is out of the scope of this study but should be investigated.

5 Summary and discussion

We have developed a linear response theory for coupled oscillator systems by directly solving the linearized equation of continuity around the nonsynchronized state. One considerable advantage of the proposed theory is that it is applicable to systems that have general coupling functions, phase lags, and a time delay, while the previous two methods of the self-consistent equation and the Ott–Antonsen reduction are restricted to the systems that have a single sinusoidal coupling function. The theoretical predictions have been successfully verified by numerical simulations.

The reference state is assumed to be nonsynchronized as the first step for a general theory. Another interesting topic would be to construct a linear response theory for other types of reference states: partially synchronized states, cluster states [38], chimera states [39, 40], glassy states [41], chaos states [42], and so on.

Finally, we comment on a possible application of the linear response theory for the identification problem. This linear response theory gives macroscopic responses from microscopic details, and the identification problem can be formulated as the inverse problem. A simple identification using the proposed theory was demonstrated in Section 4.2 to distinguish between the phase lag and the time delay by surveying the dependence on the external frequency wtih knowledge of the coupling strength and the natural frequency distribution. A more systematic and precise identification method will be discussed elsewhere.

This work was supported by the Special Postdoctoral Research Program at RIKEN (Y.T.) and JSPS KAKENHI Grand Numbers 19K20365 (Y.T.) and 16K05472 (Y.Y.Y.).

Appendix A Analytic continuation

Function $I_{n}(s)$ is first defined in ${\rm Re}(s)>0$ in (30), which is the domain of the Laplace transform (26). We continue this function into the entire complex $s$ plane, which is necessary to obtain $I_{n}(in\omega_{\rm ex})$ included in the susceptibility $\chi_{n}(\omega_{\rm ex})$. We have $I_{0}(s)=0$ and assume $n\neq 0$.

In the definition of $I_{n}(s)$ (30), the integral contour with respect to $\omega$ is the real axis and the contour does not meet the singular point $\omega=-is/n$ because $s$ is restricted in the right-half plane, ${\rm Re}(s)>0$. However, in the limit of ${\rm Re}(s)\to+0$, the pole arrives on the real axis from the lower (upper) side of the complex $s$ plane for $n>0$ ($n<0$). To avoid this pole, we smoothly modify the integral contour to the upper (lower) side and continue this modification for ${\rm Re}(s)<0$ so that we obtain the continued function. This continuation gives the explicit form of function $I_{n}(s)$ as

[TABLE]

where the second term for ${\rm Re}(s)\leq 0$ due to the residue at the pole $\omega=-is/n$ and the positive (negative) sign corresponds to $n>0$ ($n<0$).

Appendix B Stability analysis using the Nyquist diagram

The stability of the nonsynchronized state $F_{0}(\omega)$ is analyzed by turning off the external force, $H(\theta,t)=0$. Roughly speaking, the initial perturbation, $f(\theta,\omega,0)$, plays the role of the external force in the response formula (32), and the stability of $F_{0}(\omega)$ is determined by the zero points of $G_{n}(s)$, which is analytically continued via the procedure in A. If there is a zero point whose real part is positive, then $F_{0}(\omega)$ is unstable.

We suppose that $G_{n}(s)$ is a mapping from the complex $s$ plane to the complex $G_{n}$ plane. We focus on the boundary of the stability, the imaginary axis $s=iy~{}(y\in\mathbb{R})$. In the limit of $y\to\pm\infty$ we have $G_{n}\to 1$. Thus, the mapped imaginary axis forms an oriented closed curve and the boundary of the unstable region ${\rm Re}(s)$ can be identified on the complex $G_{n}$ plane by the closed curve and its orientation. This consideration implies that the nonsynchronized state, $F_{0}(\omega)$, is unstable if the mapped unstable region, which is the inside of the closed curve, contains the origin. The Nyquist diagrams for $G_{1}(s)$ are shown in Figure 4 for the Sakaguchi–Kuramoto model and the time-delayed Daido–Kuramoto model.

Appendix C The Kuramoto model with an asymmetric natural frequency distribution

We consider the Kuramoto model by setting the coupling function $\Gamma(\theta)$ as

[TABLE]

and the parameters as $\alpha=\tau=0$ in (4). The susceptibility for the order parameter, $z_{1}$, is written as

[TABLE]

which is obtained through the general expression (40).

We employ an asymmetric natural frequency distribution, $g(\omega)$, to observe the non-divergence of the susceptibility at the critical point for the static external force and show the qualitative difference from the other two types of asymmetry, the phase lag and the time delay. The asymmetry is achieved by the family of $g(\omega)$ as

[TABLE]

where $\Omega\geq 0,\gamma_{1},\gamma_{2}>0$ and the normalization constant $c$ is given by

[TABLE]

as in ref. [43]. The distribution is symmetric if $\gamma_{1}=1$ or $\Omega=0$ and tends to be bimodal with large $\Omega$. We use the parameter set $(\gamma_{1},\gamma_{2},\Omega)=(0.6,1,0.6)$ to realize the asymmetric unimodal, $g(\omega)$.

First, we consider the static external force case with $\omega_{\rm ex}=0$. The susceptibility does not diverge even at the critical point, $K_{\rm c}$, as shown in Figure 5. This nondivergence is in good agreement with the first remark mentioned in Section 3.3.

The divergence recovers if we choose the external frequency, $\omega_{\rm ex}$, appropriately. The divergence appears under the condition $J(\omega_{\rm ex})=0$ at the point of

[TABLE]

This theoretical prediction is confirmed numerically in Figure 6, where we can observe their discrepancy with larger $h$ because of the nonlinearity of the response as explained in Section 4.1.

In the Kuramoto model, the divergence of the susceptibility appears only in the real part, which is shown in Eq. (54). In other words, the response must be parallel to the direction of the external force. In contrast, in the Sakaguichi–Kuramoto model and the time-delayed Daido–Kuramoto model, the phase lag and the time delay permit divergences in both real and imaginary parts of the susceptibility, and the direction of the response is not always parallel to the external force, i.e., we may observe a phase-gap between the external force and the linear response. The coexistence of the divergence and the phase gap thus reveals the different roles of the two groups of asymmetry; asymmetry by the phase lag or the time delay permits the coexistence, while asymmetry in $g(\omega)$ does not.

Appendix D Calculation of susceptibility with Lorentzian natural frequency distributions

In our numerical simulations, we use a Lorentzian or its product as the natural frequency distributions in the main text. In that case the principal value, $J(\omega_{\rm ex})$ (38), can be computed explicitly by using the residue theorem. We introduce the integral

[TABLE]

The pole is at $\omega=\omega_{\rm ex}-i\epsilon$ and the integral is well defined for $\epsilon\neq 0$. We assume $\epsilon>0$. Moreover, we consider the single Lorentzian (44) for simplicity. Adding the upper half-circle contour and using the residue theorem by picking up the pole $\omega=\Omega+i\gamma$, we are able to compute the value of $J_{\epsilon}(\omega_{\rm ex})$ explicitly as

[TABLE]

Taking the limit $\epsilon\to+0$ in (59) and applying the continuation technique presented in A to the right-hand side of (58), we have

[TABLE]

and

[TABLE]

The above idea is applicable to the multiplicative Lorentzian (55).

Appendix E Ott–Antonsen reduction in the case with Lorentzian natural frequency distributions

In the numerical confirmation of the Sakaguchi–Kuramoto model, we used the Ott–Antonsen ansatz [19, 20], which reduces the equation of continuity to a real two-dimensional system describing the dynamics of the order parameter $z_{1}$. This reduction is limited to a class of models, but the reduced equations are useful for examinining our theory because the reduced system corresponds exactly to the large population limit and has no finite-size fluctuations. Here, we derive the reduced equations in the Sakaguchi–Kuramoto model with the natural frequency (55).

The Ott–Antonsen ansatz introduces the form of the probability density function $F$ as

[TABLE]

where the complex-valued function $a(\omega,t)$ satisfies the condition $|a\left(\omega,t\right)|<1$ and is regular on the complex $\omega$-plane. Substituting the ansatz (62) into the equation of continuity for the Kuramoto model with the first harmonic external force having the frequency of $\omega_{\mathrm{ex}}$, we obtain the equation for $a\left(\omega,t\right)$ as

[TABLE]

where the order parameter $z_{1}(t)$ and $a(\omega,t)$ are related through

[TABLE]

which is obtained from (62). The integral over $\omega$ in the right-hand side of (64) is performed by adding the large upper-half circle, which has no contribution to the integral, and picking up the two poles of $g(\omega)$, (55), at $\omega=\Omega+i\gamma_{1}$ and $\omega=-\Omega+i\gamma_{2}$. The residues give

[TABLE]

where complex variables $A$ and $B$ are defined by

[TABLE]

and the time-independent coefficients are given by

[TABLE]

Finally, in (63), by setting $\omega$ as $\omega=\Omega+i\gamma_{1}$ or $\omega=-\Omega+i\gamma_{2}$, we have the reduced equations

[TABLE]

Similarly, for the Sakaguchi–Kuramoto model with the Lorentzian $g(\omega)$, which is discussed in Section 4.1, we have the reduced equation for $z_{1}$:

[TABLE]

It is worth noting that the Ott–Antonsen reduction is also applicable to the time-delayed Kuramoto model [19].