Phase reduction theory for hybrid nonlinear oscillators

# Phase reduction theory for hybrid nonlinear oscillators

Sho Shirasaka Graduate School of Information Science and Engineering, Tokyo Institute of Technology, O-okayama 2-12-1, Meguro, Tokyo 152-8552, Japan    Wataru Kurebayashi Faculty of Software and Information Technology, Aomori University, Kobata 2-3-1, Aomori, Aomori 030-0943, Japan    Hiroya Nakao School of Engineering, Tokyo Institute of Technology, O-okayama 2-12-1, Meguro, Tokyo 152-8552, Japan
July 3, 2019
###### Abstract

Hybrid dynamical systems characterized by discrete switching of smooth dynamics have been used to model various rhythmic phenomena. However, the phase reduction theory, a fundamental framework for analyzing the synchronization of limit-cycle oscillations in rhythmic systems, has mostly been restricted to smooth dynamical systems. Here we develop a general phase reduction theory for weakly perturbed limit cycles in hybrid dynamical systems that facilitates analysis, control, and optimization of nonlinear oscillators whose smooth models are unavailable or intractable. On the basis of the generalized theory, we analyze injection locking of hybrid limit-cycle oscillators by periodic forcing and reveal their characteristic synchronization properties, such as ultrafast and robust entrainment to the periodic forcing and logarithmic scaling at the synchronization transition. We also illustrate the theory by analyzing the synchronization dynamics of a simple physical model of biped locomotion.

###### pacs:
Valid PACS appear here

## I Introduction

Hybrid dynamical systems have been used to describe physical processes that exhibit sudden qualitative changes or abrupt jumps during otherwise continuous evolution. Some examples are the collision of particles, refraction and reflection of waves, spiking of neurons, switching of gene expression, limb-substrate impacts in legged robots and animals, human-structure interaction, switching of elements in electric circuits, and breakdown of nodes or links in networked systems Aranson2006patterns (); Wolf99geometry (); Aihara2010theory (); Holmes2006dynamics (); Macdonald09lateral (); Andronov1987theory (); helbing2003modelling+Hespanha2001hybrid+Susuki09hybrid (). Because such discontinuous events are found in many areas of science and engineering Makarenkov2012dynamics (); Bernardo08Piecewise (), it is important to develop theoretical frameworks to analyze hybrid dynamical systems.

Many hybrid dynamical systems exhibit stable rhythmic activities, for example, periodic spiking of neurons, rhythmic locomotion of robots, oscillations in power electric circuits, and business cycles in economic models coombes2012nonsmooth (); McGeer90+Zimmermann07 (); Banerjee01 (); jenkins2013self (), which are typically modeled as nonlinear limit-cycle oscillations. Synchronization of such rhythmic activities may play important functional roles in biological and engineered systems, e.g., in locomotor rhythms, vibro-impact energy harvesters, and wireless sensor networks Ijspeert98central+Strogatz05 (); Dunnmon11power+Moss10low (); tanaka2009self ().

One of the fundamental mathematical frameworks for analyzing limit-cycle oscillations is the phase reduction theory winfree2001geometry (); kuramoto2003chemical (); hoppensteadt1997weakly (), which gives approximate reduced description of the dynamics of a weakly perturbed limit-cycle oscillator using a simple one-dimensional phase equation. The phase reduction theory is well established for stable limit-cycle oscillations of smooth dynamical systems and has successfully been applied to the analysis of rhythmic spatiotemporal dynamics in chemical and biological systems winfree2001geometry (); kuramoto2003chemical (). Methods for optimizing and controlling synchronization of limit-cycle oscillators have also been developed on the basis of the phase reduction theory Dorfler14synchronization+Harada10optimal+Zlotnik13optimal+Zlotnik16phaseselective ().

In phase reduction theory for smooth dynamical systems, a weakly perturbed limit-cycle oscillator described by is considered, where is the oscillator state, is a continuously differentiable vector field representing the dynamics of the oscillator, denotes external perturbation applied to the oscillator, and is a small parameter representing the intensity of the perturbation. A system without perturbation is assumed to possess a stable limit-cycle solution of period , and a phase of the oscillator state is introduced, which increases with a constant frequency and takes the same value on the isochron winfree2001geometry (); hale69 (); guckenheimer75 (), i.e., a codimension-one manifold of the oscillator states that share the same asymptotic behavior.

