Paradoxical Stabilization of Forced Oscillations by Strong Nonlinear Friction
T. Zh. Esirkepov S. V. Bulanov





Kansai Photon Science Institute, National Institutes for Quantum and Radiological Science and Technology (QST), 8-1-7 Umemidai, Kizugawa, Kyoto 619-0215, Japan


Kansai Photon Science Institute, National Institutes for Quantum and Radiological Science and Technology (QST), 8-1-7 Umemidai, Kizugawa, Kyoto 619-0215, Japan
A. M. Prokhorov Institute of General Physics, the Russian Academy of Sciences, Vavilov street 38, 119991 Moscow, Russia


In a dissipative dynamic system driven by an oscillating force,
a strong nonlinear highly oscillatory friction force
can create a quasi-steady tug, 
which is always directed opposite to the ponderomotive force
induced due to a spatial inhomogeneity of oscillations.
When the friction-induced tug exceeds the ponderomotive force,
the friction stabilizes the system oscillations
near the maxima of the  oscillation  spatial amplitude. 










Paradoxical Stabilization of Forced Oscillations by  Strong Nonlinear Friction
    S. V. Bulanov
    January 27, 2017
==============================================================================






§ INTRODUCTION


In classical mechanics, Kapitza pendulum
or  Stephenson-Kapitza pendulum
<cit.> is
a statically unstable inverted pendulum
whose statically unstable equilibrium position
is stabilized by small fast vertical oscillations of the pivot point.
This induced stability initially described by A. Stephenson <cit.>
has been first explained by P. L. Kapitza <cit.>.
To find theoretical reason of the  induced stability
Kapitza separated the pendulum motion into fast and slow oscillations
and, by averaging out fast ones,
found the effective potential 
which has minimum at the pendulum upper position, 
in contrast to a simple pendulum.

This approach created a new 
concept of dynamic stabilization in
mechanics <cit.>.

Chelomei's pendulum provides another well known example of dynamically stabilized mechanical system 
<cit.>.
The induced stability can be vindicated 
beyond the 
framework of the
method of fast and slow motions separation
<cit.>.


In general, a spatial inhomogeneity of an oscillating driving force
creates a quasi-steady ponderomotive force
<cit.>, 
directed against the spatial gradient of the
driving force.

In the Kapitza pendulum,
this ponderomotive force acts against gravitation and
makes the upper position stable.

A dissipation dampens oscillations
around the upper equilibrium position
thus further stabilizing it.

When no other forces present besides the oscillating driving force,
as in the case of charged particle dynamics in a standing electromagnetic wave,
the ponderomotive force always repels particles from
the maxima of the 

wave spatial amplitude <cit.>. 


Here we present a general model
of a dissipative dynamic system driven by an oscillating force,
where a strong nonlinear highly oscillatory friction
creates a quasi-steady tug,
which, quite counter-intuitively,
is always directed opposite to the ponderomotive force
and exceeds it for a sufficiently strong driving force.
This leads to a seemingly paradoxical
stabilization of the system oscillations
near the maxima of the spatial amplitude of the driving force.
It differs from the Kapitza pendulum effect
in that here the stabilization factor is a nonlinear growth of the friction 
with the driver force,
which creates a tug against the ponderomotive potential.





§ MODEL


We consider a simple one-dimensional model of a forced oscillation with a strong nonlinear friction,
given by the equation


    ẍ+𝒦(ℱ) ẋ =ℱ
    ,   𝒦=νℱ^2n
    .


Here the dot denotes differentiation with respect to time;
n is a natural number.
The friction coefficient 𝒦 is a non-negative function of the oscillating driving force ℱ,


    ℱ(x,t)=f_1(x)cosω t+f_2(x)sinω t
      .



The model is motivated by the dynamics
of a charged particle in a strong electromagnetic field,

where accelerating particles lose energy and undergo a recoil
due to their emission of electromagnetic radiation <cit.>.
This causes a friction which nonlinearly grows with the
electromagnetic field strength.


Following the classical approach of Ref. <cit.>,
we assume that a solution to Eq. (<ref>)
can be represented as


    x(t)=X(t)+ξ(t)


with a slowly varying function X(t)
and a fast oscillating small addition ξ(t), |ξ|≪ |X|,
which has a zero time average,


    ⟨ξ⟩=(ω/2π)∫_0^2π/ωξ(t) dt =0.


Correspondingly,
⟨ x ⟩=X,
⟨ξ̇⟩=⟨ξ̈⟩ =0, < X>≈ X(t).

