Canonical-dissipative limit cycle oscillators with a short-range interaction in phase space

We consider limit cycle oscillators in terms of canonical-dissipative systems that exhibit a short-range interaction in velocity and position space as described by the Dirac delta function. We derive analytical expressions for stationary distribution functions in phase space and energy space and propose a numerical simulation scheme for the simulation of the many body system as well. We show that the short-range interaction squeezes or stretches energy distribution functions depending on whether the interaction can be regarded as attractive or repulsive. In addition to the interaction effect, we show that energy distribution functions become narrower when limit cycle attractors become stronger. Finally, energy distributions become broader when the pumping energy is increased. The latter effect however disappears in the high energy domain.


Introduction
Classical Newtonian and Hamiltonian mechanics is a fundamental topic in physics.With the increasing number of interdisciplinary applications of physics the question arises to what extent the concepts of classical mechanics can be generalized for more general classes of systems as addressed by classical mechanics.In this context, the concept of canonical-dissipative systems has been developed to address self-propagating systems, self-excited systems, and in particular limitcycle oscillators that are supplied by energy sources [1][2][3][4][5].Several studies have demonstrated that the canonical-dissipative approach can yield helpful insights into this kind of open systems [6][7][8][9][10][11][12][13][14][15][16][17][18].
Let us briefly review some of the key properties of a typical canonical-dissipative system [1-5] that will be essential for our present study.While the transient dynamics of a canonical-dissipative system differs qualitatively from the dynamics of a classical Newtonian and Hamiltonian system, in the long term limit the trajectory of a canonical-dissipative system converges to the trajectory of a classical Newtonian and Hamiltonian system.The reason for this behavior is that the invariants of motions of the corresponding classical, dissipation-free Hamiltonian system are not constant for the canonical-dissipative system.Rather, these invariants evolve as functions of time during the aforementioned transient period and only in the long term limit converge to particular stationary values.The evolution of the invariants is caused by several possible environmental impacts in terms of energy, matter, or information supply.Taking a dynamic systems perspective, these impacts may be simply classified into two categories: damping and pumping.That is, we may say that the evolution of the invariants is determined by the interplay of the damping and pumping forces that reflect the interaction of the system under consideration with its environment.In the long time limit, a balance between damping and pumping is established, the invariants become constant and the canonical-dissipative system evolves on a trajectory which is consistent with a classical Hamiltonian system.The canonical-dissipative approach is not only concerned with the deterministic evolution of open dynamic systems but also is concerned with the relationship between damping forces and fluctuating forces.In doing so, stochastic system properties can be addressed.For example, distribution functions of the Hamiltonian energy of canonical-dissipative systems in general exhibits one or several peaks at energy values larger than the minimal energy [1,4,5,8,10,17,[19][20][21][22][23].
The objective of the present study is to examine such energy distributions for oscillatory canonical-dissipative systems that are subjected to a short-range interaction described by the Dirac delta function.Such a short-range interaction has been proposed by Haissinski [24,25] and has been studied extensively in accelerator physics [26][27][28][29][30][31][32][33][34][35][36][37][38].Note that recently it has been shown that the Tsallis entropy [39,40] is also closely related to this type of short-range interaction [41].From previous research it is known that due to the short-range interaction, systems can become more concentrated or less concentrated at the points in phase space that are most probably occupied.In other words, the short-range interaction is known to be attractive or repulsive and is known to squeeze or stretch distributions [24,34].In section 2.1 we define an ensemble of canonicaldissipative oscillators subjected to the aforementioned short-range interaction.Analytical results for the distribution function of the oscillators will be derived in sections 2.2 and 2.3.They will be compared with numerical results in section 3. The impact of the interaction and the model parameters related to pumping and damping will be discussed in section 4. Some generalizations will be addressed in section 5.

The oscillator model
To begin with, we consider a single oscillator.The oscillator is assumed to evolve in one spatial dimension and exhibits a Hamiltonian function H(x, v), where x is the position (or elongation) and v the velocity of the oscillator.The Hamiltonian oscillator equations read where I = (I x , I v ) denotes the conservative drift vector and t is time.Let P (x, v, t) denote the probability density of the oscillator.Then, P satisfies the Liouville equation with ∇ = (∂/∂x, ∂/∂v).Taking dissipation and fluctuations into account we consider a free energy functional F [P ] that determines the evolution of the probability density P (x, v, t) such that the Liouville equation becomes a Fokker-Planck equation of the form [5,42] with γ > 0. The expression δF/δP denotes the variational derivative of the functional F with respect to the probability density P .We assume that the free energy is bounded from below.It can be shown that in this case the probability density P as a function of time t converges to a stationary function.That is, the underlying Markov process related to equation (4) becomes stationary (H-theorem) [5,43,44].For the free energy functional F we may put where S BGS denotes the Boltzmann-Gibbs-Shannon entropy or information measure defined by S BGS = − ln P (here and in what follows we put the Boltzmann constant equal to unity) [45].
For oscillators in thermodynamic equilibrium, T reflects the oscillator temperature.In our context of an oscillator operating far from thermodynamic equilibrium, we consider T as an effective oscillator temperature related to the strength of an external fluctuating force acting on the oscillator (we return to this issue below).Substituting equation ( 5) into equation ( 4), we can show that trajectories x(t), v(t) are determined by the Langevin equation The function Γ(t) is a Langevin force normalized like Γ(t)Γ(t ) = 2δ(t − t ), where δ(•) denotes the Dirac delta function.We obtain Brownian motion for the special case of a free particle with H = v 2 /2 (and unit mass).Due to damping, the Hamiltonian energy of the oscillator decays as a function of time.In order to study a self-excited oscillator exhibiting a limit cycle, we need to consider an oscillator pumped by an energy source.We study energy-uptake in the framework of canonical dissipative oscillators and replace equation ( 5) by where E is a constant.In this case, the underlying Markov process exhibiting the free energy ( 8) is defined by the Langevin equation which can be shown by substituting equation ( 8) into equation ( 4).In the deterministic case (T = 0), from equations ( 9) and ( 10) it follows that Consequently, for E > 0 the Hamiltonian H as a function of time t converges to the constant E.
We refer to E as oscillator fixed point energy.In the stochastic case (i.e. for T > 0) the stationary distribution of the Hamiltonian energy P st (H) exhibits a maximum at the fixed point energy E.
Next we consider an ensemble of identical canonical-dissipative oscillators.We label individual oscillators by the index k.That is, we assume that equations ( 9) and (10) describe the evolution of oscillators k = 1, . . ., N when they are isolated from each other: The oscillator distribution ρ is formally defined by We assume that there is an interaction between the oscillators that can be expressed in terms of the energy with For κ > 0 the expression κJ/2 can be interpreted as a heat energy with J corresponding to the Tsallis entropy S T,2 [39] with nonextensivity index 2: J = S T,2 − 1.In general (i.e. for κ < 0 and κ > 0), the term κJ/2 can be interpreted as mean field energy emerging from a local (i.e.short-range) interaction W = κδ(x − x )δ(v − v ) in the two-dimensional phase space spanned by x and v.As mentioned in the introduction, this type of interaction has been proposed by Haissinski [25].In the original Haissinski model the parameter κ describes attractive interactions for κ < 0 (distributions are squeezed) and repulsive interactions for κ > 0 (distributions are stretched) [24,33].In line with the mean field theory [46][47][48], we consider a large ensemble of oscillators (i.e. the limiting case N → ∞) and replace the oscillator distribution ρ by the distribution P of an individual representative oscillator such that equation (15) becomes Taking the interaction terms ( 14) and ( 17) into account, the free energy (8) of an individual oscillator becomes Substitution of equation ( 18) into the free energy Fokker-Planck equation ( 4) yields a closed description for the canonical-dissipative oscillators under consideration.We obtain with κ = γκ and Q = γT .Since the interaction term J is quadratic in P , the closed description given by equation ( 19) is nonlinear with respect to P and exhibits the form of a so-called nonlinear Fokker-Planck equation [5,42,46,47].The oscillator trajectories are determined by a self-consistent Langevin equation [5] We see that Q can be interpreted as strength of the external fluctuating force acting on the oscillator.In the literature, the parameter Q is also referred to as noise amplitude.
In what follows, we will study in detail the interacting oscillators when H corresponds to the Hamiltonian of a harmonic oscillator with angular oscillation frequency ω and unit mass: In this case, the nonlinear Fokker-Planck equation (19) reads and the Langevin equations ( 20) and (21) We will briefly address Hamiltonian functions different from equation ( 22) in section 5.

Stationary probability density in the phase space
In the stationary case ∂P st (x, v)/∂t = 0, we let We have (see equation (23)) Further, we try to solve the equation (26).Expanding the second term by taking partial derivative with respect to v the equation becomes Considering the partial derivative of f (H) with respect to x and v, we have and We put these terms in the equation ( 27) and rearrange, so we have After that, we integrate both sides and choose the constant of integration equal to zero (because of natural boundary conditions [44]).Then equation (30) becomes We expand the equation ( 31) and use equation (29).Then, we get the equation that depends on H as follows: Let H = H − E so that f (H) = f ( H). Equation ( 32) becomes Comparing equation (33) with the classical Haissinski equation [31,33], the solution of equation (33) can be found in the form where the function g( H) is In equation ( 34) the function LW(•) denotes the Lambert-W function [49] and z > 0 corresponds to an integration constant.From f( H) we obtain the original function f (H) as After back substitution in terms of x,v,H,and E, we get and 13001-5 where the constant z satisfies That is, we will interpret z in what follows as a normalization constant.In a more concise form, the stationary probability density P st (x, v) defined by equations ( 37) and ( 38) can be written as

Stationary probability density of the Hamiltonian energy
From the probability density P st (x, v) we can compute the probability densityP st (H) which satisfies ∞ 0 P st (H)dH = 1.We have We put y = ωx and dx = dy/ω: Next we introduce ỹ = y/ √ 2 and ṽ = v/ √ 2 and transform ỹ and ṽ into polar coordinates r and ϕ.We obtain (v 2 + y 2 )/2 = ṽ2 + ỹ2 = r 2 , dydv = 2dỹdṽ = 2rdrdϕ and We use s = r 2 , ds = 2rdr and obtain That is, P st (H) reads explicitly

Limiting case of vanishing oscillator-oscillator interaction
Recall that the Lambert-W function LW(•) is implicitly defined by LW(x) exp{LW(x)} = x and we have LW(0) = 0, see [49].Differentiating the implicit definition of LW with respect to x yields d LW(x)/dx[1 + LW(x)] exp{LW(x)} = 1.This implies that at x = 0 the slope equals 1: d LW(0)/dx = 1.Having said that we can consider the limiting case of a vanishing interaction: κ → 0. In this case, the numerator and denominator in equations ( 40) and ( 46) vanish.Using the rule of de l'Hopital, we differentiate numerator and denominator with respect to κ and consider the limiting case κ → 0 again.For example, for equation (40) we calculate the limiting case of vanishing oscillator-oscillator interaction as follows: Likewise, for P st (H) we obtain with z = ωz/(2π).Equations ( 47) and ( 48) recover the well-known results for canonical-oscillators for J = 0 (i.e., without interaction) [1][2][3][4][5].However, these limiting cases not only hold for noninteracting oscillators.They also hold when the fluctuating force dominates the interaction: κ/Q → 0. That is, repeating the same calculations as before, we obtain and Note that these limiting cases include the cases considered in equations ( 47) and ( 48) as special cases but hold for more general situations.For example, we may consider oscillator systems with the amount of interaction between the oscillators as measured by κ does not vanish.

Boundaries for the critical value of the parameter κ
For systems involving attractive interactions (i.e., interaction parameters κ < 0), it is known that stationary probability densities exist only up to a critical parameter value κ c equation [31,34,35,38,41].In this section, we determine an interval in which this critical value can be found.As we will see, the precise value of the critical parameter satisfies an implicit equation that has to be solved numerically.
In this section, we consider the oscillator system from the perspective of the Hamiltonian energy H, again.We start from equation (32), which we write in the form: We need to consider the case κ < 0, for which equation ( 51) reads As shown in [31,[33][34][35]41], positive definite, continuous and normalizable solutions of equations of the form (52) only exist if Q − |κ|f (H) > 0 holds for all H. Furthermore, in this case f is maximal at Consequently, we define the critical parameter κ c by which implies f M = Q/|κ c |.By separating variables and using f M = Q/|κ c |, equation ( 52) for κ = κ c can be written in the form Integrating the right hand side from E to an arbitrary value H and accordingly the left hand side from f M to f (H), we get From equation ( 56) we obtain The minus sign holds for the monotonously increasing part of the graph f (H) on the interval [0, E], where we have H − E 0. Likewise, the plus sign holds for the decreasing part of f (H) for H ∈ [E, ∞), where we have H − E 0. Substituting equation (57) into equation (55), we eliminate the bracket (H − E) in equation ( 55) and obtain Next, we exploit the normalization condition We split the integral into two parts with E * = E for E 0 and E * = 0 otherwise.That is, for E 0 the first integral vanishes.Substituting dH (see equation ( 58)) into equation (60), yields Note that for E 0 we have f M = f (0) (see section 4.1) and the first integral vanishes (which is consistent with the comment on equation (60) made above).Note also that for the first integral we have used the minus sign occurring in equation (58) because the first integral concerns the interval [0, Using Both integrals yield positive numbers.Recall that the critical parameter is negative.In the limit E → ∞ we have f (0) → 0 and consequently both integrals approach the same value.In this case, we have the smallest possible critical value with In contrast, for E 0 the first integral vanishes and we have the maximal critical value For E > 0 we therefore can identify the interval in which the critical value of the parameter κ can be found as κ c ∈ [2κ * , κ * ].Alternatively, we say that the interval [2κ * , κ * ] contains all critical values for a fixed set of parameters ω, γ, and Q when all possible values E ∈ lR are considered.The precise value of κ c for a given energy E can be determined numerically by solving equation (64) when taking into account that both f (0) and f M depend on κ c .That is, κ c satisfies the implicit equation As mentioned above, for large enough parameters E the probability to find energies H = 0 becomes negligibly small and κ c equals 2κ * .A numerical value for the integral in equation ( 66) is 1.0964 such that for large parameters E we have κ c ≈ −2.2πω −1 2Q 3 /γ.

Numerical simulation scheme for the stationary case
In order to solve equations ( 24) and ( 25) numerically for the stationary case, we first realize that in the stationary case the probability density P st (x, v) becomes a function of the Hamiltonian H. Consequently, we have Substituting this result into equations ( 24) and ( 25), we obtain with H(t) = [v 2 (t) + ω 2 x 2 (t)]/2.The Langevin equation involves the gradient of the probability density P st (H).We compute this gradient analytically rather than numerically.From our discussion in section 2.2 it follows that Consequently, we compute the derivative like We discretize time on a grid with time step ∆t.Using an Euler forward method [44] that also holds for self-consistent Langevin equations [5], we transform equations (70, . . ., 72) into with In the simulation scheme above, e i is a Gaussian random number with variance 2. As indicated in equation (77), the value P i denotes the function value of the probability density of P st (H) at H = H i .In order to determine P st (H) numerically, we compute simultaneously an ensemble of N oscillators.That is, we solve equations (74, . . ., 77) simultaneously for N subsystems.At every iteration step we compute from the N observed H i values the probability density P st (H) using the MATLAB kernel density estimator.This implies for example that the velocities v i+1 of all N oscillators depend on the same function P st (H).However, for every individual oscillator the function P st (H) is evaluated at the appropriated, individual Hamiltonian value of that oscillator.We checked the consistency of the numerical simulation scheme with the analytical result derived in section 2.2.To this end, we plotted the P st (H) as computed from equation ( 46) and as obtained numerically in the same figure.Figures 1 and 2 show the results obtained for attractive (κ < 0) and repulsive (κ > 0) interactions between the oscillators.

Impacts of model parameters
We will study the impact of the model parameters on the distribution P st (H) of the Hamiltonian energy H as defined by equation (46).The Lambert-W-function LW(•) is a monotonously increasing function.Consequently, from equation ( 46) it follows that P st (H) is unimodal (i.e.exhibits only one maximum or peak).Let H peak and P max denote the position and the height of this peak, respectively.At issue is to study how the model parameters affect the peak position and height.To this end, we consider the functions H peak = H peak (E, γ, Q, κ, ω) and P max = P max (E, γ, Q, κ, ω).From equation (46) it is obvious that the distribution function P st (H) only depends on the ratios γ/Q and κ/Q such that it is sufficient to consider the functions H peak = H peak (E, γ/Q, κ/Q, ω) and P max = P max (E, γ/Q, κ/Q, ω).

Impact of the net pumping force described by the fixed point energy E
For E > 0 the pumping force exceeds the friction force such that the net pumping force is positive.For E < 0 the friction dominates the pumping.From equation ( 46) it follows that the peak position depends exclusively on the fixed point energy E. We have This relationship is illustrated in figure 3.In addition, figure 4 shows two distributions P st (H) computed from equation ( 46) for E < 0 and attractive (κ < 0) as well as repulsive (κ > 0) interaction between the oscillators.Figure 4 clearly demonstrates that for E < 0 the maximum can be found at the minimal possible energy H = 0.
The probability maximum P max in general depends on the fixed point energy E. Since P st (H) can be cast into the form P st (H) = r(H − E), where r(•) is a function of the difference H − E, we can regard E as a parameter that shifts the coordinate system of H.If E increases then we shift the distribution from the left to the right.When considering P st (H) as a mathematical function (not as a probability density) then we can plot P st (H) on the whole real line H ∈ lR.However, when considering P st (H) as a probability density, then the area under P st (H) for H 0 must be equal to 1.If we increase E and shift the mathematical function P st (H) for fixed normalization constant z from the left to the right along the coordinate axis H, then the area under the mathematical functions on the semi-line H 0 increases.Consequently, if we increase E and adjust the normalization constant z such that we are dealing with a normalized probability 13001-11 density defined on H 0, then the maximum of the P st (H) decreases.We have For relatively large values of E we are dealing with functions as illustrated in figures 1 and 2. The probability to find oscillators with H = 0 is relatively small.In the limit E → ∞ this probability can be neglected and a shift of E does not change the normalization constant z.That is, we have ∂z/∂E = 0 for E → ∞ and lim

The relative interaction parameter κ/Q
The ratio κ/Q can be interpreted as relative interaction parameter.It measures the strength of the interaction relative to the strength of the fluctuating force.The ratio κ/Q assumes positive values for repulsive interactions and negative values for attractive interactions.As mentioned above the peak position H peak does not depend on κ/Q.However, κ/Q crucially affects the concentration of the oscillators energies H around the most probable energy H peak .That is, P max depends on κ/Q. Figure 5 illustrates the function P max (•, κ/Q, •) for fixed values of E, γ/Q and ω.We see that P max decreases monotonously.In other words, for repulsive interaction (κ > 0) an increase of the relative interaction strength makes the distribution of H broader and results in a decrease of the peak value P max .In contrast, for attractive interactions (κ < 0), increasing the amount of the relative interaction strength |κ|/Q implies that the attractive interaction becomes stronger and the oscillator energies become more concentrated at H peak .The peak height P max increases and the energy distribution P st (H) becomes narrower.(46).Fixed parameters: E = 5, Q = 2, ω = 3.For the graph with γ/Q = 0.5 (γ = 1) the parameter κ was varied in the range of −9, . . ., 25 such that κ/Q varied in the range of −4.5, . . ., 12.5.For the graph γ/Q = 1.5 (γ = 3) the parameter κ was varied in the range of −2.5, . . ., 12.5.(11) it is clear that in the deterministic case the parameter γ determines how fast the Hamiltonian energy H relaxes to the fixed point energy E. Therefore, we may interpret γ as a measure for the attractor strength of the limit cycle of an individual canonical-dissipative oscillator.Note that the relaxation time depends not only on the parameter γ but also on other parameters involved in the function H.However, if we divide both sides of equation (11) with γ and substitute t with a new time variable t = γ t then we see that γ defines a time scale relevant for the relaxation dynamics.Therefore, it is sensible to consider γ as a key relaxation parameter.

13001-12
Furthermore, recall that attractor strength usually is defined in terms of relaxation times towards attractors.In our case, it is plausible to develop the notion that the attractor strength of the limit cycle of an individual oscillator is proportional to γ.The larger γ, the stronger is the limit cycle attractor.The ratio γ/Q describes the relative limit cycle attractor strength and can be used to answer the question: how strong is the attraction towards the limit-cycles relative to the erratic perturbations that push the oscillator away from the limit cycle and are due to the fluctuating force with amplitude Q?
As discussed in section 4.1 the parameter γ rel = γ/Q does not affect the peak position H peak .From equation (46) it follows that if γ/Q increases then the width of distribution becomes smaller.Consequently, the peak value P max increases monotonously as a function of γ rel = γ/Q: This relationship is illustrated in figure 6.Note also that figure 5 illustrates that the peak value increases as a function of γ rel .For a fixed value of κ/Q we can compare in figure 5 the values of P max for γ rel = 0.5 and γ rel = 1.5.We see that P max (γ rel = 1.5, κ/Q) > P max (γ rel = 0.5, κ/Q) for fixed values of ω and E.  (46).Fixed parameters: E = 5, Q = 2, ω = 3.For both graphs κ/Q = 0.5 (κ = 1) and κ/Q = 1.5 (κ = 3) the parameter κ was varied in the range of 1, . . ., 25 such that κ/Q varied in the range of 0.5, . . ., 12.5.

Stationary phase space probability densities
In section 2.2 we derived analytical expressions for the stationary probability densities P st (x, v) and P st (H) for canonical-dissipative oscillators involving the Hamiltonian function ( 22) of a harmonic oscillator.We consider next Hamiltonian functions that involve a potential function V (x) which is continuous and globally attractive V (x → ±∞) → ∞ but does not necessarily correspond to a parabolic potential.That is, we put A detailed analysis of the analytical steps carried out in section 2.2 shows that the derivation of P st (x, v) holds for all kinds of potentials V (x) and does not exploit the explicit form V (x) = ω 2 x 2 /2 considered in section 2.2.Consequently, P st (x, v) can be expressed in terms of the function f (H) (defined by equation ( 36)): In particular, the stationary distribution P st (x, v) for oscillators involving the Hamiltonian function ( 82) is again defined by equation (40).

Stationary distributions of Hamiltonian energy
In order to derive the stationary distribution P st (H) of the Hamiltonian energy H, we need to determine the proportionality factor β that relates P st (H) to f (H).For the harmonic potential V (x) = ω 2 x 2 /2 we found in section 2.2 that β = 2π/ω and P st (H) = 2πf (H)/ω.That is, we found that for parabolic potentials β is independent of H.However, in general the factor β can depend on the Hamiltonian energy H.
In order to determine β(H) and P st (H), we proceed as in section 2.2.We substitute H = v 2 /2 + V into f (H) (defined by equation (36)) such that f = f (v 2 /2 + V (x)) and subsequently solve the integral We first evaluate the integration with respect to v. To this end, we define the integral Recall that for the Dirac delta function that involves an arbitrary continuously differentiable function u the relation holds, where v i are the roots of u: u(v i ) = 0. From equation (84 which can equivalently be expressed as From equations (83), (84), and (87) we conclude that As indicated in equation ( 89 In the context of classical Hamiltonian systems, the integral bm am ( 2[H − V (x)]) −1 dx is wellknown.Recall that from the so-called first integral v = dx/dt = ± 2([H − V (x)] we obtain

13001-14
Interacting canonical-dissipative oscillators dt = ( 2[H − V (x)]) −1 dx with the plus sign in the case of a motion for which x(t) increases as function of time (i.e. for v > 0).Therefore, the integral bm am ( 2[H − V (x)]) −1 dx corresponds to the duration T am→bm of a particle with unit mass that travels from a m to b m and is subjected to a classical Hamiltonian dynamics.We define Using T am↔bm (H), we can express the equations (89, . . ., 91) as follows: From equation (93) it is clear that the factor β reflects the sum of all oscillation periods (roundtrip times) for oscillations around those potential minima of V (x) that exhibit potential energies smaller than H.Each individual oscillation period and the sum of all periods depend in general on H. Consequently, as anticipated in the beginning of this section, the factor β in general depends on H.
In closing this section, let us return to the harmonic oscillator potential V (x) discussed in section 2.2.For V = ω 2 x 2 /2 we have only one potential minimum (i.e.we put M = 1 in equation ( 93)) and the oscillation period T (a ↔ b) is independent of H and simply corresponds to the inverse of the oscillation frequency f osc = ω/(2π).We have T (a ↔ b) = 2π/ω.Consequently, the factor β is independent of H and we recover the result β = 2π/ω derived in section 2.2.

Conclusions
Limit cycle oscillators are at the heart of nonlinear physics and play a key role in nonlinear sciences.Unfortunately, the derivation of analytical results for stochastic limit cycle oscillators is usually mathematically involved.Canonical-dissipative limit cycle oscillators provide an exception to this rule because they frequently can be studied analytically.In view of this property of canonicaldissipative oscillators, we proposed to derive an analytical expression for the distribution function of canonical-dissipative limit cycle oscillators that are coupled with each other by means of a short-range interaction in phase space.
We discussed canonical-dissipative oscillators involving a parabolic potential in detail.We were able to determine the impacts of various forces by studying how the distribution function of the Hamiltonian energy is affected by these forces.That is, we demonstrated that stochastic properties of the canonical-dissipative oscillators exhibiting short-range interaction can conveniently be studied when we consider the distribution function of the Hamiltonian energy rather than the distribution function in the phase space spanned by position and velocity.We found that the peaks of energy distribution functions become more pronounced and consequently that energy distribution functions become narrower when limit cycle attractor forces and interaction forces become stronger.This holds for attractive interactions between oscillators.However, a peculiarity of the interaction considered in our study (and motivated by early work of Haissinski) is that the oscillator-oscillator interaction cannot only be attractive but also repulsive.For a repulsive interaction, we found that peaks become less pronounced and energy distributions become broader when the magnitude of the interaction is increased.The relative impact of the limit cycle attractor as indexed by the parameter γ rel is the same for attractive and repulsive inter-oscillator interactions.From this observation we conclude that repulsive inter-oscillator interactions can be counter-balanced by the intrinsic (within-oscillator) attractor forces of limit cycle attractors.In other words, when the relative magnitude κ/Q of the repulsive inter-oscillator interactions is increased and at the same time the relative limit cycle attractivity γ rel = γ/Q is increased appropriately, then the impacts of the two manipulations may cancel out each other such that there is a zero net effect.
In section 3, we presented a numerical simulation scheme which can be used to compute trajectories of individual oscillators of the oscillator ensemble.We showed that the results of the simulation scheme are consistent with the analytically derived results.At this stage, however, we should note that the simulation scheme can also be used to numerically compute the second-order statistical properties (e.g.autocorrelation functions, mean square displacements, and so on) which we have not addressed in our theoretical considerations.Consequently, the proposed model provides a complete description of the oscillator systems under consideration.That is, all statistical quantities of interest are well defined and can at least be computed numerically.
In section 5, we briefly illustrated that our considerations can be generalized to oscillator systems that involve potentials V different from parabolic potentials.We would like to point out that the simulation scheme proposed in section 3 can be modified such that we can deal with this more general situation.One can show that for systems addressed in section 5 involving potentials V (x) equations (74, . . ., 77) become with The realm of canonical-dissipative systems is much broader than the systems addressed in our study [1][2][3][4].The key quantity of canonical-dissipative systems are the invariants of classical mechanical systems featuring a Newtonian or Hamiltonian dynamics.The Hamiltonian energy is one of the most important invariants but other invariants may play important roles as well in certain applications.Therefore, the present study should consider a first step in exploring the properties of canonical-dissipative systems under various types of short-range interactions.The results summarized above and the analytical techniques developed in our work may be used as guidelines in investigating canonical-dissipative systems that focus on invariants different from the Hamiltonian energy.

Figure 3 .
Figure 3. H peak as a function of E drawn from equation (78).
), the integral with respect to x over the function w essentially defines the factor β. Let a m and b m with m = 1, . . ., M denote the positions at which the potential energy V equals H (i.e.we require V (a m ) = H and V (b m ) = H) with a 1 < b 1 < a 2 < b 2 < . . . .Then the integration with respect to x is only carried out over the intervals [a m , b m ]:

)
The round-trip time form a m to b m and back to a m is 2 times the duration for the one-way travel because classical Hamiltonian dynamics is time reversible.Let T am↔bm (H) = 2T am→bm (H) denote the round-trip time.Note that the round-trip time corresponds to the oscillation period of oscillations around a potential minimum (or several potential minima) of V (x) located in [a m , b m ].