Synchronization limit of weakly forced nonlinear oscillators

Nonlinear oscillators exhibit synchronization (injection-locking) to external periodic forcings, which underlies the mutual synchronization in networks of nonlinear oscillators. Despite its history of synchronization and the practical importance of injection-locking to date, there are many important open problems of an efficient injection-locking for a given oscillator. In this work, I elucidate a hidden mechanism governing the synchronization limit under weak forcings, which is related to a widely known inequality; Hölderʼs inequality. This mechanism enables us to understand how and why the efficient injection-locking is realized; a general theory of synchronization limit is constructed where the maximization of the synchronization range or the stability of synchronization for general forcings including pulse trains, and a fundamental limit of general m : n phase locking, are clarified systematically. These synchronization limits and their utility are systematically verified in the Hodgkin–Huxley neuron model as an example.


Introduction
Entrainment, which adjusts the frequencies of oscillators to that of an external forcing (signal) above a critical forcing amplitude, is a fundamental phenomenon of wide interest with a long history and a large variety of applications [1][2][3][4][5]. Generally, the ratio of these two frequencies is locked to m : n in the entrainment, which is called m : n frequency locking. Likewise, if the oscillation phase of the oscillator (= θ) and that of the external forcing (= θ ext ) satisfy n θ θ − m ext = const, it is called m : n phase locking. And, 1 : 1 phase locking is termed synchronization [4]. As opposed to conventional single-oscillator entrainment, entrainment has significance even in coupled oscillatory elements, thanks to recent studies on networks of coupled oscillators [6][7][8][9]. Simultaneously, in many branches of science and engineering, engineering entrainment (injection-locking) with an efficient forcing (injection signal) has become more important, and methods for efficient entrainment have been developed in recent years [10][11][12][13][14][15][16][17], reflecting advances in experimental techniques for observing and manipulating such oscillators (see [16][17][18][19][20] for instance).
Despite its history, and the advances made in and the wide utility of synchronization to date, it has been open problems of practical and theoretical significance to realize the synchronization limit. Part of our results realizing the synchronization limit is schematically illustrated as in figure 1, which solve the three basic open problems, posed in P1, P2, and P3 below.
P1: What sort of mechanism determines the best (i.e., global optimal) power-reduced (average square of its waveform is bounded) forcing that produces the maximum entrainment range or the most stable synchronization? And are the optimal power-reduced forcings obtained in [16,17] really the best ones if a more general class (function space) of forcings is considered? Figure 1. Best forcings revealed by the mechanism of synchronization limits when maximizing the range of entrainment: (a) power-reduced forcing in P1, (b) areareduced forcing in P2, and (c) magnitude-reduced forcing in P2, obtained for (d) a given phase response function Z (taken from the Hodgkin-Huxley phase model [17] as an example). f p opt, ( = ∞ p 1, 2, and ) here is determined algorithmically in the main text.
P2: In conjunction with P1, what determines the best area-reduced (total absolute value of its waveform is bounded) forcing or magnitude-reduced (its maximum amplitude is bounded) forcing?
P3: Do the answers to P1 and P2 regarding 1 : 1 phase locking remain valid for general m : n phase locking? If so, how is the best forcing for m : n phase locking related to that for 1 : 1 phase locking? Or, if not, what limits an ideal (most efficient) forcing for m : n phase locking?
Recent studies have approached these problems, from different formulations using calculus of variations ( [16] for optimizing the locking range, [17] for optimizing the stability of phase locking, and [21] for both of them in m : n phase locking, for a power-reduced case). However, regarding P1, these conventional studies only suggest the existence of a synchronization limit by finding (possibly local optimal) forcings in each particular situation, and regarding P2 and P3, methods used in these conventional studies are not applicable, and hence a fundamental limit of entrainability has not been clarified. Toward this end, in this paper, we find an underlying mechanism in the above three problems leads us to a unified, global view of synchronization limits and their constructions.

