Phase and frequency linear response theory for hyperbolic chaotic oscillators

TL;DR

This paper develops a linear response theory for hyperbolic chaotic oscillators, extending phase response concepts from limit cycle oscillators to chaotic systems, with practical methods for measuring phase shifts.

Contribution

It introduces a shadowing-based phase and frequency response framework for hyperbolic chaotic flows, enabling analysis of phase shifts under perturbations.

Findings

Defines phase sensitivity function via an adjoint equation

Provides a method to construct phase and shadowing trajectories explicitly

Identifies the limits of linear response regime in chaotic oscillators

Abstract

We formulate a linear phase and frequency response theory for hyperbolic flows, which generalizes phase response theory for autonomous limit cycle oscillators to hyperbolic chaotic dynamics. The theory is based on a shadowing conjecture, stating the existence of a perturbed trajectory shadowing every unperturbed trajectory on the system attractor for any small enough perturbation of arbitrary duration and a corresponding unique time isomorphism, which we identify as phase, such that phase shifts between the unperturbed trajectory and its perturbed shadow are well defined. The phase sensitivity function is the solution of an adjoint linear equation and can be used to estimate the average change of phase velocity to small time dependent or independent perturbations. These changes of frequency are experimentally accessible giving a convenient way to define and measure phase response curves for chaotic oscillators. The shadowing trajectory and the phase can be constructed explicitly in the tangent space of an unperturbed trajectory using co-variant Lyapunov vectors. It can also be used to identify the limits of the regime of linear response.

Equations

$$\dot{\phi }=1+\epsilon \vec{Z}(\phi )\cdot \vec{p}(\phi ,t).$$

$$\frac{d{\vec{x}}_0}{d\phi }=\vec{f}({\vec{x}}_0).$$

$$\dot{\vec{x}}=\vec{f}(\vec{x})+\epsilon \vec{p}(\vec{x},t)$$

$$\frac{\Delta \phi }{\epsilon }\to Z_{\Delta }({\phi }_0)={\int }_{{\phi }_0-\tau }^{{\phi }_0+\tau }\vec{Z}({\vec{x}}_0(\phi ))\cdot {\partial }_{\mu }\vec{f}\cdot \Delta \mu (\phi )d\phi .$$

$$\frac{d}{d\phi }\vec{Z}=-J_f^{\top }(\phi )\vec{Z}$$

$$\frac{dt}{d\phi }=\frac{1}{1+\epsilon \vec{Z}\cdot \vec{p}}=1-\epsilon \vec{Z}(\phi )\cdot \vec{p}(\phi ,t)+O({\epsilon }^2).$$

$$\frac{T_{\epsilon }}{T_0}=1-\epsilon {⟨\vec{Z}(\phi )\cdot \vec{p}(\phi )⟩}_{\phi }$$

$${⟨\vec{Z}(\phi )\cdot \vec{p}(\phi )⟩}_{\phi }=\underset{\epsilon \to 0}{\lim }\frac{1}{\epsilon }\frac{T_0-T_{\epsilon }}{T_0}=\underset{\epsilon \to 0}{\lim }\frac{1}{\epsilon }\left(\frac{{\nu }_{\epsilon }}{{\nu }_0}-1\right).$$

$$z_k={⟨\vec{Z}(\phi )\cdot {\vec{p}}_k(\phi )⟩}_{\phi }.$$

$$Z_{\Delta }(\vartheta )=\underset{\epsilon \to 0}{\lim }\frac{1}{\epsilon }(T_0-T_{\epsilon })$$

$${\vartheta }_{\sigma }(\vec{x})={\vartheta }_0(\vec{x})+\sum _k{\sigma }_k{q}_k(\vec{x})\text{mod}2\pi .$$

$$\vec{\nabla }{\vartheta }_{\sigma }\cdot {\vec{v}}^{(l)}=\vec{\nabla }{\vartheta }_0\cdot {\vec{v}}^{(l)}+\sum _k{\sigma }_k\vec{\nabla }{q}_k\cdot {\vec{v}}^{(l)}={\omega }_0{\delta }_{l0}+{\eta }_l.$$

$$\frac{d}{d\phi }\vec{h}={\text{J}}_f(\phi )\vec{h}.$$

$$\frac{d}{d\phi }{h}_k={\lambda }^{(k)}(\phi )h_k.$$

$$\dot{\phi }=1+\epsilon \vec{Z}({\vec{x}}_0(\phi ))\cdot \vec{p}({\vec{x}}_0(\phi ),t)$$

$$\frac{d}{d\phi }{\vec{v}}^{(k)}$$

$$\frac{d}{d\phi }{\vec{u}}^{(k)}$$

$$\frac{d}{dt}{\vec{x}}_{\epsilon }=\vec{f}\left({\vec{x}}_{\epsilon }\right)+\epsilon \vec{p}\left({\vec{x}}_{\epsilon },t\right).$$

$$\frac{d}{d\phi }\vec{h}$$

$$\frac{d{h}_k}{d\phi }={\lambda }^{(k)}{h}_k+\epsilon {\vec{u}}^{(k)}\cdot \vec{p},\text{for }k\ne 0$$

$$\frac{d{h}_0}{dt}=\frac{dt}{d\phi }-1+\epsilon {\vec{u}}^{(0)}\cdot \vec{p},\text{for }k=0.$$

$$\frac{dt}{d\phi }=1-\epsilon {\vec{u}}^{(0)}\cdot \vec{p}.$$

$$\dot{E}$$

$$\dot{U}$$

$$\dot{W}$$

$$g(E)=2.5e^{-(E+1{)}^2}+0.01e^{\frac{1}{2}(E+6)}.$$

$${\vartheta }_{\sigma }={\vartheta }_0+\sum _{k=0}^6\sum _{l=0}^3\sum _{m=0}^3{\sigma }_{klm}^{\pm }{q}_{klm}^{\pm }$$

$$q_{klm}^{+}=\cos (k{\vartheta }_0)R^l{z}^m,q_{klm}^{-}=\sin (k{\vartheta }_0)R^l{z}^m,$$

$${\dot{\vartheta }}_{\sigma }={\dot{\vartheta }}_0+\sum _k{\sigma }_k{\dot{q}}_k={\omega }_0+{\eta }_0.$$

$$\dot{x}$$

$$\dot{y}$$

$$\dot{q}$$

$${\dot{x}}_1$$

Full text

††preprint: AIP/123-QED

**Phase response curves are a powerful tool to predict and analyze synchronization of weakly forced or coupled oscillators. The state of chaotic oscillators, however, is not characterized by a unique geometric phase. Even if a geometric phase is imposed, the phase difference between two identical chaotic oscillators is not asymptotically constant or even bounded, whereas phase response is commonly measured as the asymptotic phase shift caused by a single pulsed perturbation. In this report, we reinterpret phase as a time isomorphism rather than a geometric angle. This allows us to generalize linear phase response theory to chaotic oscillators, as well as to predict and measure the phase response via experimentally accessible frequency shifts. **