Substituting (<ref>) into Eq. (<ref>) and 
expanding the functions ℱ and 𝒦 in powers of ξ
as

    ℱ≈ F+ξ∂_X F ,    F=ℱ(X,t) ;
    
       𝒦≈ K+ξ∂_X K ,    K=𝒦(ℱ(X,t)),


we obtain


    Ẍ+ ξ̈+K Ẋ+K ξ̇+
    ξẊ∂_X K + ξξ̇∂_X K
    =F+ξ∂_X F
      ,


where ∂_X is the partial derivative with respect to X
(the first argument of ℱ).

The time derivatives ξ̈ and ξ̇
are not small, being proportional to ω^2 and ω, respectively.
They are assumed to be much greater than Ẍ and Ẋ.

The friction coefficient K defined in Eq. (<ref>)
is not necessarily small; it has a 
time-independent component

    κ = ⟨ K ⟩ >0 .


In Eq. (<ref>),
slowly varying and fast oscillating terms should cancel out separately.
Neglecting the time derivatives of X,
in the zeroth order approximation 
with respect to ξ
we find  for the fast oscillating term


    ξ̈+K ξ̇=F   .


Here X as an argument of functions F and K is assumed to be constant.

The forced oscillation solution of Eq. (<ref>)
can be cast in the form

    ξ=e^-κ t-Δ(t)∫_0^t e^κτ+Δ(τ) F(τ) dτ
    ,
    
    Δ(t) = ∫_0^t [ K(τ) - κ ] dτ
    .


The first term in the expansion 
with respect to Δ of the dependence given by Eq.  (<ref>)  


    ξ=(κ f_1-ω f_2) sinω t
    -(κ f_1+ω f_2) cosω t/ω(κ^2+ω^2)


approximates the first harmonic of the solution.


Averaging Eq. (<ref>) over time
and  taking into account that < F(X,t)>≈ 0 for nearly constant X(t), we obtain


    Ẍ+(κ+⟨ξ∂_X K⟩) Ẋ=
    ⟨ξ∂_X F ⟩ -
    ⟨ξξ̇∂_X K⟩  .


Using here the expression for ξ, given by Eq.  (<ref>),
and the definitions for F and K formulated above,
we obtain
the equation for the slowly varying function X(t)
and the average friction coefficient:


    Ẍ   +κẊ=
    -∂ _X (f_1^2+f_2^2)/4(κ^2+ω^2)
    +
    n^2κ^2 ∂ _X (f_1^2+f_2^2)/(n+1)(κ^2+ω^2)^2
    +
    κ (f_2 ∂_X f_1 - f_1 ∂_X f_2)/2ω(κ^2+ω^2)[
    2n/n+1(
    κ^2-ω^2/κ^2+ω^2)
    -1
    ]
    ,
    
       κ=
    2^-2n2nnν
    (f_1^2+f_2^2)^n
    ,


where 2nn is the binomial coefficient.
The first term on the r.h.s. of Eq. (<ref>) 
corresponds to the classical ponderomotive force <cit.>, 
modified due to the friction,
the last terms represent the tug induced by the friction.
In the case of f_2=0, 
the third term in the r.h.s. of Eq. (<ref>) vanishes.
Then the tug becomes always directed opposite to the ponderomotive force.
It can even exceed the latter in magnitude,
when

    ν > 
    2^2n (n!)^2 (n+1)^1/2/(2n)! (4n^2-n-1)^1/2ω/f_1^2n  .

In this case, trajectories of the system described by the model (<ref>)
drift to the local maximum of the driving force spatial amplitude.

If the system were non-dissipative,
the oscillations near that maximum would be destabilized
by the ponderomotive force, so that
the corresponding trajectories would drift 
against the spatial gradient of the driving force
and would eventually escape to regions of a lower spatial amplitude 
of the driving force <cit.>.

Sufficiently strong friction makes the trajectories to drift 
along the spatial gradient of the driving force,
provided that these trajectories have already got to a region with sufficiently high driving force.

This causes a seemingly paradoxical
stabilization of the oscillations
near the local maximum of the driving force spatial amplitude.





§ NUMERICAL SIMULATIONS


The effect of the stabilization of oscillations due to strong friction
is further demonstrated by the numerical integration of the model equation (<ref>)
with the friction coefficient taken in the form 𝒦=νℱ(x,t)^4
in two cases.

In the first case, shown in Fig. <ref> (a,b,c),
the driving force amplitude is bell-shaped,
ℱ(x,t)=f_0exp(-(x/l_0)^2) cosω t, 
with the width, oscillation frequency, amplitude, and  
friction factor equal to
l_0=10, ω=1, f_0=3, and ν=0.2, 
respectively.

