Synchronization and anchoring of two non-harmonic canonical-dissipative oscillators via Smorodinsky-Winternitz potentials

Two non-harmonic canonical-dissipative limit cycle oscillators are considered that oscillate in one-dimensional Smorodinsky-Winternitz potentials. It is shown that the standard approach of the canonical-dissipative framework to introduce dissipative forces leads naturally to a coupling force between the oscillators that establishes synchronization. The non-harmonic character of the limit cycles in the context of anchoring, the phase difference between the synchronized oscillators, and the degree of synchronization are studied in detail.

non-harmonic case, we consider self-oscillations in Smorodinsky-Winternitz potentials [21][22][23]. These potentials play an important role in physics as confinement potentials [24,25]. In the context of the CD approach, the four-dimensional, two particle Smorodinsky-Winternitz systems should be considered as benchmark systems because they feature three invariants rather than only two. The first two invariants are the particle Hamiltonian energy functions. The third invariant is an appropriately adjusted angular momentum [21].
Recall that the standard CD oscillator in a one-dimensional space with coordinate q and momentum p is defined by [8,10,11,26] where m denotes mass, k is the spring constant and H = p 2 /(2m)+ kq 2 /2 corresponds to the Hamiltonian energy in the conservative case in which the function g is neglected. The function g describes the dissipative mechanism. The parameter γ 0 is the coupling parameter of oscillator with the dissipative mechanism. Note that −∂g/∂p = −γp(H − B)/m such that the dissipative mechanism is composed of a negative friction term (i.e., pumping mechanism) +γBp/m with pumping parameter B 0 and a nonlinear friction term −γpH/m. The amplitude dynamics can be obtained using standard techniques. We put q = A(t) exp(iωt)+c.c. with ω 2 = k/m, where A denotes the complex-valued oscillator amplitude that is related to the real-valued amplitude r(t) and the oscillator phase Here and in what follows, c.c. denotes the complex conjugate expression. Assuming that γ is a small perturbation parameter, by means of the slowly varying amplitude approximation and the rotating wave approximation [27], we obtain, in lowest order of γ, the following amplitude dynamics Note that higher order correction terms in γ can be obtained using alternative techniques (see e.g., [28]). From equation (2) it follows that dr/dt = −γω 2 r[r 2 − 2B/(mω 2 )]/4 and dφ/dt = 0. The solution reads r(t) = 2Br 2 0 /[mω 2 r 2 0 + (2B − mω 2 r 2 0 ) exp(−γBt/m)] with r 0 = r(t = 0). Figure 1 (a) shows a simulation of the self-oscillator (2) and the analytical solution r(t). From equations (2) and the analytical solution r(t) it follows that γ (in combination with the factor B/m) determines the time scale of the amplitude dynamics A(t) and r(t), respectively. Importantly, in the long time limit, r approaches r st = 2B/(mω 2 ) and H converges to the pumping parameter B. That is, B acts as a fixed point value or target value for the energy dynamics H(t).
Let us generalize the single-oscillator case to a model of two coupled self-oscillators. Each selfoscillator oscillates in a one-dimensional space and is subjected to the force of a Smorodinsky-Winternitz potential. Let us describe the oscillator coordinates q 1 and q 2 and momenta p 1 and p 2 by means of the vectors q = (q 1 , q 2 ) and p = (p 1 , p 2 ). The Smorodinsky-Winternitz potentials read [21] with k > 0 and α, β 0. For α = β = 0, the potentials correspond to parabolic potentials and k can be interpreted as a spring constant. For α, β > 0, the potentials exhibit minima at q 1 = ±(α/k) 1/4 and q 2 = ±(β/k) 1/4 and exhibit repulsive singularities at q 1 = 0 and q 2 = 0. By contrast, for q 1 → ±∞ and q 2 → ±∞ they increase like parabolic potentials. Therefore, for α, β > 0, the potentials are asymmetric with respect to their minima. Let us define the dynamics of the two oscillators in the conservative case by means of the Hamiltonian dynamics In equation (4), the functions H 1 , H 2 and H tot correspond to the Hamiltonian energy functions of the individual oscillators and the total energy of the two-oscillators system. It can be shown that the 44001-2  (14) numerically (stochastic EF [29]). Averages of 10 trials are shown. Trajectories of 10000 (circles) and 30000 (squares) time units were used in each trial. The maximal cross correlation scores for γ 3 = 0 are by-chance values that decay to zero when trajectory length goes to infinity. Parameters in a.u.: D = 0.02, all other parameters except for γ 3 as in panels (d) and (e). τ = 0.001. dynamics (4) exhibits three invariants S j with j = 1, 2, 3 given by [21] where L denotes the angular momentum L = p 2 q 1 − p 1 q 2 . In line with the CD oscillator (1), we define the CD case like where γ j 0 are the coupling constants. B 1,2 0 are the pumping parameters. B 3 is a target value for S 3 . Importantly, since S j are invariants of the conservative dynamics (4), for the dissipative dynamics (6) it follows that dg tot /dt = −(∂g tot /∂p) 2 0. In view of the boundedness of g tot (i.e., g tot 0), g tot is a Lyapunov function that becomes stationary in the long term limit. This in turn implies that ∂g tot /∂p = 0 such that equation (6) reduces to equation (4). In total, for t → ∞, the system (6) converges to an attractor that corresponds to a solution of the conservative system (4).