I Introduction

Synchronization, the adaptation of frequencies of self-sustained oscillators to a driving force, plays a vital role in many systems, ranging from biological and chemical systems to artificial devices Pikovsky et al. (2003); Winfree (2001); Kuramoto (2003); Glass (2001), and its understanding is essential for prediction and control of collective behavior. Synchronization can manifest in many forms, weakly as a resonance in periodically forced stochastic oscillators Tönjes and Kori (2011), or more strongly as a locking of oscillation frequencies, phases, complete or generalized synchronization Rulkov et al. (1995). Phase synchronization in weakly coupled or weakly forced, autonomous limit cycle oscillators can be understood by linear phase response theory, which describes the evolution of a phase $\varphi=\varphi(t)$, defined on a circle with the perimeter of its natural period $T_{0}=2\pi/\omega_{0}$, in linear order of a perturbation $\varepsilon\vec{p}(\varphi,t)$ as

[TABLE]

Note that in this convention, phase $\varphi$ has the dimension of time. While Eq. (1) describes the change of phase velocity in linear order of $\varepsilon$, the equation is nonlinear in $\varphi$ and even small perturbations can aggregate to nonlinear synchronization effects. Equations like (1) are sometimes referred to as Winfree type phase equations in recognition of his unifying works in mathematical biology Winfree (1967, 2001). The function $\vec{Z}(\varphi)$ is called phase sensitivity function and its components are proportional to phase response curves (PRCs). The PRCs essentially determine a system’s synchronization behavior and are used in mathematical modeling of weakly coupled oscillators across scientific disciplines from biology, in particular neuroscience Hoppensteadt and Izhikevich (2012) and chronobiology Winfree (2001), chemistry, ecology to electrical engineering and many others Pikovsky et al. (2003). Based on the PRCs it is possible to design perturbation protocols that can stabilize or destabilize various collective modes in ensembles of oscillators including complete synchronization, clustering and the asynchronous state Kiss et al. (2007) or perform other control tasks in an optimal way Zlotnik et al. (2013). In this paper we will discuss if and in what sense Eq. (1) can be used for more general dynamics than limit cycle oscillators. The key is to note that phase in Eq. (1) has the dimension of time and evolves as time in an unperturbed system Freitas et al. (2018). Thus, instead of interpreting phase as a geometric angle-like variable we can reinterpret phase as a time isomorphism $\varphi=\varphi(t)\in\mathbb{R}$ defined by Eq. (1) which parameterizes a typical trajectory $\vec{x}_{0}(\varphi)$ on a hyperbolic attractor. Indeed, in the following we will adopt the viewpoint that phase is time in the unperturbed system, i.e.

[TABLE]

For stable limit cycle oscillators the distance between a perturbed trajectory $\vec{x}(t)$ with

[TABLE]

and the phase shifted unperturbed trajectory $\vec{x}_{0}(\varphi(t))$ is bounded by $O(\varepsilon)$ for all times and arbitrary perturbations. Then Eq. (1) with $\vec{Z}(\varphi)=\vec{Z}(\vec{x}_{0}(\varphi))$ and $\vec{p}(\varphi,t)=\vec{p}(\vec{x}_{0}(\varphi),t)$ predicts the phase velocity in linear order of $\varepsilon$. Throughout the paper we assume $\vec{x}_{0}=\vec{x}_{0}(\varphi)$ to be a solution of the unperturbed system (2) evolving on an invariant set, e.g. a limit cycle or a chaotic attractor. Vector fields, such as the phase sensitivity can be expressed as functions of space or of time $\vec{Z}=\vec{Z}(\vec{x}_{0})=\vec{Z}(\varphi)$ with respect to the points of the trajectory. Note, that the scalar $\varepsilon$ in (3) quantifies to the linear order the strength of any perturbation. Such a perturbation does not need to be additive but can be applied to a system parameter, as well. E.g. with $\vec{f}=\vec{f}(\vec{x},\mu)$ and $\mu=\mu_{0}+\varepsilon\Delta\mu$ (3) takes the form $\dot{\vec{x}}=\vec{f}(\vec{x},\mu_{0})+\varepsilon\Delta\mu\partial_{\mu}\vec{f}(\vec{x},\mu_{0})$.

In Sec. II we review the classic experimental and numerical methods to obtain the phase sensitivity for autonomous limit cycle oscillations. In Sec. III we generalize these methods to hyperbolic chaotic oscillators. We show in Sec.III.1 how our re-interpretation of phase as time in the unperturbed system can be used to define phase sensitivity from the frequency response of an oscillatory system. In Sec. III.2 we improve on a well established linear least squares method Schwabedal et al. (2012) to define approximate Isochrons for chaotic oscillators. The main contribution of this paper in Sec. III.3 is the proposal to use covariant Lyapunov vectors Ginelli et al. (2007) to define the phase sensitivity function for hyperbolic chaotic oscillators. We test this proposal in numerical examples in Sec. IV.

II Phase response functions of limit cycle oscillators

There are three common and equivalent approaches to obtain the linear PRCs of autonomous limit cycle oscillators as described in the works of Winfree, Kuramoto and Malkin. These are based respectively (i) on asymptotic phase or time shifts caused by single impulses at a prescribed phase Winfree (1967), (ii) on calculating the gradient of isochrons which are parameterized by the periodic phase Kuramoto (2003), and (iii) on the solution to an adjoint linearized equation Malkin (1949); Izhikevich (2007). The direction and amplitude of $\vec{Z}(\vec{x}_{0})$ at a point $\vec{x}_{0}$ of a limit cycle follow from two geometric considerations : $\vec{Z}$ must be perpendicular to the stable invariant manifold since perturbations on this manifold do not lead to phase shifts. Secondly, a perturbation in the direction of the flow advances the phase by an amount inversely proportional to the flow velocity, i.e. $\vec{Z}(\vec{x})\cdot\vec{f}(\vec{x})=1$.

II.1 Measuring time shifts

The first method is an experimental approach and requires no mathematical model of the system dynamics. Deviations from a stable limit cycle caused by a small, single pulsed perturbation $\vec{p}=\Delta\vec{x}\delta(t-t_{0})$ at a phase $\varphi_{0}=\varphi(t_{0}^{-})$ decay exponentially fast. The instantaneous phase shift $\varphi(t_{0}^{+})-\varphi(t_{0}^{-})=\Delta\varphi=\varepsilon\vec{Z}(\varphi_{0})\cdot\Delta\vec{x}$ according to Eq. (1) remains constant afterwards and can be measured as a permanent time shift between the perturbed and an unperturbed system signal. $Z_{\Delta}(\varphi_{0})=\lim_{\varepsilon\to 0}\Delta\varphi/\varepsilon$ is a phase response function. The index $\Delta$ stands for any experimentally realizable pulsed perturbation, either in the dynamic variables or in the system parameters. E.g., kicking the system in a single component of the state variable $\vec{x}$, i.e. replacing $x$ by $x+\varepsilon\Delta x$, will result in a time shift $\Delta\varphi/\varepsilon\to Z_{x}\Delta x=Z_{\Delta}$. In control problems the system state may not be directly accessible and a system parameter $\mu$ may only vary within practical limits. In this case it is impossible to apply a delta kick and a localized parametric forcing $\mu=\mu_{0}+\varepsilon\Delta\mu(\varphi)$ over a finite time interval $[\varphi_{0}-\tau,\varphi_{0}+\tau]$ and finite strength must be applied which results in a time shift