When the perturbation is sufficiently small, phase reduction theory enables us to systematically approximate the original multidimensional system using a simple one-dimensional reduced phase equation of the form where is the oscillator phase and gives the phase of the oscillator state . The range of the phase is identified with a one-dimensional torus . The function , which is the gradient of the isochron and is called the phase sensitivity function in this study, quantifies the linear response of the oscillator phase to perturbations given at phase on . It is known that can be obtained as a -periodic solution to the adjoint linear problem of the system, with a normalization condition , where denotes the Jacobi matrix of and its transpose malkin1949methods (); hoppensteadt1997weakly (); ermentrout2010mathematical ().

Recently, the phase reduction theory has been extended to nonconventional cases such as stochastic Yoshimura08+Teramae09+Goldobin10 (), delay-induced Novicenko12+Kotani12 (), collective kawamurai2008collective+kori2009collective (), spatially extended kawamura2013+nakao2014 (), and strongly modulated kurebayashi2013phase () oscillations. Similar reduction methods that rely on sets of initial conditions characterized by the same long-term behavior have also been developed for heteroclinic orbits shaw2012phase (), limit tori demir2010+kawamura2015+mauroy2012 () and stable fixed points Ichinose98+Rabinovich99+Mauroy13 (). However, application of the phase reduction theory to oscillatory hybrid dynamical systems has so far been limited to low-dimensional systems, or to a specific class of systems whose phase sensitivity function is obtained from adiabatic approximation coombes2012nonsmooth (); izhikevich2000phase (); park2013infinitesimal (). To the best of our knowledge, no systematic phase reduction theory for oscillators of high-dimensional () systems with discontinuity in has been developed. This is mainly because the nonsmoothness of the vector fields at the jumps prevents straightforward utilization of the adjoint equation.

In this study, we develop a systematic phase reduction theory for a general class of autonomous limit-cycle oscillators in hybrid dynamical systems. This paper is organized as follows: in Sec. II, limit-cycle oscillations in hybrid dynamical systems are introduced. In Sec. III, the phase reduction theory for hybrid limit cycles is developed. In Sec. IV, the theory is illustrated by analyzing synchronization dynamics of two examples of hybrid limit-cycle oscillators, that is, an analytically tractable Stuart-Landau-type oscillator and a physical model of biped locomotion. Section V summarizes the results, and Appendices A-I provide mathematical details of the main results presented in Sec. II-IV.

## Ii Hybrid Limit Cycles