44001-3
Let us consider the harmonic case defined by α = β = 0. Using q k = A k (t) exp(iωt) + c.c. with ω 2 = k/m again and A k (t) = r k (t) exp[iφ k (t)]/2 and assuming that γ j are small perturbation parameters, we obtain with . In terms of the real-valued amplitudes r 1 and r 2 and the phase difference ψ = φ 1 − φ 2 , we obtain d dt with U(ψ, r 1 , r 2 ) = mωr 1 r 2 sin(ψ)(L 2 /m − B 3 ), L = mωr 1 r 2 sin(ψ), and A detailed stability analysis based on equations (8)-(10) for the case that both oscillators are excited like B 1 , B 2 > 0, shows that the target level B 3 defines the location of the aforementioned attractor. By rescaling B 3 we obtain the location parameter η = mω 2 B 3 / (4B 1 B 2 ). It can be shown that for B 3 0 ⇒ η 0, the two-oscillators system exhibits a stable attractor characterized by S 3 = 0 ⇒ ψ = 0 • ∨ ψ = 180 • and H 1,2 = B 1,2 . This implies that g tot → γ 3 B 2 3 /2 > 0 for t → ∞. For B 3 0, we distinguish between two cases. If η ∈ [0, 1], then the limit cycle attractor is characterized by H 1,2 = B 1,2 again but with S 3 = B 3 . The latter relation implies that the phase difference is given by ψ = arcsin(  Figure 1 (b) illustrates the attractor location in terms of the phase difference ψ as a function of the location parameter η. Figure 1 (c) shows, for a representative simulation, the trajectories q 1 and q 2 . The trajectories demonstrate that the two-oscillators system converges to a stable periodic pattern (limit cycle).
Let us consider the non-harmonic case with α, β > 0. First we note that the individual oscillators (i.e., γ 1,2,3 = 0) for small amplitudes experience a linearized force of f (q j ) = −k lin q j with k lin = 4k irrespective of α and β. Consequently, the oscillation frequency is two times the oscillation frequency of the harmonic case and the period is half the period of the harmonic case. This implies that when removing the singularities in the potentials V j by putting α = β = 0, then the oscillation frequency drops in a discontinuous fashion from 2ω 0 to ω 0 with ω 2 0 = k/m. In general, the oscillation period for the oscillators j = 1, 2 can be computed from the integral being the initial energy of oscillator j. The integration limits q min,max are the turning points defined by V j = H j (t = 0). Numerical computations show that T j is independent of H j (t = 0) and α and β. In line with the small amplitude oscillation case, we obtain T j = T 0 /2 ∀ H j (t = 0) H j,min , α, β > 0 with T 0 = 2π/ω 0 , where H j,min denote the minimal energy values H 1,min = √ αk and H 2,min = √ βk. Let us consider the case γ 1,2,3 > 0. For the sake of brevity, we consider only the case in which all three invariants of the conservative dynamics, H 1,2 and S 3 , converge to their respective target values (i.e., H 1,2 → B 1,2 and S 3 → B 3 ) such that g tot → 0. Our first objective is to show that the shapes of the oscillator limit cycles are distorted compared to the harmonic case. Panels (d) and (e) of figure 1 show the phase portraits q 1 , p 1 and q 2 , p 2 obtained from a numerical simulation. The trajectories (solid lines) converge to "egg shaped" limit cycles. The limit cycles are defined by the limit cycles of the corresponding conservative oscillators with γ 1,2,3 = 0. For the example shown in panels (d) and (e), the limit cycles of the corresponding conservative oscillators are illustrated by circles. Importantly, the limit cycles reveal an anchoring phenomenon. The dynamics slows down (more pronounced as in the harmonic case) when the oscillators swing to the right of their potential minimum locations q 1 = (α/k) 0.25 and q 2 = (β/k) 0.25 , see figures 1 (d) and (e). By contrast, the dynamics speeds up when the oscillators swing to the left of their potential minimum locations. This is because the forces are relatively weak on the right-hand sides [where V j (q j ) are approximatively parabolic potentials] and relatively strong on the left-hand sides [where V j (q j ) exhibit singularities].
Our second objective is to address the synchronization of the oscillators. On the limit cycles, the oscillators oscillate with the same oscillation frequency of ω = 2ω 0 , see above. From H 1,2 → B 1,2 and S 3 → B 3 it follows that Using two first equations of (11), the function L occurring in W (defined above) can be expressed as with m, n ∈ {0, 1}. This implies that the last equation of (11) can be written as B 3 = W(q 1 , q 2 ). The synchronized state is then described by That is, the coordinate q 2 is given by a nonlinear (implicit) mapping from q 1 to q 2 . Importantly, this mapping is stable against perturbation because the two-oscillators system is attracted to the state with g tot = 0 and S 1,2,3 = B 1,2,3 . Therefore, the two oscillators are synchronized. Note that the same argument holds in the opposite direction. Considering the second oscillator as independent oscillator, the coordinate q 1 of the first oscillator is given by a nonlinear (implicit) mapping from q 2 to q 1 . Let us illustrate the synchronization of the two oscillators by considering the CD oscillator model (6) under the impact of fluctuating forces. Using a standard approach for introducing noise terms into CD systems [6][7][8]10], equation (6) becomes where Γ j (t) are independent Langevin forces [29] normalized to 2 units. The parameter D 0 is the diffusion constant. The Langevin equation (14) exhibits a Fokker-Planck equation that can be cast into the form of a free energy Fokker-Planck equation [26]. The stationary probability density P(q 1 , q 2 , p 1 , p 2 ) can then be expressed in terms of a Boltzmann function of g tot as P = exp(−g tot /D)/Z 0 , where Z 0 is a normalization factor [19,26]. Considering γ j as small perturbation parameters, we may introduce the smallness parameter γ 0 and put γ j = c j γ 0 with c j 0. Then, P = exp(−g tot /θ)/Z 0 withg tot = g tot /γ 0 holds, where θ = D/γ 0 can be considered as a non-equilibrium temperature. In this form, the analogy to equilibrium systems becomes obvious [6]. Importantly, without coupling between the oscillators, that is, for γ 3 = 0 ⇒ c 3 = 0, we have P(q 1 , q 2 , p 1 , p 2 ) = P 1 (q 1 , p 1 )P 2 (q 2 , p 2 ) with P j = exp(−g j /D)/Z j for j = 1, 2, where Z j are normalization factors again. That is, the probability density P(q 1 , q 2 , p 1 , p 2 ) factorizes. By analogy, the transition probability density (conditional probability density) factorizes. Therefore, for γ 3 = 0, there are no cross-correlations between the oscillators at any time lag. Let us show with the help of a stochastic CD model (14) that for γ 3 > 0, the two-oscillators model exhibits a stable synchronized state. To this end, we solved numerically equation (14) for a fixed value of γ 3 and calculated the cross-correlation coefficients Corr q 1 (t), q 2 (t − τ) for different time lags τ ∈ [0, T = T 0 /2]. We determined the maximal coefficient. Subsequently, we varied γ 3 . In doing so, we 44001-5 obtain the maximal cross-correlation coefficient as a function of the coupling parameter γ 3 . Figure 1 (f) summarizes the simulation results. For γ 3 = 0, there was a finite by-chance value for the maximal crosscorrelation coefficient that decayed when the simulation duration was increased. As far as the impact of γ 3 is concerned, figure 1 (f) demonstrates that the maximal cross-correlation increased as a function of γ 3 -as predicted. This increase of the maximal cross-correlation coefficient was taken as the evidence that for γ 3 , the two oscillators were to some degree synchronized.
Future studies may focus, in particular, on the stochastic aspects of the proposed CD two-oscillators model. For example, it has been suggested to use the analytical solution for the short-time propagator to define maximum likelihood estimators that can be used to estimate the model parameters of CD systems [30]. In fact, for single CD oscillator models in a series of studies, the CD theory has been applied to the experiments on human rhythmic motor behavior and model parameters have been estimated from experimental data both in the harmonic [12][13][14] and non-harmonic case [15].