In the second case, Fig. <ref> (d,e,f),
the driving force is spatially periodic,
ℱ(x,t)=f_0cos^2(2π x/l_0) cosω t,
l_0=10, ω=1, f_0=8, and ν=0.25.

In the case of bell-shaped driving force shown in Fig. <ref> (a,b,c),
the trajectories, starting at t=0 from locations 
where the driving force spatial amplitude is relatively weak,
exhibit several oscillations and then escape
being pushed away by the ponderomotive force.

In the region of a high spatial amplitude of the driving force,
the friction-induced tug overcomes the ponderomotive force,
therefore the trajectories started from this region drift towards the
maximum of the driving force spatial amplitude.
As they get closer to that maximum,
their drift becomes slower and their oscillation amplitude decreases.

In the case of spatially periodic driving force (see Fig. <ref> (d,e,f)),
the trajectories, initially oscillating near the 
minima of the driving force spatial amplitude
and having enough large oscillation amplitude,
reach the region of higher driving  force,
where the trajectories are caught by the friction-induced tug.
Eventually all such trajectories are reduced to 
oscillations near the maxima of the driving force spatial amplitude.

In both cases the drift is slowing down
near the maxima of the driving force spatial amplitude 
in agreement with Eq. (<ref>),
because the gradient ∂_X f
vanishes at the maximum of f.







§ LIMIT CYCLE


On a trajectory asymptotically turning into periodic oscillations
seen in Fig. <ref> (b,c,e,f),
the driving force is almost constant,
which simplifies theoretical consideration allowing 
more detail description of the driven oscillations.
This case is described by the approximation of 
ℱ(x,t)=f_0 cosω t. 
We change variables to (τ,y)
and introduce the friction parameter σ:

    t = τ/ω,    x(t) = (f_0/ω^2)y(τ)  ;
    
    σ=ν f_0^2n/ω .


Then Eq. (<ref>) becomes

    y”(τ) + σcos^2n(τ) y'(τ) = cosτ  .