[TABLE]

II.2 Isochrons

Isochrons (or isophases) $I_{\varphi}$ are invariant manifolds under the system propagation over one oscillation period $T_{0}$. They intersect the limit cycle in one point $\vec{x}_{0}(\varphi)$ which is an attracting fixed point of the time $T_{0}$ forward map on $I_{\varphi}$. All points on an isochron have the same phase $\varphi(I_{\varphi})=\varphi(\vec{x}_{0})$ and the same phase velocity $\dot{\varphi}(I_{\varphi})=\dot{\varphi}(\vec{x}_{0})=1$. Thus phase is defined everywhere in the basin of attraction of the limit cycle as a scalar field $\varphi=\varphi(\vec{x})$. Phase response is not restricted to small perturbations of a system close to the limit cycle Kuramoto (2003); Wilson and Ermentrout (2018). The phase sensitivity is given as the gradient $\vec{Z}(\vec{x})=\vec{\nabla}\varphi(\vec{x})$, which is orthogonal to the isochrons and with $\dot{\varphi}=1$ follows $\dot{\varphi}=\vec{\nabla}\varphi\cdot\dot{\vec{x}}=\vec{Z}\cdot\vec{f}=1$ everywhere.

II.3 Malkin’s adjoint method

Malkin’s method considers deviations from the limit cycle only to the linear order. Here isochrons are linear subspaces in the tangent space at each point $\vec{x}_{0}$ of the limit cycle. Vectors $\vec{h}$ in the tangent space evolve under the periodic action of the system Jacobian matrix $\textrm{J}_{f}(\varphi)=\textrm{J}_{f}(\vec{x}_{0}(\varphi))$ with $(\textrm{J}_{f})_{ij}=\partial f_{i}/\partial x_{j}$ along the limit cycle as $d\vec{h}/d\varphi=\textrm{J}_{f}\vec{h}$ . Invariance under system propagation over one period means that an invariant subspace is spanned by Floquet vectors. Perturbations in the stable directions do not change the phase, whereas perturbations in the direction $\vec{f}$ of the flow do not decay. The co-vectorfield $\vec{Z}(\vec{x}_{0})$ which is the unique solution of the adjoint linear equation

[TABLE]

on the limit cycle normalized to $\vec{Z}\cdot\vec{f}=1$ is orthogonal to the stable invariant subspace (see Sec.III.4) and therefore equal to the linear phase sensitivity Malkin (1949); Izhikevich (2007). Malkin’s adjoint method is the standard way to obtain the phase sensitivity numerically, when the linearization $\textrm{J}_{f}(\varphi)$ of the dynamics at the limit cycle $\vec{x}_{0}(\varphi)$ is available.

III Phase response for chaotic oscillators

Since the discovery of chaotic phase synchronization Rosenblum et al. (1996) many heuristic approaches have been suggested to generalize phase response theory to autonomous chaotic oscillators and to define PRCs or phase coupling functions Pikovsky et al. (1997); Josić and Mar (2001); Beck and Josić (2003); Schwabedal et al. (2012); Rusin et al. (2011); Kurebayashi et al. (2014). The main difficulty is, that due to mixing and chaotic phase diffusion, usually no globally differentiable isochrons exist in chaotic oscillators. Phase shifts caused by perturbations are not asymptotically constant and can therefore not be measured in a unique way. All three methods must and can be modified if one wants to apply them to chaotic oscillators.

In the following we will distinguish time like phase $\varphi$ from an angle like geometric phase $\vartheta(\vec{x})$ which parameterizes a periodic foliation of the state space into Poincare sections $P_{\vartheta}=P_{\vartheta+2\pi}$ and is increasing monotonously ($d\vartheta/d\varphi>0$) along a trajectory $\vec{x}_{0}(\varphi)$. In general a geometric phase $\vartheta_{0}$, e.g. reconstructed from a time series by Hilbert-transform or some other embedding technique, does not evolve uniformly. Such a geometric phase is called a protophase. For limit cycles a simple rescaling from an arbitrary protophase $\vartheta_{0}$ to a uniformly evolving geometric phase $\vartheta$ is always possible Kralemann et al. (2008).

III.1 Measuring frequency response

Eq. (1) describes a time isomorphism $\varphi=\varphi(t)$. Conversely, time as a function of phase evolves to the linear order in $\varepsilon$ as

[TABLE]

For perturbations $\vec{p}=\vec{p}(\varphi)=\vec{p}(\vec{x}_{0}(\varphi))$ without explicit time dependence we can take the average over $\varphi$ corresponding to an average along an unperturbed trajectory $\vec{x}_{0}(\varphi)$ and obtain

[TABLE]

where $T_{0}=1/\nu_{0}$ and $T_{\varepsilon}=1/\nu_{\varepsilon}$ are the average periods, $\nu_{0}$ and $\nu_{\varepsilon}$ the frequencies of the unperturbed and of the perturbed system, respectively. Instead of measuring an asymptotic time shift caused by a single perturbation pulse it is also possible to measure the shift of the average oscillation period or frequency in linear response to a perturbation that only depends on the position on the attractor

[TABLE]

The phase sensitivity may be expanded into a set of vector fields $\vec{Z}=\sum_{k}z_{k}\vec{p}_{k}$ which are orthonormal under the scalar product on the left hand side such that

[TABLE]

Or the system is kicked with $\vec{p}=\Delta\vec{x}\sum_{i}\delta(t-t_{i})$ everytime $t_{i}$ a Poincare section $P_{\vartheta}$ is crossed after one oscillation. Then from Eq. (1) follows that the average PRC on that Poincare section is

[TABLE]

if the limit exists, i.e. the chaotic system does have a linear response to the perturbation $\vec{p}$. All propositions for a phase sensitivity $\vec{Z}$ must be judged by comparing the predicted frequency shifts to measurements. The works Rusin et al. (2011); Kurebayashi et al. (2014) use frequency response to define such average or effective PRCs on Poincare sections $P_{\vartheta}$ constructed from the $T_{0}$ forward map.

III.2 Optimizing a geometric phase