The state of a hybrid dynamical system that we consider in this study is represented by a pair of the discrete state for some ( is allowed) and the continuous state . We denote the set of pairs as , which is a collection of all possible transitions from discrete state to . As in Ref. khan2011sensitivity (), we describe a hybrid dynamical system by the following hybrid automaton:

 ˙X(t)=F(I(t),X(t)),\ \ if\ \  I(t)=i and X(t)∉Πij for any j, (1)
 X(t+0)=Φ((i,j),X(t)), I(t+0)=j,\ % \ if\ \ I(t)=i and X(t)∈Πij for some j, (2)
 Πij={{X|L((i,j),X)=0},\ \ if\ \ (i,j)∈G.empty set, \ \ otherwise. (3)

Here, Eq. (1) describes the smooth dynamics of the continuous state when the discrete state is , Eq. (2) the jump of when the discrete state switches from to , and Eq. (3) represents a plane in the space of continuous state on which the switching from to takes place. In Eq. (1), is the vector field of the system. In Eq. (2), “” indicates the moment just after the switching of the discrete state at , the transition function gives the new continuous state after the switching of the discrete state from to , and is an dimensional zero-level surface of the function on which the switching takes place. It is assumed that the functions , , and are continuously differentiable with respect to and do not depend explicitly on time.

Suppose there exists a periodic solution of period of Eqs. (1)-(3). As in Ref. Akhmet2005on+akhmet2010principles (), we make several assumptions on the system (see Appendix A for details) so the continuous part of the solution is piecewise continuously differentiable with respect to the initial continuous state on and linear stability analysis of the solution can be performed. Let be an initial condition of Eqs. (1)-(3) and , be the moments of switching of the discrete state, where . For convenience of notation, we also define .

To simplify the expression of the periodic orbit , we hereafter use the following notation. By distinguishing the discrete states visited more than once in one period, and by renumbering the state indices, we introduce a set of discrete states where , such that the discrete state is switched in numerical order as at (see Fig. 1). Here, is finite because the period is finite and the assumption (C2) in Appendix A assures that the system stays in each discrete state for some nonzero duration. We also introduce the following simplified notations for the discrete state transitions on the periodic orbit :

 Lk(X0(t))=L((k,k+1),X0(t)), (4) Φk(X0(t))=Φ((k,k+1),X0(t)). (5)

We call the periodic solution a hybrid limit cycle if it is linearly stable (see Appendix B for the linear stability analysis of the periodic solution).

## Iii Phase reduction

The aim of phase reduction is to describe the dynamics of the system state around the hybrid limit cycle by using a scalar phase . We first introduce the phase function on , which gives the phase value of the state on as

 Θ(s0(t+nT))=Θ(I0(t+nT),X0(t+nT))=t (mod T), (6)

where is an arbitrary positive integer. Namely, we identify the time (mod ) as the oscillator phase , which increases with a constant frequency on , i.e.,

 ˙θ(t)=˙Θ(I0(t),X0(t))≡1. (7)

In the following, we will denote a system state with phase on also as as a function of .

Phase can also be introduced in a neighborhood containing within its basin of attraction by introducing an equivalence relation to initial conditions in whose asymptotic behaviors are the same. Namely, we introduce the isochron of by assigning the same phase value to the set of states in that eventually converge to the same state on . Suppose that and are taken from , where is on at phase , i.e., . If and are asymptotically equivalent, we define the phase of as

 Θ(s1)=Θ(s2)=θ. (8)

Note that the convergence concept of the solutions in hybrid dynamical systems demands somewhat careful attention (see Appendix C for details). Some properties of the isochron and the phase function on are discussed in Appendix D.

The above definition of the phase guarantees that the following relation holds for almost all (excluding the Lebesgue measure zero set of the moments of switching) for an unperturbed oscillator:

 ˙θ(t) =˙Θ(I(t),X(t))=∇Θ(I(t),X(t))⋅F(I(t),X(t))=1, (9)

where represents the gradient of with respect to . Namely, the phase rotates on a circle at a constant frequency .

When a sufficiently small perturbation with is introduced to the oscillator as

 ˙X(t)= F(I(t),X(t))+ϵp(I(t),X(t),t), (10)

we can obtain the following approximate phase equation closed in at the lowest order:

 ˙θ(t) =˙Θ(I(t),X(t)) (11) =1+ϵ∇Θ(I0(θ),X0(θ))⋅p(I0(θ),X0(θ),t)+O(ϵ2) (12) =1+ϵZ(θ)⋅p(I0(θ),X0(θ),t)+O(ϵ2), (13)

where we defined the phase sensitivity function

 Z(θ)=∇Θ(I0(θ),X0(θ)) (14)

characterizing the linear response property of the oscillator phase to perturbations. We consider that the phase evolves as a solution of a suitably regularized, multivalued system of Eq. (13), such as the Filippov system Filippov2013 (); Cortes2012discontinuous (). (We do not consider impulsive perturbation at the moment of switching in this study, which requires special treatment.) The ideas underlying the phase approximation Eq. (13) and some notes on the notion of the solution of it are given in Appendix E.

Thus, once we obtain the phase sensitivity function , the dynamics of a weakly perturbed hybrid limit cycle described by Eq. (10) can be reduced to a single phase equation (13). Using the reduced phase equation, we can analyze various synchronization dynamics of hybrid limit cycles in detail. As we derive in Appendix F, is given by a periodic solution to the following adjoint system:

 ˙Z(t)=−A†(k,t)Z(t)\ \ for% \ \ t (mod T)∈(τk−1(s∗),τk(s∗)), (15) Z(t)=(Ck)†Z(t+0)\ \ at\ % \ t (mod T)=τk(s∗), (16)

which is normalized to satisfy Eq. (9) on , that is,

 Z(t)⋅F(I0(t),X0(t))=1. (17)

Here, is the Jacobi matrix of estimated on , and is a “saltation matrix” Bernardo08Piecewise () given by

 Ck =DΦk(X0(τk(s∗))) (18) −[DΦk(X0(τk(s∗)))˙X0(τk(s∗))−˙X0(τk(s∗)+0)]⊗(∇Lk(X0(τk(s∗)))∇Lk(X0(τk(s∗)))⋅˙X0(τk(s∗))), (19)

where is the Jacobi matrix of and represents a tensor product of two vectors. represents expansion or contraction of small deviations from by the mapping at the switching , where the second term on the right-hand side takes into account the shift in the switching time caused by the perturbation. In general, the above adjoint system can be integrated only backward in time because can be singular.

In numerical calculations, we integrate these adjoint equations backward in time with occasional renormalization of so Eq. (17) is satisfied. Then, reflecting the linear stability of , all modes except the neutrally stable periodic solution decay and is eventually obtained. This is a standard procedure for calculating of ordinary limit cycles and is called the adjoint method after Ermentrout ermentrout2010mathematical ().

## Iv Examples

As an application of the phase reduction theory for hybrid limit cycles that we developed, we analyze injection locking of hybrid limit-cycle oscillators, i.e., synchronization of the oscillator to a periodic external signal kuramoto2003chemical (). We apply a weak periodic signal to the hybrid limit cycle described by Eqs. (1) and (2). Using the phase reduction theory, the state of the perturbed oscillator, described by Eq. (10), can be approximately represented by its phase , which obeys the reduced phase equation (13).

To analyze the synchronization dynamics, we consider the phase difference between the oscillator and the periodic signal,

 ψ=θ−TTextt, (21)

where is the phase of the hybrid limit cycle, is the natural period of the hybrid limit cycle, and is the period of the external periodic signal. The frequency mismatch between the oscillator and the signal is given by . As shown in Appendix G, using the standard averaging approximation for weakly perturbed oscillators kuramoto2003chemical (); hoppensteadt1997weakly (); ermentrout2010mathematical (), the dynamics of can be derived from the reduced phase equation (13) as

 ˙ψ=ϵ[Δ+Γ(ψ)] (22)

where the -periodic phase coupling function is given by

 Γ(ψ)=1Text∫Text0Z(TTextt+ψ)⋅p(t)dt. (23)

As mentioned previously for the phase equation, we consider that the phase difference evolves as a solution of Eq. (22) in a regularized sense, if necessary. See Perestyuk2013averaging+Perestyuk2011differential () for the averaging approximation in nonautonomous systems with jumps and multivalued righthand sides.

Synchronization dynamics of the oscillator can easily be understood from the phase coupling function . If Eq. (22) has a stable fixed point, then the phase difference converges to this point and the oscillator is phase locked to the periodic signal. If there exist multiple stable fixed points, then the oscillator can be phase locked to the periodic signal at multiple phase differences depending on the initial condition. If Eq. (22) does not have a fixed point, continues to increase or decrease and phase locking does not occur.

### iv.1 Glued Stuart-Landau oscillator

As the first example, we introduce an analytically tractable model of a hybrid limit-cycle oscillator, which is constructed by gluing two Stuart-Landau oscillators (normal forms of the supercritical Hopf bifurcation guckenheimer83nonlinear ()) of different amplitudes. The glued Stuart-Landau oscillator has two discrete states and a two-dimensional continuous state variable , where denotes the transpose of a matrix. The dynamics is described by

 F(1,X)=(x−ay−(x2+y2)(x−by)ax+y−(x2+y2)(bx+y)), F(2,X)=(x−ay−α2(x2+y2)(x−by)ax+y−α2(x2+y2)(bx+y)), Φ1(X)=(xαy), Φ2(X)=(αxy), Π1,2={X | (L1(X)=0) ∩ (x≤0)},Π2,1={X | (L2(X)=0) ∩ (x≥0)}, L1(X)=y, L2(X)=−y, (27)

and the parameters are set as and . With these parameter values, this system has a stable limit cycle of period . We take the origin of the phase at and , i.e., .

The periodic orbit is depicted on  (Fig. 2(a)), which satisfies

 (I0(θ),X0(θ))=(1,(−sin(2πθ),cos(2πθ))) (28)

for , and

 (I0(θ),X0(θ))=(2,(−0.5sin(2πθ),0.5cos(2πθ))) (29)

for , where and are domains of the phase.

The phase sensitivity function can be obtained by solving the adjoint linear problem analytically and is given by

 Z(θ)=−12π(cos2πθ−sin2πθ,sin2πθ+cos2πθ)† (30)

for , and

 Z(θ)=−1π(cos2πθ−sin2πθ,sin2πθ+cos2πθ)† (31)

for .

For the waveform of the periodic injection signal , we consider rectangular waves:

 p1(t)=−c (if t mod Text∈D),0 (otherwise), (32)

and

 p2(t)=0 for all t, (33)

where the constant is set so that the squared mean of the injection signal becomes unity, i.e., unless otherwise specified, and is the time domain where the forcing takes place.

In Figs. 2(b) and 2(c), the phase sensitivity function obtained by analytically solving the proposed adjoint systems is compared with the result of the direct numerical simulation (see Appendix H for details). The results are in good agreement and show the validity of the proposed adjoint method. The discontinuities in are characteristic of a hybrid limit-cycle oscillator.

Figure 3(a) displays the averaged dynamics of the phase difference for several initial values, where the result of phase reduction is compared with direct numerical simulations. It can be seen that the asymptotic phase differences and their dependence on initial conditions are well predicted from the phase coupling function . Figure 3(b) shows the boundaries of the region where the injection locking takes place, called the Arnold tongue hoppensteadt1997weakly (). Results of the numerical simulation also agree well with the analytical prediction by the phase reduction theory. Thus, the injection locking of hybrid limit-cycle oscillators by weak periodic input can be theoretically predicted by using obtained by the adjoint method.

Here, we emphasize one peculiar property of the hybrid oscillator. For smooth oscillators, the period of the phase slipping near the critical point generally obeys the inverse square-root scaling law, , where is the critical value of  kuramoto2003chemical (); izhikevich2010dynamical (), because the phase coupling function generally has quadratic maximum and minimum. In contrast, as shown in Fig. 3(c), the hybrid oscillator with discontinuous can possess nonsmooth with sharp nonquadratic maximum or minimum when subjected to impulsive signals. For such , we can show that the period of phase slipping obeys

 ϵTslip∼−ln|Δ−Δc| (34)

at the leading order in the vicinity of the critical value , where and is the extremum at the corner of , and the semiderivatives and are nonzero (see Appendix I for the derivation). Since nonsmooth corners in cannot exist in smooth systems, singular scaling law of this kind is characteristic of hybrid oscillators.

In Fig. 3(d), the transient dynamics of for two different types of rectangular wave input, one with a low-duty ratio (impulsive) and the other with a mediate one (mild, nonimpulsive) , is compared. Here the magnitude of the input signal is normalized so that the uniform norm of becomes unity, i.e., , for each case. For the impulsive case, approaches the stable phase difference with a nonzero angle, while in the nonimpulsive case, the approach is tangential. This implies that the decay of the deviation from is faster than exponential in the impulsive case and the time required to establish entrainment is drastically shorter. Moreover, variations in the input period only slightly changes for the impulsive input, while shows significant change for the nonimpulsive input. This ultrafast entrainment and robustness of the stable phase difference can be attributed to the existence of the region where is extremely steep.

Note that the discontinuity in is necessary for the existence of such a region , because is given by the convolution (23) of and ; when the input is an ideal impulse, the slope of in can be infinite. Therefore, these interesting synchronization properties are distinctive feature of the hybrid oscillators driven by impulsive periodic input. Such a type of very fast (or finite-time) synchronization has been studied for simple neuron models with discontinuity whose can be obtained analytically, as well as in some fast-slow models in the fast-relaxation limit izhikevich2000phase (); veryfast (); veryfast2 (). Our argument based on the phase reduction theory for hybrid limit cycles is general and can be applied to high-dimensional systems where the nonsmoothness is not the result of adiabatic approximation.

### iv.2 Passive bipedal walker

Next, we analyze a physical example of hybrid limit-cycle oscillator, namely a two-link model of a passive walker walking down a slope garcia1998simplest (), proposed as a simple model of biped locomotion. Figure 4(a) shows a schematic diagram of the model, where is the gravitational acceleration; is the length of the legs; and are the masses of the hip and the foot, respectively; and specify the angles of the swing and support legs; is the angle of the slope; and is a periodic torque applied to the ankle of the support leg. It is assumed that , i.e., the hip mass is much larger than the foot mass, the tip of the support leg does not slip along the ground, and the collision of the foot with the ground is perfectly inelastic (no slip and no bounce). This model exhibits a stable limit-cycle oscillation for appropriate parameter values that corresponds to periodic movements of the legs. This is a four dimensional hybrid dynamical system with impacts, hence it cannot be dealt with by the conventional methods coombes2012nonsmooth (); izhikevich2000phase (); park2013infinitesimal () mentioned above.

The model has a one discrete state and a continuous state variable . The dynamics is described by

 F(1,X)=⎛⎜ ⎜ ⎜ ⎜ ⎜⎝˙ϕ1sin(ϕ1−γ)˙ϕ2sin(ϕ1−γ)+˙ϕ21sinϕ2−cos(ϕ1−γ)sinϕ2⎞⎟ ⎟ ⎟ ⎟ ⎟⎠, Φ1(X)=⎛⎜ ⎜ ⎜ ⎜⎝−ϕ1˙ϕ1cos2ϕ1−2ϕ1˙ϕ1cos2ϕ1(1−cos2ϕ1)⎞⎟ ⎟ ⎟ ⎟⎠, Π1,1={X | (L((1,1),X)=0)∩(ϕ2<−δ)}, L1(X)=2ϕ1−ϕ2, (37)

where we have rescaled time by , and is a small positive constant (we set ), which is introduced to avoid foot scuffing (contact of the swing leg with the ground in the middle of the swing). The parameter is set as . Note that Eq. (IV.2) is an equation of motion representing continuous dynamics of the walker during the single-leg support phase, in which the walker stands on the support leg and moves the swing leg, where and are the angular coordinates of the swing and support legs. See Ref. goswami1996compass () for a detailed derivation of the above type of equations in nearly the same setting. Note also that the main effect of the inclined ground to the walking dynamics, i.e., the ground reaction force, is already included in the dynamical model. This effect is not considered a perturbation and therefore need not be weak, as long as the model exhibits stable rhythmic walking. If there exist additional small effects from the flat inclined ground, such as slight up and down, they can be incorporated into the reduced phase model perturbatively. Finally, although the motion of Eq. (IV.2) during the single-leg support phase appears to be conservative, the collision of the leg with the ground, described Eq. (IV.2) and Eq. (37), is perfectly inelastic (plastic), so the impact of the leg with nonzero velocity relative to the ground causes energy dissipation. This energy loss is compensated by the gravitational potential energy, which is supplied to the system at each moment of the collision of the leg with the ground. Thus, a stable limit cycle can arise in this hybrid dynamical system.

Figure 4(b) shows the stable periodic orbit of the model. Using the shooting method developed in khan2011sensitivity (), a point on , which we define as the origin of the phase, and the period of the orbit can be obtained as

 s∗=(1,(0.009000,−0.05869,−0.0009629,−0.3432)†),T=3.882. (38)

Figure 4(c) shows the phase sensitivity function with discontinuity in the middle, which is obtained numerically by the proposed adjoint method. The result agrees well with the one obtained by the direct method, thus the proposed adjoint method also works nicely for this model.

Using the reduced phase equation, we study the injection locking of the passive walker to the periodic ankle torque actuation. That is, we apply a weak periodic torque to the walker, where the frequency of the torque is close to that of the natural frequency of the walker, and analyze whether the walker synchronizes with the applied weak torque. We introduce periodic actuation of the ankle torque as the injection signal , where

 τ(t)=ce−0.5[t (mod Text)]/Textsin(4πt/Text). (39)

The magnitude of the waveform is determined to satisfy the normalization condition . As in the case of the glued Stuart-Landau oscillator, we can obtain the phase coupling function from the phase sensitivity function and the injected periodic signal , and predict the dynamics of the phase difference between the oscillator and the signal.

Figure 5(a) shows the dynamics of . The reduced phase equation predicts that there are two stable fixed points of , and direct numerical simulations of the original model from different initial conditions confirm that the phase difference of the passive walker is actually attracted to either of the stable fixed points. Figure 5(b) plots the Arnold tongue showing the region where the phase locking of the passive walker to the injected signal takes place. The results obtained by the phase reduction theory agree well with the results obtained by direct numerical simulations of the original model.

Thus, the proposed phase reduction theory is also useful in analyzing realistic physical systems, even when the hybrid limit-cycle oscillator has a high-dimensional continuous state.

## V Summary

We formulated a phase reduction theory for a general class of hybrid limit-cycle oscillators and derived the adjoint equation for the phase sensitivity function. The proposed theory provides precise phase sensitivity functions and the derived phase equation accurately predicts the injection locking properties of hybrid oscillators. We illustrated synchronization properties characteristic to hybrid oscillators, such as ultrafast entrainment to periodic signal and negative logarithmic scaling at the synchronization transition, and explained them by using the reduced phase equation with discontinuous phase sensitivity functions.

The phase reduction theory developed in this study would serve as a powerful tool for investigating synchronization phenomena in complex nonsmooth systems and for finding various applications in controlling distributed interacting nonlinear oscillators Ijspeert98central+Strogatz05 (); Dunnmon11power+Moss10low (); tanaka2009self ().

###### Acknowledgements.
We acknowledge the valuable comments on this work of Jaap Eldering, Shigefumi Hata and Hiroshi Kokubu. This study was supported by Japan Society for the Promotion of Science (Japan) KAKENHI (Grants No. 22684020, No. 25540108, No. 26103510, No. 26120513 and No. 15J12045), the Core Research for Evolutional Science and Technology (Japan) Kokubu project of Japan Science and Technology Agency (Japan), and the FIRST Aihara project of Japan Society for the Promotion of Science (Japan).

## Appendix A Assumptions for the periodic solution

In this section, we introduce the assumptions that are necessary for the periodic solution to be piecewise continuously differentiable with respect to the initial condition Akhmet2005on+akhmet2010principles (). Hybrid dynamical systems can exhibit pathological behaviors, which do not occur in smooth dynamical systems, such as the grazing, livelock, sliding, and Zeno phenomena due to the effect of discrete switching Bernardo08Piecewise (); Schaft00 (). The grazing phenomenon Bernardo08Piecewise () occurs when the orbit becomes tangent to the switching surface. This condition can be written as

 ∇L((i,j),X0(t))|L=0⋅˙X0(t)=0, (40)

where is the gradient of with respect to and denotes inner product of vectors. The livelock, sliding, and Zeno phenomena Schaft00 () can arise when the points are allowed to be accumulation points of the switching surfaces. These conditions can lead to infinite sensitivity to the initial conditions jeffrey2011nondeterminism (). In this study, we do not consider such pathological situations, namely, we assume that (C1) the orbit is always transversal to the switching plane and that (C2) each continuous state right after the discrete state transition has a neighborhood that is disjoint from the switching surfaces.

## Appendix B Linear stability of the hybrid limit cycle

In this section, we formalize the linear stability of the periodic solution. Let () be the th initial-condition sensitivity vector khan2011sensitivity () with respect to an initial state on at , defined as

 ξα(t)=limϵ→0(X(t;(I0(0),X0(0)+ϵeα))−X0(t)ϵ). (41)

Here, the second argument of represents an initial state that is slightly perturbed in the th direction ( is the th unit vector) from on . We introduce a sensitivity matrix as the collection of sensitivity vectors in all directions. Note that

 Ξ(0)=(e1,...,eN)=I (42)

where is the identity matrix, because .

In Refs. Akhmet2005on+akhmet2010principles (); khan2011sensitivity (), the linear variational system for has been derived as

 ˙Ξ(t)=A(k,t)Ξ(t)\ \ for\ \ t (mod T)∈(τk−1(s∗),τk(s∗)), (43)
 Ξ(t+0)=CkΞ(t)\ \ at\ \ t (mod T)=τk(s∗), (44)

where is the Jacobi matrix of estimated on , and is the so-called saltation matrix Bernardo08Piecewise () (Eq. (LABEL:eq._jumpmatrix) in the main article).

We define a monodromy matrix from the sensitivity matrix as . If this has one simple eigenvalue equal to 1 and all other eigenvalues lie strictly inside the unit circle on the complex plane, the periodic solution is linearly stable khan2011sensitivity (). We call such a stable isolated periodic solution a hybrid limit cycle. It can be shown that the eigenvalues and their algebraic multiplicities of monodromy matrices do not depend on the choice of the initial state. We can also consider initial-condition sensitivity vectors with respect to the state on the limit cycle , instead of , as

 ξα(t;θ)=limϵ→0(X(t;(I0(θ),X0(θ)+ϵeα))−X0(t+θ)ϵ). (45)

We denote by a matrix the collection of the sensitivity vectors in all directions, and introduce a monodromy matrix of the linear variational system, which satisfies

 ˙Ξ(t;θ)=A(k,t+θ)Ξ(t;θ)\ \ for\ \ t+θ (mod T)∈(τk−1(s∗),τk(s∗)), (46)
 Ξ(t+0;θ)=CkΞ(t;θ)% \ \ at\ \ t+θ (mod T)=τk(s∗), (47)
 Ξ(0;θ)=I. (48)

There exist a unique state transition matrix of the linear variational system (43)-(44) satisfying

 Ξ(t;θ)=Hk(t+θ,s+θ) Ξ(s;θ)\ \ for\ \ t+θ,s+θ∈[τk−1(s∗)+nT,τk(s∗)+nT), (49)

which is a solution to

 ˙Hk(t,s)=A(k,t) Hk(t,s)\ \ % for\ \ t,s∈[τk−1(s∗)+nT,τk(s∗)+nT], (50)

with the initial condition .

Suppose that is in the interval , where . Using , the monodromy matrix can be expressed as , where

 M1=H1(T,τm0(s∗)) ( m0∏∏+l=m∗Hl+1(τl+1(s∗),τl(s∗))Cl)⋅Hm∗(τm∗(s∗),θ), (51)

and

 M2 =Hm∗(θ,τm∗−1(s∗))⋅m∗−1∏∏+k=1CkHk(τk(s∗),τk−1(s∗)). (52)

Here, we denote by the ordered product of matrices. With these matrices and , one can also show that , where , which follows from , is used. Therefore, and has the same set of eigenvalues with the same algebraic multiplicities, because the Jordan blocks with nonzero eigenvalues of the products of the matrices and are identical (horn2012, , Th. 3.2.11.1.).

## Appendix C Asymptotic equivalence of initial conditions in hybrid dynamical systems

In this section, we introduce an asymptotic equivalence of the initial conditions that we use for defining the isochrons, which differs from the one used in smooth systems. Suppose and are the continuous parts of the solutions to a system with initial conditions and , respectively. In smooth systems, the asymptotic equivalence relation of the initial conditions and is defined as the convergence of the error in the Euclidean topology. That is, if where is the Euclidean norm, then and are asymptotically equivalent. In hybrid dynamical systems, the moments of switching of the two solutions for these initial conditions generally do not coincide in a finite time. Hence, the Euclidean norm of the error of the continuous part causes some kind of “peaking behavior” biemond2013tracking (); can be larger than a constant (of order ) for some for an arbitrarily large due to the continuous state jumps. This violates the definition of convergence in Euclidean topology. Therefore, we need to consider convergence in some other suitable topology to define the asymptotic equivalence notion in hybrid dynamical systems.

Various topologies suitable for hybrid dynamical systems have been proposed in the literature, such as the Skorohod topology broucke2002continuous () originally designed as a tool to analyze stochastic processes, the topology of graphical convergence that is based on set-valued analysis goebel2006solutions (), and the quotient topology generated on the hybrifold simic2005towards (), which is, roughly speaking, the manifold constructed by identifying the switching surface with its image of the transition function .

In this study, we adopt an asymptotic equivalence that is based on convergence in B-topology Akhmet2005on+akhmet2010principles (), which is defined as follows: if for any , there exist such that, in the time domain , every moment of switching of the solution lies in some neighborhood of the moment of switching of the solution , and for all , which are outside the neighborhoods of the moments of switching of , holds, then we call the initial conditions and asymptotically equivalent. The benefit of this definition is intuitively clear; the error is evaluated outside of the neighborhoods of the points of discontinuity, and thus we can avoid the effect of the peaking behavior, and the error that occurs at the moment of switching is guaranteed to disappear.

## Appendix D Some properties of the isochron and the phase function

Using the following equivalence relation (Eq.  in the main article)

 Θ(s1)=Θ(s2)=θ, (53)

we can introduce the “conditional” isochron , i.e., the set of continuous states sharing the same phase value for each discrete state . We can then define the isochron of with phase as the union . We note that the notion of the isochron in hybrid dynamical systems has been proposed in burden2014hybrid (), but the phase dynamics of weakly perturbed hybrid oscillators has not been discussed so far.

We can show that, in a domain , where is a neighborhood of such that the solution starting from the point in is piecewise continuously diffentiable with respect to the initial condition, the phase function is totally differentiable with respect to the continuous state and that it is one-sided differentiable on the switching surfaces as follows.

We consider two slightly different initial conditions, and , where , and . Note that, when , should be taken from the subset of whose elements point in the opposite direction as in the tangent space Lee2002 () of the switching boundary. From the differentiability with respect to initial conditions, the following relation holds for all outside the neighborhoods of the moments of switching,

 X2(t)=X1(t)+ϵΞs1(t)h+O(ϵ2), (54)

where and are the continuous part of the solutions with the initial conditions and , respectively, and denotes the initial condition sensitivity matrix with the initial state . Since and are taken from , which is a subset of the basin of attraction of the hybrid limit cycle , the solutions relax to , which we denote as in this section. One can easily see that the following relation holds:

 X0,2(t)=X0,1(t)+∫θ(s2)−θ(s1)0F(k,X0,1(t+t′))dt′. (55)

Since is not an equilibrium, we can assume that the first entry of is nonzero without loss of generality and denote its absolute value as . Considering the continuity of the continuous part of the solution with respect to the initial condition and the continuity of the vector field with respect to the continuous variable , there exists such that for any , the inequality