Entrainment modeled by the phase equation
Here we introduce the well-known phase equation for the weakly forced nonlinear oscillators [1][2][3][4][5]. The entrainment process of a limit-cycle oscillator in the weak forcing limit can be modeled by ψ ω ϵ ψ Ω ), Z is the phase response (sensitivity) function, and ω and Ω are the natural frequency of the oscillator and the frequency of the weak forcing ϵ Ω f t ( ), respectively, following the notation in [2]. In general, m : n phase locking occurs when ≈ ω Ω m n is satisfied for positive relatively prime integers m and n. In this situation, the above equation is further simplified by the method of averaging (after setting ϵ to 1) to the following phase equation is determined by f and Z as 1 2 , in which T = π Ω n 2 (n times the natural period of the oscillator), θ ∈ π π − [ , ] represents Ωt n , and 〈 〉 · denotes the integration over its period π 2 : · d . When considering the case of = = m n 1, i.e., 1 : 1 entrainment, we will abbreviate Γ m n as Γ, for simplicity. Now, we define the synchronizability of the oscillator by equation (1), which is composed of the forcing f and the phase sensitivity Z, and we introduce some practical constraints on f and Z.
We first consider a general class of periodic functions θ f ( ) as the weak forcing, namely, those satisfying the following constraints , namely, that f is an L p -function on π π ≡ − S [ , ]. Henceforth, we assume ⩾ p 1, due to the following physical interpretation of the constraints (2). First, for p = 2, average square of f (the power) is fixed at M 2 , which is the case considered in [16,17,21], and this case corresponds to the power-reduced forcing in P1 1 . (An example is shown in figure 1(a).) For p = 1, = f M || || p corresponds to the area-reduced forcing in P2: figure 1(b).). On the other hand, the case of = ∞ p implies the magnitude-reduced forcing in P2 ( figure 1(c) Thus, the constraints (2) continuously cover various situations in a natural way. Besides in equation (2), i.e., a charge-balance constraint [22,23], is introduced, because it is required in practical situations where total injection (injected charge) should be 0.
On the other hand, for the phase sensitivity θ Z ( ), we assume a general class of Z being twice differentiable for the case of < ⩽ ∞ p 1 , and Z being locally Lipschitz continuous for p = 1, which is required to prove existence of the optimal forcing; detailed information is given in appendix A of the supplementary material (available at stacks.iop.org/jpa/47/402002/ mmedia). We note that these (mild) assumptions are normally satisfied for the oscillators in real environments.
Phase locking occurs when the phase difference is locked, i.e., ϕ The range of frequency difference Δω, where the solution for a stable steady state exists for ϕ, defines the locking range R f [ ] for a certain fixed forcing waveform [4]. Therefore, the locking range is the difference between the maximum (at ϕ ϕ = + ) and minimum (at ϕ ϕ

Fundamental limits of synchronization
Now, we formulate the optimal synchronization problem: the global optimal forcing waveform maximizes the locking range R f [ ] under the constraints (2), which gives the maximum width of the Arnold tongue for the weak forcing, as in figure 1. If we assume f to satisfy is introduced, where λ is a Lagrange multiplier. Moreover, J f [ ] is rewritten as the following inner product of f and g where where ϕ * denotes a stable fixed point for equation (1). Note this ϕ * is set to 0 on the new coordinate.) if Δω = 0 in (1) [17]. Thus, for the case of maximizing the stability S f [ ], the same arguments for maximizing the locking range R f [ ] are possible simply by replacing θ (3). Henceforth, we restrict our argument to the case of maximizing locking-range in the proof below, for the sake of simplicity.  (2), and hereafter we denote the global optimal forcing for a given p as f p opt, . For this optimization problem, local optimal forcings can be captured by using the calculus of variations, such as the Euler-Lagrange equation for > p 1 (or the bang-bang principle [24] for = ∞ p ). But, this approach has a limitation: it is intrinsically local and heuristic, and its result lacks global information. Namely, we cannot understand 'how' and 'when' the global optimal forcing is realized. Furthermore, it is impossible using this approach to show that a certain arbitrarily tall pair of pulses realize the entrainment limit for p = 1 in P2, as explained later. Hence, by using only the calculus of variations, it is hard or impossible to answer questions P1, P2, and P3 regarding the physical limit of synchronizability.
However, if we realize that equation (3)  , then answers to the basic questions P1, P2, and P3 are systematically obtained as theorems and their proofs are obtained, as follows. holds if and only if there exist constants r and s, not both 0, such that θ θ = r f sg | ( )| | ( )| p q for almost all θ ∈ S [26]. Having this equality condition in mind, a candidate of the optimal f * for J f [ ] under the constraints (2) is given as 1 with σ θ ( ) being any function having either ±1 values for θ ∈ S, by assuming the second equality holds in the following general relationship: . Further, by assuming the first equality in the above relationship, we can narrow the above f * uniquely to where sgn is the signum function defined as sgn Here we call this f * as an ideal forcing since it realizes a possible ideal entrainment of the maximum locking range, and we now assume such an f * to exist (which is later verified from equation (5)). Under this assumption, for any given Z, fg . Then, the ideal locking range J f [ * ] is a function of Δϕ and λ, for a given Z and p. In order to maximize J f [ * ], the function θ λ + Z |¯( ) | q should be maximized by tuning the two parameters Δϕ and λ under the constraints (2). For this purpose, we define the following functions: ( , ) ( , ) is introduced, where μ is a Lagrange multiplier. Thus, the optimal entrainment problem is reduced down to the finite-dimensional optimization of Δϕ λ H ( , ), and the above argument clarifies the mechanism of how optimal forcings are realized.
Straightforward calculations show that the optimal solutions Δϕ λ ( * , * ) to Δϕ λ H ( , ) are determined from the following equations , and  H ( ) represents the bordered Hessian matrix of H; detailed information on  H ( ) is given in appendix B of the supplementary material (available at stacks.iop.org/jpa/47/402002/mmedia). Note, for every Δϕ λ ( * , * ), the associated f * indeed maximizes the associated Γ at ϕ + and minimizes Γ at ϕ − by straightforward calculation, although its detail is omitted here.
To determine Δϕ λ ( * , * ) of the ideal forcing f * in equation (4), we have numerically solved equations (5a) and (5b), and checked whether the obtained Δϕ λ ( , ) satisfies equation (5c); the results for the example in figure 2 are listed in appendix C of the supplementary material (available at stacks.iop.org/jpa/47/402002/mmedia). Since all possible Δϕ λ ( * , * ) can be obtained numerically, and from the above argument concerning the equality condition and the associated equation (4), all global or local optimal forcings are captured in L s ( ) p 2 . Among them, the global optimal (i.e., the best) forcing f p opt, is identified, simply by comparing the associated locking ranges R f [ ], as shown in figures 2(a) and (b) for = ∞ p 1, 1.01, 2, 5, for θ Z ( ) shown in figure 2(c), respectively. The above steps constitute the algorithm for realizing the global optimal forcing, i.e., the fundamental limit of injection-locking 3 . Thus, the answer to P1 has been obtained. Overview of all optimal forcings for various ∈ ∞ p [1, ] obtained for the Hodgkin-Huxley (HH) neuron phase model [17]. Black, red, and blue curves represent the global optimal, the second optimal, and Z, respectively. Panel (a), (b), respectively, show 1 : 1 and 1 : 2 phase locking optimal forcings for the HH neuron model [17], having the associated . This implies that these optimals belong to a broader class of functions, compared with the one considered in the calculus of variations. More precisely, in the calculus of variations, we usually assume the function space to consist of piecewise-smooth functions (or absolutely continuous functions, at best), which is a subset of L p . Namely, the result here is stronger than the one obtained by the conventional calculus of variations. 3 Note, for the special case of p = 2, equation (5b) gives λ = 0, equation (5a) becomes θ θ Δϕ ) 0, and hence the result for the power-reduced forcing in [16] is naturally recovered. This implies the optimal forcing obtained in [16] (as well as in [17]) is nothing other than the global optimal forcing in the more general function space of forcings, i.e., L S ( ) 2 (cf footnote 2).
3.2. 1 : 1 phase locking for p = 1 and p ¼ ∞ there is only one exceptional case: m : 1 phase locking, where an 'asymptotically' best forcing is constructed 6 . Thus, the answer to P3 has been obtained. Figure 2(b) shows the global optimal 1 : 2 phase locking forcings of the Hodgkin-Huxley (HH) neuron phase model [17] for various p. We observe that these forcings have simpler waveforms (compared with the 1 : 1 cases in figure 2(a)), uniformly spaced bipolar pulses (p = 1), nearly sinusoidal (p = 2), and almost uniformly spaced rectangles ( = ∞ p ), which is a typical feature for 1 : n global optimal forcing. The reason is as follows. In 1 : n phase locking, the above θ Z ( ) n determines the optimal forcings (through the algorithms related to P1 and P2); only the multiples of n-th harmonics in the original θ Z ( ) affect the optimal forcings. In addition, as we observed the Fourier coefficients of Z in figure 2 caption, the second (n = 2) harmonic dominates other higher harmonics, resulting in a virtually sinusoidal oscillation (analogous to the one near the Hopf bifurcation point). This situation becomes more typical if we consider a larger n, since the magnitude of higher multiples of n-th harmonics in θ Z ( ) decays sufficiently fast as the multiples of n becomes large for generic limit-cycle oscillators. In fact, this insight explains the reason for the more sinusoidal-like best forcings systematically obtained for a larger n in the general m : n case (under a different but similar setting) [21].