A possible heuristic approach is to define isochrons as a family of Poincare sections $P_{\vartheta}$ parameterized by a geometric phase $\vartheta\in[0,2\pi)$ and optimize these surfaces under a set of constraints such that the variance of the return time is minimized Schwabedal et al. (2012). Here, instead of the time domain, we perform the optimization in the frequency domain which has some advantages, as we will see. As in Schwabedal et al. (2012) we expand a geometric phase $\vartheta_{\sigma}(\vec{x})$ around a proto-phase $\vartheta_{0}(\vec{x})$ in the neighborhood of the attractor into an appropriate set of non-constant, differentiable functions $q_{k}(\vec{x})$

[TABLE]

e.g. Laguerre polynomials and spherical harmonics in spherical coordinates or Taylor polynomials and Fourier components in cylindrical coordinates. Given the vector fields $\vec{v}^{(l)}(\vec{x}_{0})$ in the stable, unstable and $\vec{v}^{(0)}=\vec{f}$ neutrally stable directions on the attractor, we require the gradient $\vec{\nabla}\vartheta_{\sigma}$ to be orthogonal to the stable and unstable directions $\vec{\nabla}\vartheta\cdot\vec{v}^{(l\neq 0)}\approx 0$ and $\vec{\nabla}\vartheta\cdot\vec{f}=\dot{\vartheta}\approx\omega_{0}$. Indeed, such a vector field $\vec{Z}\parallel\vec{\nabla}{\vartheta}$ exists and is uniquely determined by the vectors $\vec{v}^{(l)}$. It can be used as phase sensitivity function in some sense, as we will discuss in the next section. As a finite sum of differentiable functions $\vartheta_{0}$ and $q_{k}$, the gradient $\vec{\nabla}\vartheta_{\sigma}$ of the geometric phase $\vartheta_{\sigma}$ is a differentiable approximation of $\vec{Z}$ such that $\vec{\nabla}\vartheta_{\sigma}\approx\omega_{0}\vec{Z}$ and

[TABLE]

Applying the method of linear least squares to Eq. (12) the coefficients $\sigma_{k}$ can be found which minimize the square norm of the deviations $\eta_{l}(\vec{x}_{0})$ over all points on the attractor. Choosing $l=0$, i.e. $\vec{v}^{(0)}=\vec{f}$, we can include points, and calculate $\vec{f}$ there, which are close but not exactly on the attractor. Then $\vec{\nabla}\vartheta_{\sigma}\cdot\vec{f}=\dot{\vartheta}_{\sigma}\approx\omega_{0}$ will evolve approximately uniformly in the neighborhood of the attractor. The advantages over the method Schwabedal et al. (2012) of Schwabedal et al. are that the phase velocity of points which are not on the attractor can easily be calculated in contrast to the return times and that we can include additional linear constraints if the stable and unstable directions are available. Note that the lengths of the vectors $\vec{v}^{(l)}$ with $l\neq 0$ are arbitrary. Choosing them, e.g. in some relation to the flow velocity $\vec{f}$ makes Eq. (12) a weighted linear least squares problem. Secondly, including all Lyapunov vectors in the linear least squares problem essentially results in a smooth geometric phase with a phase gradient that on the attractor approximates the theoretical phase sensitivity $\frac{1}{\omega_{0}}\vec{\nabla}\vartheta\approx\vec{Z}$.

III.3 Using co-variant Lyapunov vectors

Measuring the frequency response will give some approximation of $\vec{Z}(\vec{x}_{0})$, or rather projections of $\vec{Z}$ on the chosen perturbations. However, the response of chaotic systems is often not differentiable. In Fig.1a we show the mean period as a function of system parameter $b$ in the chaotic Roessler oscillator (see Sec.IV.2). At bifurcation points where the system attractor changes non differentiably the mean period is also not differentiable. Then the limits Eq. (8) and Eq. (10) (with $\varepsilon\sim\Delta b$) may not exist and the frequency measurements can give contradicting results for different $\varepsilon$. For uniformly hyperbolic chaos, on the other hand, linear response has been proven Ruelle (1997). In this class of systems the stable, unstable and neutrally stable manifolds intersect in each point $\vec{x}_{0}$ of the attractor and are nowhere tangential. The tangent space at $\vec{x}_{0}$ is spanned by the co-variant Lyapunov vectors $\vec{v}^{(k)}(\vec{x}_{0})$ in the stable, unstable and neutrally stable directions. These vector fields on a hyperbolic chaotic attractor have their correspondence in the Floquet vectors on a limit cycle. In the unperturbed system a small deviation $\vec{h}(\varphi)$ from a tajectory $\vec{x}_{0}(\varphi)$ is evolved by the aperiodic Jacobian matrix $\textrm{J}_{f}(\varphi)=\textrm{J}_{f}(\vec{x}_{0}(\varphi))$

[TABLE]

In a co-moving (co-variant) base of Lyapunov vectors, the dynamics of the components $h_{k}$ of $\vec{h}(\varphi)=\sum_{k}h_{k}\vec{v}^{(k)}$ decouple as

[TABLE]

The $\lambda^{(k)}$ are local Lyapunv exponents and the averages $\Lambda^{(k)}=\langle\lambda^{(k)}\rangle_{\varphi}$ are the Lyapunov exponents on the system attractor. If $\vec{f}=\vec{f}(\vec{x})$ is time independent, one Lyapunov exponent $\Lambda^{(0)}=\lambda^{(0)}=0$ is zero, globally and locally, and the corresponding Lyapunov vectorfield is $\vec{v}^{(0)}=\vec{f}$. Shifts $h_{0}$ in the direction of $\vec{f}$ result in time shifts which do not grow or decay. The vector field $\vec{Z}(\vec{x}_{0})$ which is orthogonal to the Lyapunov vectors in the stable and unstable directions and is normalized to $\vec{Z}\cdot\vec{f}=1$ is the obvious generalization of the phase sensitivity function to chaotic oscillators (see Fig.2). However, any perturbation with components in the unstable directions will lead to an exponentially growing deviation of a perturbed trajectory from the unperturbed trajectory. The distance between these trajectories is not bounded as $O(\varepsilon)$, they have no well defined phase relationship $\varphi(t)$ and shifts in any geometric angle-like phase due to pulsed perturbations are not asymptotically constant. Nevertheless, $\vec{Z}(\vec{x}_{0})$ does have all necessary properties for a phase sensitivity function for one particular perturbed trajectory, which depends on the perturbation $\vec{p}(\vec{x},t)$, and shadows the unperturbed trajectory $\vec{x}_{0}(\varphi)$.