The general solution can be represented as


    y'(τ)=
    e^-σ S_n(τ)
    [ y'(0) - Y_n(0) ] +  Y_n(τ)
      ,
     
       Y_n (τ)  =
    e^-σ S_n(τ)/e^πασ+1∫_-π/2^π/2
    e^σ S_n(ζ+τ+π/2)sin(ζ + τ)dζ
    ,
    
       
    S_n(τ) = 
    ατ+
    ∑_m=1^n
    2nn+msin(2mτ)/2^2nm
    ,
      α=2^-2n2nn
    .


As one can see, any solution at τ→∞ 
tends to the limit cycle described by the function Y_n.


The amplitude of the limit cycle in terms of the derivative y',
A_n=max|Y_n(τ+π)-Y_n(τ)|, is

    A_n=
    1/cosh(πασ/2)∫_-π/2^π/2
    e^σ S_n(τ)cos(τ)dτ
    .

For large σ, it decreases  as a negative power of σ,
e.g., A_1 ≈ 1.88σ^-2/3.

Correspondingly, the oscillation amplitude decreases
when the trajectory drifts towards the driving force maximum,
as seen in Fig. <ref>.

In the case of n=1,
the function describing the limit cycle, Eq. (<ref>),
can be represented as a Fourier series
in terms of odd harmonics of the driving force


    Y_1(τ)=
    ∑_m=1^∞[
    
    C_me^i(2m-1)τ+C_m^*e^-i(2m-1)τ]
      ,


where the asterix denotes complex conjugation.
The sequence C_1, C_m, m≥2,
representing the frequency spectrum of Y_1,
is expressed in terms of modified Bessel functions
of the first kind, I_k:

    C_1=
    ∑_k=-∞^∞(-1)^k+1/4k+2-iσ[
    i I_k^2(σ/4) + I_k(σ/4)I_k+1(σ/4)
    ]
      ,
    
    
    C_m =
    i^m
    I_m-1/2-iσ/4(σ/4)/I_3/2-iσ/4(σ/4)[
    (2+4i/σ)C_1 +C_1^* - 2/σ]
    .


The amplitude of the first harmonic, C_1, as a function of σ can be approximated by

    C_1(σ)≈σ/2/4+σ^2-δ
    -
    i[
    σ/3-2/4+σ^2-δ
    +
    16+3σ^2/(4+σ^2-δ)^2]

with δ≈ 0.098.

The spectral density of the 
function Y_1(τ),
Eq. (<ref>),
is |2C_m|^2;
it is shown in
Fig. <ref>.

As one can see, 
the frequency spectrum contains high order harmonics according to Eqs. (<ref>-<ref>).






§ CONCLUSION


In conclusion, in contrast to known dynamical destabilization 
under the action of dissipation known as dissipation-induced instabilities 
(see review article <cit.> and references therein) we show that
a strong nonlinear friction can cause
a seemingly paradoxical stabilization of forced oscillations
near the maxima of the driving force spatial amplitude.

In particular, such a friction
occurs in the dynamics of charged particles
in ultra-strong electromagnetic fields
due to radiation reaction <cit.>.

The threshold for the described stabilization of forced oscillations
corresponds to the criterion of the importance of the radiation reaction effects
<cit.>.


In a standing electromagnetic waves,
which can be formed in multiple high power laser configurations <cit.>,
the stabilization due to a strong radiation reaction
is manifested in
an anomalous electron bunching near the electric field maxima <cit.>.

As shown in Refs. <cit.>,
electrons  can be captured for many laser periods due to
radiation friction impeding the ponderomotive force.

When radiation reaction dominates,
the electron motion in a standing wave 
evolves to limit cycles and strange attractors <cit.>.

In the case of a circularly polarized standing wave,
analytical expressions exist for the limit cycles
near the electric field maxima <cit.>.

A collision of multiple ultra-intense electromagnetic waves
creates structurally determinate patterns 
in the electron phase space <cit.>
due to a counterplay of the ponderomotive force
and the friction-induced tug.


Although in the present work we have been motivated 
by the intention to build up the theory of the radiative electron dynamics  
in the field of extremely high intensity lasers, we believe that the formulated above 
concept of dissipative stabilization of nonlinear dynamic systems will be useful for 
applications well beyond the framework of the 
laser-matter interaction physics <cit.>,
remembering a saying of William Thomson (Lord Kelvin) 
“I never satisfy myself until I can make a mechanical model of a thing”
<cit.>.






99

Kapitza1
P. L. Kapitza, 

Soviet Phys. JETP 21, 588 (1951).

Kapitza2
P. L. Kapitza,

Usp. Fiz. Nauk 44, 7 (1951).

Stephenson
A. Stephenson, 

Philos. Mag., Ser. 6, 15, 233 (1908).

Blekhman
I. I. Blekhman, Vibrational Mechanics
(World Scientific Publishing, 2000).

Chelomei V. N. Chelomei, Soviet Physics Doklady 28, 387 (1983). 


Butikov
E. I. Butikov,

J. Phys. A: Math. Theor. 44, 295202 (2011).

LLM
L. D. Landau and E. M. Lifshitz, Mechanics (Pergamon, New York, 1976), chap. 30.

Licht-Lieb
A. J. Lichtenberg and M. A. Lieberman,
Regular and chaotic dynamics

(Springer-Verlag, New York, 1992), p. 144.

Destabilization
R. Krechetnikov and J. E. Marsden, Rev. Mod. Phys. 79, 519 (2007).

Rad-Dom
A. Di Piazza, C. Müller, K. Z. Hatsagortsyan, C. H. Keitel,

Rev. Mod. Phys. 84, 1177 (2012).




Survey
S. V. Bulanov, T. Zh. Esirkepov, S. S. Bulanov, J. K. Koga, Z. Gong, et al.,



arXiv:1701.03349.

Multi-las
S. S. Bulanov, V. D. Mur, N. B. Narozhny, J. Nees, V. S. Popov,
Phys. Rev. Lett. 104, 220404, (2010).

Anom-bunch
A. Gonoskov, A. Bashinov, I. Gonoskov, C. Harvey, A. Ilderton, et al.,

Phys. Rev. Lett. 113, 014801 (2014).

Fedotov
A. M. Fedotov, N. V. Elkina, E. G. Gelfer, N. B. Narozhny, H. Ruhl,

Phys. Rev. A 90, 053847 (2014).

Attr
T. Zh. Esirkepov, S. S. Bulanov, J. K. Koga, M. Kando, K. Kondo, et al.,


Phys. Lett. A 379, 2044 (2015).

Jirka
M. Jirka, O. Klimo, S. V. Bulanov, T. Zh. Esirkepov, E. Gelfer, et al.,

Phys. Rev. E 93, 023207 (2016).

Kirk
J. G. Kirk,

Plasma Phys. Control. Fusion 58, 085005 (2016).

Vranic
M. Vranic, T. Grismayer, R. A. Fonseca, L. O. Silva,

Plasma Phys. Control. Fusion 59, 014040 (2017).

Gong
Z. Gong, R. H. Hu, Y. R. Shou, B. Qiao, C. E. Chen, et al.,


Phys. Rev. E 95, 013210 (2017).

MTB 
G. A. Mourou, T. Tajima, S. V. Bulanov, Rev. Mod. Phys. 78, 309 (2006).

Marklund
M. Marklund and P. Shukla,

Rev. Mod. Phys. 78, 591 (2006).

Kelvin
W. Thomson (Lord Kelvin),
 Notes of Lectures on Molecular Dynamics and the Wave Theory of Light.
Delivered at The Johns Hopkins University, Baltimore, by Sir William Thomson, Professor in the University of Glasgow. Stenographically Reported by A.S.Hathaway, Lately Fellow in Mathematics of The Johns Hopkins University. (Baltimore: Johns Hopkins, 1884) 270-271.







equationsubsection




§ APPENDIX



Here we present mathematical derivations for some formulae shown above.
Notations for variables are the same as in the main text.
The equation numbering is preserved for those equations which appear in the main text;
auxiliary formulae are numbered within sections.





 §.§ Equations (14) and (15)

In this section we derive Eqs. (<ref>) and (<ref>) from Eq. (<ref>).
We use the following formulae explicitly written or assumed in the main text:

    F = f_1 cos(ω t) + f_2 sin(ω t),
    
    
    K = ν[ f_1 cos(ω t) + f_2 sin(ω t) ]^2n,
    
    
    12ξ=(κ f_1-ω f_2) sinω t
    -(κ f_1+ω f_2) cosω t/ω(κ^2+ω^2)  ,
    
    
    13Ẍ+(κ+⟨ξ∂_X K⟩) Ẋ=
    ⟨ξ∂_X F ⟩ -
    ⟨ξξ̇∂_X K⟩  .


The time-averaged friction coefficient, defined in Eq. (8) in the main text, is

    κ=⟨ K⟩=
    ω/2π∫_0^2π/ων[ f_1 cos(ω t) + f_2 sin(ω t) ]^2n dt=
    
    1/2π∫_0^2πν[ f_1-i f_2/2e^i τ + f_1+i f_2/2e^-i τ]^2n dτ=
    
       =ν/2π∑_m=0^2n2nm(f_1+i f_2/2)^m(f_1-i f_2/2)^2n-m∫_0^2π
    e^2(n-m)iτ
    dτ=
    ν/2π2nn(f_1+i f_2/2)^n(f_1-i f_2/2)^n2π
    
    .

Here 2nm=(2n)!/[m!(2n-m)!] is the binomial coefficient.
The last integral (marked with red) is zero except the case m=n,
for which it equals 2π (therefore the sum contains the only nonzero term,
for the index m=n).
In this way we obtain Eq. (<ref>) from the main text:

    15κ=
    2^-2n2nnν
    (f_1^2+f_2^2)^n
    .


In Eq. (<ref>),
the term ⟨ξ∂_X K⟩ 
is zero because
it contains only odd harmonics,

    ⟨ξ∂_X K⟩ = 
    ω/2π∫_0^2π/ω
    ξ(t) ×
    2n ν[ f_1 cos(ω t) + f_2 sin(ω_t) ]^2n-1[ ∂_X f_1 cos(ω t) + ∂_X f_2 sin(ω t) ]
    dt
    =0.


The next averaged term is obtained by a simple integration

    ⟨ξ∂_X F ⟩ = 
    ω/2π∫_0^2π/ω(κ f_1-ω f_2) sinω t-(κ f_1+ω f_2) cosω t/ω(κ^2+ω^2)[ f_1 cos(ω t) + f_2 sin(ω_t) ]
    dt=
    
    
    =
    -f_1 ∂_X f_1 + f_2 ∂_X f_2/2(κ^2+ω^2)
    -κ(f_2 ∂_X f_1 - f_1 ∂_X f_2)/2ω(κ^2+ω^2)
    .


The last averaged terms is

    -⟨ξξ̇∂_X K⟩   = 
    -ω/2π∫_0^2π/ωξ(t)ξ̇(t)×
    2n ν[ f_1 cos(ω t) + f_2 sin(ω_t) ]^2n-1[ ∂_X f_1 cos(ω t) + ∂_X f_2 sin(ω t) ]
    dt=
    
       =
    -
    2nν/8ω(κ^2+ω^2)^2∑_m=0^2n-12n-1m(f_1+i f_2/2)^m(f_1-i f_2/2)^2n-1-m×
       ×{
    (κ - iω)^2(f_1-i f_2)^2
    	[
    	(∂_X f_2 -i∂_X f_1) δ_m,n
    	-i(∂_X f_1 -i∂_X f_2) δ_m,n+1]
    +
    .
    
       .
    +
    (κ + iω)^2(f_1+i f_2)^2
    ×[
    	(∂_X f_2 + i∂_X f_1) δ_m,n-1
    	+i(∂_X f_1 + i∂_X f_2) δ_m,n-2]
    }
    =
    
       =
    
    n κ/(n+1)ω(κ^2+ω^2)^2[
    2nκω(f_1 ∂_X f_1 + f_2 ∂_X f_2)
    +(κ^2-ω^2)(f_2 ∂_X f_1 - f_1 ∂_X f_2)
    ].

Here δ_m,n is the Kronecker delta;
it equals one for m=n and zero otherwise.
In the last line ν is expressed in terms of κ using Eq. (<ref>).

Combining the results  into Eq. (<ref>)
we obtain Eq. (<ref>) from the main text:


    14Ẍ+κẊ=
    -∂ _X (f_1^2+f_2^2)/4(κ^2+ω^2)
    +
    n^2κ^2 ∂ _X (f_1^2+f_2^2)/(n+1)(κ^2+ω^2)^2
    +
    κ (f_2 ∂_X f_1 - f_1 ∂_X f_2)/2ω(κ^2+ω^2)[
    2n/n+1(
    κ^2-ω^2/κ^2+ω^2)
    -1
    ]
    .






 §.§ Equations (20-22)


In this section we derive Eqs. (<ref>),
(<ref>) and (<ref>),
representing the solution of Eq. (<ref>):

    19
    y”(τ) + σcos^2n(τ) y'(τ) = cosτ  .

We seek the solution for the first derivative y' in the form

    y'(τ) = M(τ)exp(-σ∫_0^τcos^2n(η)dη)
    .

Substituting this ansatz into Eq. (<ref>),
we obtain 

    M'(τ) = exp(σ∫_0^τcos^2n(η)dη)cos(τ)
    .

Integrating this equation for M and substituting the result into Eq. (<ref>),
we obtain a general solution of Eq. (<ref>):

    y'(τ) = exp(-σ∫_0^τcos^2n(η)dη)
    [y'(0) +
    ∫_0^τexp(σ∫_0^ζcos^2n(η)dη)
    cos(ζ)
    dζ]
    .


The term in the exponent can be expanded as

    S_n(τ) 
       = ∫_0^τcos^2n(η)dη
    = ∫_0^τ(e^iη+e^-iη/2)^2ndη
    =∫_0^τ∑_m=0^2n1/2^2n2nm
    e^i (2n-m)ηe^-i mηdη
    =
    
       =1/2^2n∑_m=0^2n2nm∫_0^τ e^2i (n-m)ηdη
    = 1/2^2n2nnτ +
    1/2^2n∑_m=0
     m≠n^2n2nme^2i (n-m)τ - 1/2i(n-m)
    =
    
       =1/2^2n2nnτ +
    1/2^2n∑_m=0
    m→ n-m'^n-12nme^2i (n-m)τ - 1/2i(n-m)
    +
    1/2^2n∑_m=n+1
    m→ n+m'^2n2nme^2i (n-m)τ - 1/2i(n-m)
    =
    
       =1/2^2n2nnτ +
    1/2^2n∑_m'=1^n2nn-m'e^2i m' τ - 1/2i m'
    +
    1/2^2n∑_m'=1^n2nn+m'e^- 2i m'τ - 1/-2i m'
    =
    
       =1/2^2n2nnτ +
    1/2^2n∑_m=1^n2nn+me^2i m τ - e^-2i m τ/2i m
       =1/2^2n2nnτ +
    1/2^2n∑_m=1^n2nn+msin(2m τ)/m
    .


Here we rearranged the sums changing to new indices (marked with red)
and used the identity
2nn-m=2nn+m for n≥ m≥ 0.

In this way we obtain Eq. (<ref>) in the main text:

    S_n(τ) = 
    ατ+
    ∑_m=1^n
    2nn+msin(2mτ)/2^2nm
    ,
      α=2^-2n2nn
    .
    22


Below we use the following properties of the function S_n:

    S_n(-τ)=S_n(τ),   
    S_n(τ+π)=S_n(τ)+πα,   
    S_n(π/2)=πα/2.


The limit cycle is a periodic solution Y_n of Eq. (<ref>),
therefore for any τ

    Y_n(τ) = Y_n(τ+2π).

Using (<ref>) we rewrite this expression as

    Y_n(τ)
       =
    e^-σ S_n(τ+2π)[y'(0) +
    ∫_0^τ+2π
    e^σ S_n(ζ)cos(ζ)
    dζ]
    =
    
       =
    e^-2πασe^-σ S_n(τ)[
    y'(0) +
    ∫_0^τ
    e^σ S_n(ζ)cos(ζ)
    dζ
    +
    ∫_τ^τ+2π
    e^σ S_n(ζ)cos(ζ)
    dζ]
    =
    
       =
    e^-2πασ[
    Y_n(τ) +
    e^-σ S_n(τ)∫_τ^τ+2π
    e^σ S_n(ζ)cos(ζ)
    dζ]
    ,

which gives the formula for Y_n(τ)

    Y_n(τ) 
       =
    e^-σ S_n(τ)/e^2πασ-1∫_τ^τ+2π
    e^σ S_n(ζ)cos(ζ)
    dζ
    =
    e^-σ S_n(τ)/e^2πασ-1∫_0^2π
    e^σ S_n(ζ+τ)cos(ζ+τ)
    dζ=
    
       =e^-σ S_n(τ)/e^2πασ-1[
    ∫_0^π
    e^σ S_n(ζ+τ)cos(ζ+τ)
    dζ
    +∫_π^2π
    e^σ S_n(ζ+τ)cos(ζ+τ)
    dζ]=
    
       =
    e^-σ S_n(τ)/e^2πασ-1[
    ∫_0^π
    e^σ S_n(ζ+τ)cos(ζ+τ)
    dζ
    -∫_0^π
    e^σ S_n(ζ+τ)+πασcos(ζ+τ)
    dζ]=
    
       =
    e^-σ S_n(τ)/e^2πασ-1
    (1-e^πασ)
    ∫_0^π
    e^σ S_n(ζ+τ)cos(ζ+τ)
    dζ
    =
    e^-σ S_n(τ)/e^πασ+1∫_-π/2^π/2
    e^σ S_n(ζ+τ+π/2)sin(ζ+τ)
    dζ
    .

In this way we obtain Eq. (<ref>) in the main text:

    Y_n (τ)  =
    e^-σ S_n(τ)/e^πασ+1∫_-π/2^π/2
    e^σ S_n(η+τ+π/2)sin(η + τ)dη
    .
    
    21


Accepting Y_n as a particular solution of Eq. (<ref>)
and rewriting Eq. (<ref>),
we obtain another form of the general solution of Eq. (<ref>),
namely Eq. (<ref>):

    y'(τ)=
    e^-σ S_n(τ)
    [ y'(0) - Y_n(0) ] +  Y_n(τ)
      .
    
    20





 §.§ Equation (23)

In this section we derive Eq. (<ref>).
The maximum difference |Y_n(τ+π)-Y_n(τ)|
is reached for τ=-π/2,
therefore

    A_n   =Y_n(π/2)-Y_n(-π/2)=
    
       =
    e^-σ S_n(π/2)/e^πασ+1∫_-π/2^π/2
    e^σ S_n(η+π)sin(η + π/2)dη
    -
    e^-σ S_n(-π/2)/e^πασ+1∫_-π/2^π/2
    e^σ S_n(η)sin(η - π/2)dη
    =
    
       =
    e^-πασ/2/e^πασ+1∫_-π/2^π/2
    e^σ S_n(η)+πασcos(η)dη
    +
    e^πασ/2/e^πασ+1∫_-π/2^π/2
    e^σ S_n(η)cos(η)dη
    =
    
       =2e^πασ/2/e^πασ+1∫_-π/2^π/2
    e^σ S_n(η)cos(η)dη
    =1/cosh(πασ/2)∫_-π/2^π/2
    e^σ S_n(η)cos(η)dη
    .

Thus we obtain Eq. (<ref>) in the main text:

    23A_n=
    1/cosh(πασ/2)∫_-π/2^π/2
    e^σ S_n(τ)cos(τ)dτ
    .






 §.§ Equations (25) and (26)

In this section we derive Eqs. (<ref>) and (<ref>)
for the case of n=1.
The function Y_1 is a periodic solution
of the equation

    Y'_1(τ) + σcos^2(τ) Y_1(τ) = cosτ  .

It is represented as a Fourier series, which
obviously should contain only odd hamonics of the driver:


    Y_1(τ)=
    ∑_m=1^∞[
    C_me^i(2m-1)τ+C_m^*e^-i(2m-1)τ]
      ,
    
    24


where symbol “*” denotes complex conjugation.

Substituting Eq. (<ref>)
into Eq. (<ref>)
we obtain

    ∑_m=1^∞[
    i(2m-1) C_me^i(2m-1)τ - i(2m-1) C_m^*e^-i(2m-1)τ]
    +
    
       +
    σ(e^2iτ/4+1/2+e^-2iτ/4)
    ∑_m=1^∞[
    C_me^i(2m-1)τ + C_m^*e^-i(2m-1)τ]
    =e^iτ/2+e^-iτ/2
    .


Rearranging the sums and collecting the Fourier coefficients of the terms corresponding to the same harmonics, we easily find the following recursive relations:

    C_2=-[
    (2+4i/σ)C_1 +C_1^* - 2/σ] , 
    
    i^-m C_m = σ/4/2(m-1/2-iσ/4)(
    i^-(m-1) C_m-1 - i^-(m+1) C_m+1).

The last expression is the same as a recurrence identity for the 
modified Bessel function of the first kind <cit.>,

    I_μ(z) = z/2μ(I_μ-1(z)-I_μ+1(z))

for z=σ/4 and μ=m-1/2-iσ/4,
and I_μ(z) ∝ i^-m C_m.
The general solution of Eq. (<ref>) is
a linear combination of the modified Bessel functions of
the first and second kinds, where 
one should cancel out an unbounded term:

    C_m = U i^m I_m-1/2-iσ/4(σ/4).

The coefficient of proportionality, U, is determined using the relation Eq. (<ref>):

    C_2 = -U I_3/2-iσ/4(σ/4)
    =-[
    (2+4i/σ)C_1 +C_1^* - 2/σ]
    .


As a result we obtain Eq. (<ref>) of the main text:

    C_m =
    i^m
    I_m-1/2-iσ/4(σ/4)/I_3/2-iσ/4(σ/4)[
    (2+4i/σ)C_1 +C_1^* - 2/σ]
    .
    
    26


The formula Eq. (<ref>) for the coefficient C_1
can be obtained using series representation for the function Y_1.
For n=1, we have α=1/2, S_1(τ)=τ/2+sin(2τ)/4,
and Eq. (<ref>) becomes

    Y_1 (τ)  =
    e^πσ/4/e^πσ/2+1exp(-σ/4sin(2τ))
    ∫_-π/2^π/2exp(σ/2η)
    exp(-σ/4sin(2η+2τ))
    sin(η + τ)dη
    .

Using the generating function for the modified Bessel function of the first kind <cit.>,

    exp[ z/2(p+1/p) ]
    =
    ∑_k=-∞^∞ I_k(z) p^k
    ,

we change exponent terms involving sin function into series
in the following way

    exp(-σ/4sin(2ζ))
    =
    ∑_k=-∞^∞ I_k(σ/4) i^k e^2ikζ
    .

Then the right-hand side of Eq. (<ref>) is transformed into a double series

    Y_1 (τ)  =
    e^πσ/4/e^πσ/2+1
    ∑_m=-∞^∞∑_k=-∞^∞ 
    i^m+k
    I_m(σ/4) 
    I_k(σ/4)
    e^2imτ
    ∫_-π/2^π/2
    e^ση/2 + 2ik(η+τ)
    sin(η + τ)dη
    .

Integrating term-by-term and rearranging the sums
by collecting terms involving the same harmonics,
it is not difficult to obtain the following formulae

    Y_1(τ) = ∑_m=-∞^∞(
    C_m e^i (2m-1) τ + C_1-m e^-i (2m-1) τ),
    
    
    C_m = i^m
    ∑_k=-∞^∞(-1)^k+1/4k+2-iσ[
    I_k(σ/4) - i I_k+1(σ/4)
    ] I_k+1-m(σ/4).

Using Eq. (<ref>) and
the symmetry property I_-k(z)=I_k(z) <cit.>,
it is easy to show that C_1-m=C^*_m.
For m=1 we obtain Eq. (<ref>) from the main text:

    C_1=
    ∑_k=-∞^∞(-1)^k+1/4k+2-iσ[
    i I_k^2(σ/4) + I_k(σ/4)I_k+1(σ/4)
    ].
    
    25


The representation Eq. (<ref>) is equivalent to
Eq. (<ref>) provided that C_1 is defined by Eq. (<ref>).

Fig. <ref> shows a comparison of the numerical solution of Eq. (<ref>)
for σ=20 with the analytical solution defined
by Eqs. (<ref>), (<ref>), and (<ref>),
where the sum in Eq. (<ref>) is cut at m=6
(which corresponds to the 2m-1=11^ th harmonic)
and the sum in Eq. (<ref>) is cut at k=± 10
(i.e. only terms with the index k from -10 to 10 are taken into account).



99


math
M. Abramowitz and I. A. Stegun,
Handbook of Mathematical Functions, With Formulas, Graphs, and Mathematical Tables
(Dover Publications, 1974) chapter 9, section 6.