Conclusion and discussion
In conclusion, we have proved a mechanism governing synchronization (injection-locking) limits, and clarified how and why the best forcing realizes the synchronization limit, by unveiling a hidden aspect, i.e., Hölderʼs inequality, behind it, for a general class of externally forced limit-cycle oscillators. Namely, synchronization limit is now characterized as the equality condition of Hölderʼs inequality. To the best of the authorʼs knowledge, no previous study has addressed this mechanism or the existence of the synchronization limit. Since the phase equation (1) appears in many areas of science and engineering, the obtained results here have direct, broad impacts, as follows.
(i) Designing the best forcing (beyond power-reduced forcings, including pulse trains) and (simultaneously) designing a better phase sensitivity (i.e., oscillators with better entrainability) is one of the direct consequences, since optimization of the forcing f and that of Z are now equivalent by permuting Z and f in equation (3). Such an application is feasible for a practical system in real environments, for instance as in [27,28]. Also, it is noted that our results for the area-reduced forcings here enables us to design the optimal pulse trains for injection-locking for the first time. (ii) In addition, designing an efficient 'coupling' between oscillators for better mutual synchronization [29] is also promising, since the results of a single oscillator with one external forcing here is modified to the case of mutual synchronization straightforwardly. (iii) Furthermore, the algorithm for realizing the optimality here should provide a unified, systematic method for optimally entraining a given oscillation pattern in an ensemble of oscillators (or excitable elements) with global coupling [6], as well as with local coupling 6 The construction is as follow. Starting from m copies with the optimal forcing with prime period T 0 for 1 : 1 phase locking, add a certain small perturbation such that the m copies of the forcing become a single forcing with prime period mT 0 while still satisfying the constraints (2). The resulting locking range becomes arbitrarily close to the ideal one (which is realized only in 1 : 1 phase locking) as the perturbation becomes smaller, since the associated Γ m/1 in equation (1) becomes arbitrarily close to the Γ 1 1 of the best forcing for 1 : 1 phase locking. [8], since this oscillation pattern is regarded as a limit-cycle oscillation in higher phase space, which is described by equation (1). (iv) Though our present framework focuses on the noiseless case, noisy oscillators can be treated in the same way as here, by virtue of the recent progress in this direction [30,31]. These applications and extensions (i)-(iv) of our findings are of high importance for further theoretical as well as experimental work, which will be reported in the near future.