Phase response conjecture for shadowing trajectories: Given a trajectory ${\vec{x}}_{0}(\varphi)$ on a uniformly hyperbolic invariant set of a flow generated by a dynamics $d{\vec{x}}_{0}/d\varphi=\vec{f}({\vec{x}}_{0})$, and without any other continuous symmetries than time-shift invariance, for any sufficiently small perturbation $\varepsilon\vec{p}({\vec{x}},t)$ of arbitrary but finite duration, i.e. $|\vec{p}({\vec{x}},t)|=0$ for $t\notin[t_{0},t_{0}+\tau]$, there exists a unique time isomorphism $\varphi=\varphi(t)$ with $\varphi(t_{0})=t_{0}$ and a unique $\varepsilon$-close trajectory ${\vec{x}}_{\varepsilon}(t)$ such that $d{\vec{x}}_{\varepsilon}/dt=\vec{f}({\vec{x}}_{\varepsilon})+\varepsilon\vec{p}({\vec{x}}_{\varepsilon},t)$ holds exactly and $\lim_{t\to\pm\infty}|{\vec{x}}_{\varepsilon}(t)-{\vec{x}}_{0}(\varphi(t))|=0$. The time derivative of $\varphi$ in linear order of $\varepsilon$ is given by

[TABLE]

where the phase sensitivity function $\vec{Z}({\vec{x}}_{0})$ is the unique vector field orthogonal to the stable and unstable manifolds at ${\vec{x}}_{0}$ and normalized to $\vec{Z}({\vec{x}}_{0})\cdot\vec{f}({\vec{x}}_{0})=1$. $\Box$

Equation (15) defines to the linear order the instantaneous time-shift of the shadowing trajectory relative to the unperturbed trajectory for arbitrary perturbations. After the perturbation is switched off the shadowing trajectory $\vec{x}_{\varepsilon}$ will converge to the unperturbed trajectory with an accumulated asymptotic phase shift $\Delta\varphi=\varphi(t)-t$. A mathematical proof of the existence of a shadowing trajectory for flows and equivalence of Lipschitz boundedness of the shadow to structural stability was given in Pilyugin (1997); Pilyugin and Tikhomirov (2010). In our conjecture, by imposing the boundary condition $\varphi(t_{0})=t_{0}$ and requiring asymptotic convergence of the shadow to the unperturbed trajectory in both temporal directions the isomorphism $\varphi=\varphi(t)$ and the shadow ${\vec{x}}_{\varepsilon}(t)$ are defined uniquely. Moreover, using co-variant Lyapunov vectors Ginelli et al. (2007), the phase and the shadowing trajectory can be constructed explicitly in linear order of $\varepsilon$. The conjecture is also valid for structurally stable invariant sets of non-hyperbolic dynamics, i.e. unstable periodic orbits (UPOs) embedded into a non-hyperbolic chaotic attractor. Phase sensitivity of UPOs has been used in Pikovsky et al. (1997) to study chaotic phase synchronization. In Fig.1b we demonstrate linear phase response by constructing the shadow of the period-1 UPO in the chaotic Roessler oscillator under periodic forcing. We chose a forcing amplitude and frequency such that the shadowing trajectory is synchronized to the forcing. Equation (15) has the same significance as Eq. (1) for periodic oscillators; it is a nonlinear equation for the phase dynamics based on linear response theory, expressing the effect of a perturbation as a product of a phase sensitivity function and the perturbation itself. This makes it, for instance, possible to use linear methods to construct perturbations that optimize the response for some purpose Kiss et al. (2007); Zlotnik et al. (2013). Using the method of linear least squares from Sec.III.2 it is possible to construct a differentiable geometric phase $\vartheta_{\sigma}(\vec{x})$ which approximates $\vec{Z}$ on the attractor as $\vec{\nabla}\vartheta_{\sigma}\approx\omega_{0}\vec{Z}(\vec{x})$.

III.4 Construction of the shadow trajectory

Let us consider a solution $\vec{x}_{0}(\varphi)$ of an autonomous dynamics Eq. (2) on a hyperbolic attractor with Jacobian matrix $(\textrm{J}_{f})_{ij}=\partial f_{i}/\partial x_{j}$ and a small deviation $\vec{h}(\varphi)=\sum_{k}\left(\vec{u}^{(k)}\cdot\vec{h}\right)\vec{v}^{(k)}=\sum_{k}h_{k}\vec{v}^{(k)}$ where the $\vec{v}^{(k)}=\vec{v}^{(k)}\left(\vec{x}_{0}\right)$ and the $\vec{u}^{(k)}=\vec{u}^{(k)}\left(\vec{x}_{0}\right)$ are co-moving bases of biorthonormal Lyapunov vectors and co-vectors following the equations

[TABLE]

The $\lambda^{(k)}=\lambda^{(k)}(\varphi)\in\mathbb{R}$ are local Lyapunov exponents of the system. Equations (16,17) conserve biorthonormality $\vec{u}^{(k)}\cdot\vec{v}^{(l)}=\delta_{kl}$ along a trajectory while the $\lambda^{(k)}$ on average compensate for the expansion or contraction in the directions of the Lyapunov vectors Ginelli et al. (2007). Equation (16) for $\lambda^{(0)}=0$ is trivially solved by $\vec{v}^{(0)}=\vec{f}$, the Lyapunov vector field corresponding to the neutrally stable direction of the flow. We now consider the evolution of $\vec{x}_{\varepsilon}(t)=\vec{x}_{0}(\varphi(t))+\vec{h}(\varphi(t))$ in a perturbed system

[TABLE]

Here we have introduced the isomorphism $\varphi=\varphi(t)$. To the linear order of $\varepsilon$ and $|\vec{h}|$ we have

[TABLE]

Multiplying Eq. (19) by $\vec{u}^{(k)}$ from the left, using biorthonormality, $\frac{d}{d\varphi}\vec{h}=\sum_{k}h_{k}\frac{d}{d\varphi}\vec{v}^{(k)}+\vec{v}^{(k)}\frac{d}{d\varphi}h_{k}$ and Eq. (16) we obtain

[TABLE]

and

[TABLE]

For the correct isomorphism $t=t(\varphi)$ the perturbed trajectory $\vec{x}(t(\varphi))$ is always contained in the subspace spanned by the stable and unstable directions at $\vec{x}_{0}(\varphi)$, i.e. $dh_{0}/dt=0$. Therefore

[TABLE]

Let the perturbation be of finite but arbitrary long duration, i.e. $|\vec{p}(\vec{x},t)|=0$ for $t\notin[t_{0},t_{0}+\tau]$. Then the sufficient conditions for convergence of the perturbed trajectory to the unperturbed trajectory for perturbations of any form are $h_{0}=0$, $h_{k}(t_{0})=0$ for $\Lambda^{(k)}<0$ and $h_{k}(t_{0}+\tau)=0$ for $\Lambda^{(k)}>0$. In other words, the perturbed trajectory is fully contained in the unstable manifold at the beginning of the perturbation, hence convergence for $t\to-\infty$, and fully contained in the stable manifold at the end of the perturbation such that the shadow converges to the unperturbed trajectory for $t\to\infty$ (see Fig.2). Using these as initial and final conditions we can integrate Eqs. (20,22) forward in time, beginning at $t_{0}$ for the components of $\vec{h}$ in the stable directions and backward in time beginning at $t_{0}+\tau$ for the components in the unstable directions. Equation (20) being linear, the distance $|\vec{h}|$ of the perturbed trajectory is always of order $\varepsilon$ on uniformly hyperbolic invariant sets, i.e. when the dynamics in the stable and unstable directions is uniformly contracting or expanding. Furthermore denoting $\vec{Z}=\vec{u}^{(0)}$ Eq. (22) is identified as Eq. (6) and can in linear order of $\varepsilon$ be rewritten as Eq. (15). It is, however, easier to integrate Eq. (22) when $\vec{x}_{0}(\varphi)$ is given at discrete time points $\varphi_{i}$. Using $\varphi(t_{0})=t_{0}$ as initial condition makes the time isomorphism unique. The phase sensitivity function $\vec{Z}(\vec{x}_{0})$ is the unique vector field solving the adjoint equation $d\vec{Z}/d\varphi=-\textrm{J}_{f}^{\top}\vec{Z}$, i.e. Eq. (17) for $k=0$, on the hyperbolic attractor normalized to $\vec{Z}\cdot\vec{f}=1$.

III.5 Discussion

This method of generating a shadowing trajectory is of equivalent accuracy as a recently proposed linear least squares method Wang et al. (2014) which, however, cannot reproduce the correct time isomorphism or the phase sensitivity. It is reasonable to assume linear response theory is valid for perturbations $\varepsilon\vec{p}$ that lead to small distances $|\vec{h}|$ in Eq. (20). The more the stable and unstable directions $\vec{v}^{(k)}$ align, the larger the Lyapunov co-vectors $\vec{u}^{(k)}$ become. This puts practical limits on the perturbation strength. If $|\vec{h}|\leq h_{\textrm{max}}$ for a given perturbation $\varepsilon\vec{p}$, then $h_{\textrm{max}}$ depends linearly on $\varepsilon$. In particular for non-uniformly hyperbolic systems $h_{\textrm{max}}$ may occasionally become very large even for small $\varepsilon$. Furthermore, while the trajectory $\vec{x}_{\varepsilon}(t)$ is an exact solution of the perturbed system, it may however not be a typical solution, i.e. time averages are not necessarily equal to averages with respect to the natural invariant measure on the perturbed attractor. If the chaotic attractor is not structurally stable then in the vicinity of a larger bifurcation, e.g. a periodic window, the attractor, and thus the oscillation period, can change discontinuously in response to the perturbation (see Fig. 1a). The measured frequency response Eqs. (8,10) may only approximately be predicted by the projection of the perturbation on $\vec{Z}$, i.e. by the averaged response of the UPOs embedded in the chaotic attractor, which are not close to a bifurcation.

Although many analytic results are valid for hyperbolic systems, physical examples of hyperbolic chaotic flows are rare Kuznetsov and Pikovsky (2007). On the other hand, the algorithm Ginelli et al. (2007) for the numerical determination of the Lyapunov vectors is quite robust against occasional near tangencies of stable and unstable manifolds along the trajectories on the chaotic attractor. Even for non-hyperbolic systems such as the Roessler system for small enough perturbations one can construct shadowing trajectories which remain close to an unperturbed trajectory for periods of time longer than expected from the largest rate of divergence given by the largest Lyapunov exponent. In Sections IV.1 and IV.2 we present examples of non-hyperbolic chaotic oscillators where our method can reliably predict the frequency response. In Section IV.3 we show that our method works with a known example of hyperbolic chaotic oscillations and in IV.4 we discuss why our method works poorly in the non-hyperbolic Lorenz system.

IV Examples

IV.1 Electrochemical oscillations

As an example we consider current oscillations during the electro-dissolution of a metal in an acidic environment. A mathematical model for such electro-chemical oscillations, which exhibits a period doubling route to non-hyperbolic chaos, was developed in Koper and Gaspard (1992) and used in Kiss and Hudson (2001) to reproduce in simulations the experimentally observed chaotic current oscillations through a nickel electrode in sulfuric acid. After an appropriate re-scaling we obtain

[TABLE]

with nonlinearity

[TABLE]

The applied Voltage $V_{a}$ and the electrode potential drop $E$ can be measured and controlled. At the parameters $V_{a}=36.7380$, $R_{s}=0.02$ and $d=0.119$ the system attractor has developed two chaotic bands around an unstable period-two orbit. The applied voltage $V_{a}$ needs to be controlled precisely since the region of chaotic oscillations in parameter space is very small. Only the phase sensitivity in the $E$ component is of experimental interest since $U$ and $W$ quantify a gradient of chemical concentrations in the solution (double layer approximation) and cannot be measured. However, for the computation of the Lyapunov vectors the full knowledge of the system state, velocity and Jacobian are assumed. We have calculated the Lyapunov exponents on the chaotic attractor as $(\Lambda^{(0)},\Lambda^{(+)},\Lambda^{(-)})=(0,0.07,-2.5)$. For Lyapunov vectors $\vec{v}^{(\pm)}$ of unit length, the local Lyapunov exponents $\lambda^{(\pm)}$ exhibit large excursions to both positive and negative values (Fig.3e). As a consequence the deviations of a shadowing trajectory in Eq. (20) can become quite large, even for small perturbations.

We define $x=E$, $y=E(t-0.3)$ and $z=E(t-0.6)$. A geometric proto-phase $\vartheta_{0}$ with positive phase velocity and an amplitude $R$ can be defined via $x=R\cos\vartheta_{0}$ and $y=R\sin\vartheta_{0}$. The phase $\vartheta=\vartheta_{\sigma}(\vartheta_{0},R,z)$ that we want to optimize is expanded into

[TABLE]

with

[TABLE]

$\sigma^{+}_{000}=0$ and $\sigma^{-}_{0lm}=0$. The choice of the cutoff values for the Fourier harmonics $k$ and polynomial orders $l$, $m$ depends on the particular geometry of a system. Lower values avoid over-fitting with large deviations at points that are not on the attractor, whereas larger values can give better results for the points on the attractor. Since the stable and unstable directions in the time delayed coordinates are not known we only use Eq. (12) with $l=0$ and $\vec{f}=\frac{d}{d\varphi}(x,y,z)$ to minimize the variance of the deviations $\eta_{0}$ in

[TABLE]

The delay embedding of the chaotic attractor with color-coded optimized phase is shown in Fig.3a. The velocity of the proto-phase and of the optimized phase as functions of $\vartheta$ are shown in Fig.3c. In both cases the mean phase velocity is $\omega_{0}=3.001$ but the standard deviation of the optimized phase velocity is at $0.03$ within 1% of $\omega_{0}$.

Next we perform a series of perturbation experiments. A small delta kick in the applied potential $V_{a}$ is executed after each full rotation when the system crosses the Poincare section $P_{\vartheta}$ at a given optimized geometric phase $\vartheta$ in the delay coordinates. The measured shift in the average period according to Eq. (10) gives the PRC $Z_{\Delta}(\vartheta)$ (Fig.3d, small markers). This PRC can be compared with the components $Z_{1}$ of $\vec{Z}(\vartheta)$ in the $E$ direction. Here $\vec{Z}=\vec{Z}(\vec{x}_{0})$ is calculated numerically from the co-variant Lyapunov vectors Ginelli et al. (2007). As a function of $\vartheta(\vec{x}_{0})$ the values of $Z_{1}(\vec{x}_{0})$ form a family of curves, shown as thin blue lines in Fig.3d. An average response $\langle Z_{1}\rangle(\vartheta)$ is calculated via a narrow Gaussian filtering of the data points $Z_{1}(\vartheta)$ (red squares). For this chaotic oscillator our linear frequency response theory predicts the measured PRC $Z_{\Delta}(\vartheta)$ very well. However the strength of the delta kicks must be very small $\left(\Delta E<5\times 10^{-4}\right)$ in order to approximately retain the structure of the chaotic attractor, and even smaller $\Delta E\approx 1\times 10^{-4}$ for a shadow trajectory which in linear order of the perturbation stays within an acceptable small distance to the unperturbed trajectory (Fig.3f).

IV.2 Roessler oscillator

The chaotic Roessler oscillator is often used as an example for chaotic phase synchronization Rosenblum et al. (1996); Pikovsky et al. (2003); Rosa et al. (1998). Chaotic phase diffusion in the Roessler system is extremely small Rosa et al. (1998), which facilitates phase synchronization in this system. Tangencies of the Lyapunov vectors occur but the directions of the co-variant Lyapunov vectors are well separated most of the time Ginelli et al. (2007). The transition to chaos occurs via period doubling with frequent bifurcations of the attractor where the response is not differentiable (Fig.1a). However, these structural changes in the attractor can be small if the main UPOs are not close to a bifurcation. We study the chaotic Roessler oscillator with a logarithmic variable $z=\exp(q)$, effectively making additive perturbations in $q$ multiplicative in $z$, ensuring that $z$ remains positive. The dynamics in these variables is given by

[TABLE]

where we have used $a=0.25$, $b=0.9$ and $c=6.0$. Unusually large values of the phase sensitivity $\vec{Z}(\varphi)$ (Lyapunov co-vector) do occur which have a strong influence on the average phase response. In the averages $\langle Z_{i}\rangle(\vartheta)$ we therefore disregard values of the phase sensitivity larger than three standard deviations. A protophase $\vartheta_{0}$ and radial distance $R$ for this system is defined as $x=R\cos\vartheta_{0}$ and $y=R\sin\vartheta_{0}$. For the optimized phase we use the same expansion and cutoff as in the previous example of the electrochemical oscillator (Sec.IV.1). However, we determine the optimized phase using the full information of the unit length Lyapunov vectors $\vec{v}^{\pm}$ and the flow direction $\vec{v}^{0}=\vec{f}$. The resulting optimized phase $\vartheta(\vec{x})$ is then used in the perturbation experiments to measure the average shift in the rotation period in response to delta kicks at the crossing of a given Poincare section $P_{\vartheta}$ after each full rotation. We compare the predicted average response by the method of Lyapunov vectors to the measured response to delta kicks of strength $\varepsilon=0.05$ in Fig.4.

IV.3 Hyperbolic chaotic oscillations

We will now demonstrate our theory at the following example of hyperbolic chaotic dynamics Kuznetsov and Pikovsky (2007)

[TABLE]

The amplitudes $a_{1}^{2}=x_{1}^{2}+y_{1}^{2}$ and $a_{2}^{2}=x_{2}^{2}+y_{2}^{2}$ of two oscillators are coupled via a negative feedback loop where the first oscillator acts as an activator and the second as an inhibitor leading to sequential switching between low and high amplitude oscillations. Through weak forcing with $\kappa=0.3$ the phase of the lower amplitude oscillator synchronizes to the phase of the high amplitude oscillator. By coupling the first oscillator to the second harmonics of the second oscillation via the product $\kappa x_{2}y_{2}$ the phases of the oscillators after each round of switching are chaotic following an expanding circle map. The system’s Lyapunov exponents are $\lambda\in\left\{-1.34,-0.97,0,0.09\right\}$. Because the amplitudes can become very small, for numerical stability we simulate (33-36) using angle and log-amplitude variables $\psi$ and $q$, i.e. $x+iy=\exp(q+i\psi)$. The phase sensitivity has components $\vec{Z}=(Z_{q_{1}},Z_{q_{2}},Z_{\psi_{1}},Z_{\psi_{2}})$ in these variables corresponding to delta Kicks in the log-amplitudes and angles or $\vec{p}_{q}=(x,y)\delta(t-t_{0})$ and $\vec{p}_{\psi}=(-y,x)\delta(t-t_{0})$ in the original variables. As Poincare sections we define the sets of geometric angles $\vartheta_{0}$ with $a_{1}^{2}-\left\langle a_{1}^{2}\right\rangle=R\cos\vartheta_{0}$ and $a_{2}^{2}-\left\langle a_{2}^{2}\right\rangle=R\sin\vartheta_{0}$. We re-parametrize these angles $\vartheta_{0}\to\vartheta$ such that $\vartheta$ is uniformly distributed over $[0,2\pi)$. This is achieved by defining $\vartheta(\vartheta_{0})$ linearly increasing with the rank of the protophases sorted over the points of the attractor. Optimization of the shapes of the Poincare sections is not necessary. In Fig.5a we show a projection of the hyperbolic chaotic attractor in the $(q_{1},q_{2})$ plane and $(a_{1}^{2},a_{2}^{2})$ in the inset. The switching dynamics can be seen in Fig.5b where $x_{1}$ and $x_{2}$ are plotted as a function of $\vartheta$. The mapping of the angle $\psi_{1}$ of the first oscillator from one crossing of the Poincare section $\vartheta_{0}=0$ to the next is shown in Fig.5c. It follows an expanding circle map. In Fig.5e,f we show the components $Z_{q_{1}}$ and $Z_{\psi_{1}}$ of the phase sensitivity function, i.e. the respective components of the Lyapunov co-vector $\vec{u}^{(0)}$ with $\vec{u}^{(0)}\cdot\vec{f}=1$, as a function of $\vartheta$, corresponding to phase shifts from perturbations in the log-amplitude or the angle of the first oscillator. Next we measured the frequency shift caused by pulsed perturbations in the log-amplitude or the angle of the first oscillator as a function of $\vartheta$. That is, after each full oscillation of the system at a Poincare section at a given geometric phase $\vartheta$ the log-amplitude or the angle of the first oscillator was increased by $\varepsilon=0.1$ and the mean period was determined by the elapsed time between $1000$ crossings. The the frequency response function at $\vartheta$ is then calculated as $Z_{\Delta}(\vartheta)=\frac{1}{\varepsilon}\left(T_{0}-T_{\varepsilon}\right).$ This measure is also shown in Fig.5e,f and it traces the correpsonding components of the Lyapunov co-vector field $\vec{Z}$ very well. Finally, by integrating (20) with kicked log-amplitude $q_{1}$ perturbations of strength $\varepsilon=0.01$ at crossings of $\vartheta=0$ we have calculated the displacements $\vec{h}(\vartheta)$ of the shadowing trajectory and plotted $h=|\vec{h}|$ as a function of $\vartheta$ in Fig.5d.

IV.4 Lorenz System

Finally we present phase and frequency response in the non-hyperbolic chaotic Lorenz system

[TABLE]

We use $\sigma=16$, $\beta=2.0$ and $\rho=28$ where the system is chaotic with Lyapunov exponents $\Lambda^{+}=0.8$ and $\Lambda^{-}=-20$. The proto-phase is defined by $R\cos\vartheta_{0}=z-z_{0}$ and $R\sin\vartheta_{0}=\sqrt{x^{2}+y^{2}}-u_{0}$ with respect to the fixed point coordinates $u_{0}=\sqrt{2\beta(\rho-1)}$ and $z_{0}=\rho-1$. We expand the optimized phase around that proto-phase as

[TABLE]

with

[TABLE]

$\sigma^{+}_{00}=0$ and $\sigma^{-}_{0l}=0$ and find the coefficients $\sigma^{\pm}_{ml}$ which minimize the variance of the phase velocity. In Fig.6a we show a projection of the Lorenz attractor to the $(x,z)$ coordinates. The points are colored according to 100 intervals of the optimized geometric phase. Red and blue shades, respectively, signify positive and negative average frequency response to perturbations in the $z$ direction, predicted by the $z$-component of the phase sensitivity $\vec{Z}$. In Fig.6b we project the chaotic oscillations to the coordinates $\left(\sqrt{x^{2}+y^{2}},z\right)$ used in the definition of the proto-phase. In numerical experiments we have performed 1000 kicked perturbations $\varepsilon\vec{p}=0.5\vec{e}_{z}\delta(t-t_{kick})$ at constant optimized geometric phase after each oscillation and measure the resulting shift of the oscillation period. In Fig.6c we compare the frequency response in the perturbation experiments to the average phase sensitivity at that geometric phase predicted by the $z$-component of $\vec{Z}$. Shown are the average phase sensitivity after convolution of $Z_{z}$ with a narrow Gaussian (black curve) and the average phase sensitivity restricted to values within three standard deviations (yellow curve). Apparently large deviations in $Z_{z}$ have a strong influence on the predicted average response. The frequency response measured in the perturbation experiments (blue crosses) follow in parts the features of both averages but can also deviate significantly from the predictions. The standard deviation $\textrm{std}(Z_{z})=0.17$ is ten times larger than the actual response, and the extreme values seem to follow a power-law over two orders of magnitude (Figs.6d and 7b). The reason for this are frequent near tangencies of the unstable Lyapunov direction and the flow, which can be seen in Fig.7a from the distribution of angles (blue histogram). The Lorenz system is an example of a non-hyperbolic chaotic oscillator where our method performs poorly.

V Conclusions

Measuring the frequency response to pulsed perturbations at a given Poincare section is a simple and experimentally viable way to define and measure phase response functions of chaotic oscillators. In this work we have presented a theoretical approach to predict these frequency shifts with the help of co-variant Lyapunov vectors. A phase sensitivity $\vec{Z}=\vec{Z}(\vec{x}_{0})$ can be constructed for the points on the chaotic attractor. Time shifts along a chaotic trajectory in response to arbitrary perturbations can be calculated to the linear order of the perturbation strength in the same way as for limit cycle oscillators with Winfree type phase equations Eq. (1). These time shifts are only exact for a certain perturbed trajectory shadowing the unperturbed trajectory. However, averaging the time shifts for time independent perturbations over the whole attractor can approximate the frequency shift for arbitrary perturbed trajectories. Given the phase sensitivity $\vec{Z}(\vec{x}_{0})$ a differentiable geometric phase $\vartheta(\vec{x}_{0})$ can be constructed with a gradient $\vec{\nabla}\vartheta$ which approximates the phase sensitivity and minimizes the variance of the phase velocity on the attractor and in its vicinity. We demonstrate our theory with a chaotic electro-chemical oscillator and the chaotic Roessler oscillator, both examples of non-hyperbolic, i.e. non structurally stable systems, where the numerically determined Lyapunov vectors can give good approximations of the linear frequency response. Because of large deviations in the Lyapunov co-vector field $\vec{Z}(\vec{x}_{0})$ Frequency response in the non-hyperbolic chaotic Lorenz system is not well predicted. We have also included an example of hyperbolic autonomous oscillations, where the Lyapunov vectors and the phase sensitivity $\vec{Z}(\vec{x})$ can be determined numerically robustly.

Acknowledgements.

We thank Z. Arai, H. Nakao, A. Pikovsky and K. Takeuchi for valuable discussions. H.K. acknowledges the financial support from MEXT KAKENHI Grant No. 15H05876 and JSPS KAKENHI Grant No. 18K11464.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Appendix A

We use the method developed by Ginelli et al. Ginelli et al. (2007) to determine the co-variant Lyapunov vectors $\vec{v}^{(k)}(\varphi)$ along a chaotic trajectory $\vec{x}_{0}(\varphi)$ evolving according to Eq. (2) on the system attractor. Both, the time step $d\varphi$ forward map $\vec{M}(\vec{x}_{0},d\varphi)=\vec{x}_{0}(\varphi+d\varphi)$

[TABLE]

and its Jacobian matrix $\textrm{J}_{M}$ with

[TABLE]

are integrated simultaneously by standard RK4 fourth order Runge-Kutta method. The Lyapunov vectors (except $\vec{v}^{(0)}=\vec{f}$) are normalized $|\vec{v}^{(k)}|=1$ so that $d\vec{v}^{(k)}/d\varphi$ and $\vec{v}^{(k)}$ are orthogonal. With that and from Eq. (16) follow the local Lyapunov exponents

[TABLE]

Convergence of the Lyapunov vectors means independence from initial conditions in both forward and backward integration. Choosing two different random initial matrices of Lyapunov vectors, convergence to the co-variant Lyapunov vectors can be monitored. Given the matrix $\textrm{V}=(\vec{f},\vec{v}^{(1)},\vec{v}^{(2)},\dots)$ of co-variant Lyapunov vectors, the matrix $\textrm{U}=(\vec{Z},\vec{u}^{(1)},\vec{u}^{(2)},\dots)$ of co-variant Lyapunov co-vectors is simply the inverse Matrix of V, i.e.

[TABLE]

The Lyapunov co-vectors do not have unit length. Because of biorthonormality, alignment of the Lyapunov vectors, brings V closer to degeneracy and results in large Lyapunov co-vectors. Where Lyapunov vectors, and thus stable, neutrally stable and unstable subspaces become tangential, U, and $\vec{Z}$ in particular, is divergent.