To appear in CHAOS, Volume 17, Issue 2, June 2007
On the existence of chaotic circumferential waves in spinning disks
Arzhang Angoshtari∗ and Mir Abbas Jalali†
Center of Excellence in Design, Robotics and Automation
Department of Mechanical Engineering, Sharif University of Technology
P.O.Box: 11365-9567, Azadi Avenue, Tehran, Iran
(Dated: August 4, 2021)
We use a third-order perturbation theory and Melnikov’s method to prove the existence of chaos
in spinning circular disks subject to a lateral point load. We show that the emergence of transverse
homoclinic and heteroclinic points respectively lead to a random reversal in the traveling direction
of circumferential waves and a random phase shift of magnitude π for both forward and backward
wave components. These long-term phenomena occur in imperfect low-speed disks sufficiently far
from fundamental resonances.
Keywords: Chaos, canonical perturbation theory, Melnikov’s method, spinning disks, traveling wave reversal
Transversal vibration modes of hard disk drives
(HDDs) are excited by the lateral aerodynamic
force of the magnetic head. Previous works [1, 2]
revealed that chaotic orbits are inevitable ingre-
dients of phase space flows when the lateral force
is large, or the disk is rotating near the critical
resonant speed. For low-speed disks, however, an
adiabatic invariant (a first integral) was found [2]
using a first-order averaging based on canonical
Lie transforms. According to the first-order the-
ory, regular vibrating modes of imperfect, low-
speed disks are independent of the angular veloc-
ity of the disk, Ω0. In such a circumstance, the
speed of circumferential waves is the natural fre-
quency of the lateral mode, ω, derived from linear
vibration analysis. HDDs are usually operated
with angular velocities smaller than ω (safely be-
low resonance). Moreover, the magnitude of the
lateral force F is very small. Given these condi-
tions, we show that it is impossible to continue
the Lie perturbation scheme [2] up to terms of
arbitrary order and remove the time variable t
from the Hamiltonian. In fact, due to the spe-
cial forms of nonlinearities in the dynamical equa-
tions of spinning disks, one can not remove t from
third-order terms. Subsequent application of a
second-order Melnikov theory reveals that trans-
verse homoclinic and heteroclinic points do exist
for all F,Ω0 6= 0. This implies chaos, or equiva-
lently, non-integrability of governing equations.
I. INTRODUCTION
Dynamics of continuum media, like fluids, rods, plates
and shells is usually formulated as a system of partial
differential equations (PDEs) for physical quantities in
∗Electronic address: arzhang˙a@mehr.sharif.edu
†Electronic address: mjalali@sharif.edu;
URL: http://sharif.edu/~mjalali
terms of the spatial coordinates x and the time t as
L(u) = 0. (1)
Here L is a nonlinear operator and u(x, t) is the vector
of dependent variables. When the boundary conditions
are somehow simple, approximate variational methods
based on modal decomposition and Galerkin’s projection
[3] can be used to reduce the order of governing equa-
tions. These methods begin with solving an auxiliary
eigenvalue problem, which is usually the variational field
equation δL(u,Λ) = 0, and build some complete basis
set Uk(x,Λk) for expanding u(x, t) in the spatial do-
main. Here Λk is an eigenvalue, which is characterized
by the vectorial index k. EachUk is called an eigenmode
or a shape function. Such a basis set should preferably
satisfy boundary conditions and be orthogonal.
Once a complete basis set is constructed, one may sup-
pose a solution of the form
u(x, t) =
Vk(t) ·Uk(x,Λk). (2)
Substituting from (2) into (1) and taking the inner prod-
L(u) ·Uk′(x)dx = 0, (3)
leaves us with a system of nonlinear ordinary differen-
tial equations (ODEs) for the amplitude functions Vk(t).
The evolution of the reduced ODEs shows the interaction
of different modes and their influence on the development
of spatiotemporal patterns.
In a series of papers, Raman and Mote [1, 4] used the
modal decomposition method to investigate transversal
oscillations of spinning disks whose deformation field is
described in terms of the displacement vector (u, v, w)
with u and v being the in-plane components. The most
important application of the spinning disk problem is in
design, fabrication and control of HDDs. The governing
PDEs for the evolution of displacement components were
first derived by Nowinski [5] and reformulated more re-
cently by Baddour and Zu [6]. Let us define (R, φ) as
http://arxiv.org/abs/0704.1406v1
mailto:arzhang_a@mehr.sharif.edu
mailto:mjalali@sharif.edu
http://sharif.edu/~mjalali
the usual polar coordinates. For a rotating disk with the
angular velocity Ω0, Nowinski’s theory assumes that the
in-plane inertias (u, v)Ω20, 2Ω0(u,t, v,t) and (u,tt, v,tt) are
ignorable againstRΩ20. This is a rough approximation for
high-speed disks and one needs to use the complete set
of equations as Baddour and Zu [6] suggest. Nowinski’s
theory, however, has its own advantages (like the exis-
tence of a stress function) that facilitate the study of the
most important transversal modes. In low-speed disks, or
disks with high flexural rigidity, one has Ω0 ≪ ω. Hence,
it is legitimate to apply Nowinski’s governing equations
in such systems.
In this paper we analytically prove the existence of
chaos, and therefore, non-integrability of the reduced
ODEs that govern the double-mode oscillations of im-
perfect spinning disks. We investigate low-speed disks
subject to a lateral point force exerted by the magnetic
head. The lateral force in HDDs is very small and its
origin is the aerodynamic force due to air flow in the gap
between the disk and the head. We show that chaotic
circumferential waves dominate some zones of the phase
space over the time scale t ∼ O(ǫ−3) with ǫ being a small
perturbation parameter. This indicates very slow evolu-
tion of random patterns, and the practical difficulties of
their identification.
The paper is organized as follows. In §II, we present
the Hamiltonian function in terms of Deprit’s [7] Lis-
sajous variables. In §III, we use a canonical perturba-
tion theory [8, 9] to eliminate the fast anomaly l from
the Hamiltonian. The action associated with l then be-
comes an adiabatic invariant. Transversal intersections
of destroyed invariant manifolds, and therefore, non-
integrability of the normalized equations, is proved by
a second-order Melnikov method in §IV. We present a
complete classification of circumferential waves in §V and
end up the paper with concluding remarks in §VI.
II. PROBLEM FORMULATION
Let us assume Um(R) as an orthogonal basis set that
represents the disk deformation in the radial direction.
The index m stands for the number of radial nodes
that Um(R) has. According to Raman and Mote’s [1]
treatment of imperfect disks, the following choice of the
transversal displacement field
w(R, φ, t) = Um(R)[x(t) cosnφ+ y(t) sinnφ], (4)
reduces Nowinski’s governing equations to a system of
ODEs for the amplitude functions x(t) and y(t) as
ẍ+ λ2x+ ǫγ
x2 + y2
x = ǫF cos(nΩ0t), (5a)
ÿ + ω2y + ǫγ
x2 + y2
y = ǫF sin(nΩ0t), (5b)
where λ, ω and γ are constant parameters that depend
on the geometry and material of the disk. ǫ is a small
perturbation parameter, F is the weighted integral of the
lateral point force, and Ω0 is the angular velocity of the
disk.
We suppose small deviations from perfect disks and
write the constant parameter of (5) as λ2/ω2 = 1 + ǫη.
We also define nΩ0 = ǫΩ with O(Ω) ∼ O(ω). Denoting
(px, py) as the momenta conjugate to (x, y), it can be
verified that equations (5) are derivable from the Hamil-
tonian function
p2x + p
x2 + y2
+ ΩP+
x2 + y2
)2 − F (x cos p+ y sin p)
. (6)
We have introduced the action P and its conjugate an-
gle p = ǫΩt to make our equations autonomous, which
is a preferred form for the application of canonical per-
turbation theories. The extended phase space has now
dimension six. Dynamics generated by (6) is better un-
derstood after carrying out a canonical transformation
(x, y, px, py) → (l, g, L,G) to the space of Lissajous vari-
ables [7] so that
x = s cos(g + l)− d cos(g − l), (7a)
y = s sin(g + l)− d sin(g − l), (7b)
px = −ω [s sin(g + l) + d sin(g − l)] , (7c)
py = ω [s cos(g + l) + d cos(g − l)] , (7d)
, d =
, L ≥ 0, |G| ≤ L.
In the space of Lissajous variables, the Hamiltonian de-
fined in (6) becomes
H = H0(L) + ǫH1(l, g, p, L,G, P ), (8)
H0 = ωL,
H1 = ΩP − F [s cos(g + l − p)− d cos(g − l − p)]
(s2 + d2)− 2sd cos(2l)
(s2 + d2) + s2 cos(2g + 2l)− 2sd cos(2g)
+ d2 cos(2g − 2l)− 2sd cos(2l)
From (8) we conclude that l is the fast angle, and g and
p are the slow ones. Therefore the long-term behavior
of the flows generated by (8) can be analyzed by averag-
ing H over l. After removing l, its corresponding action
L will be a constant of motion for the flows generated
by the averaged Hamiltonian 〈H〉l, and the phase space
dimension reduces from 6 to 4.
III. CANONICAL THIRD-ORDER AVERAGING
In order to average H over l, we use the normal-
ization procedure of Deprit and Elipe [9]. Denoting
X ≡ (l, g, p) and Y ≡ (L,G, P ), we define a Lie transfor-
mation (l, g, p, L,G, P ) → (l̄, ḡ, p̄, L̄, Ḡ, P̄ ) as
X = EW (X̄), Y = EW (Ȳ ), (9)
so that the Hamiltonian function in terms of the new
variables, K ≡ 〈H〉l, does not depend on l̄. EW is the Lie
transform generated by the function W and it is defined
EW (Z̄) = Z̄ + (Z̄;W ) +
((Z̄;W );W )
(((Z̄;W );W );W ) + · · · . (10)
In this equation, (f1; f2) denotes the Poisson bracket of
f1 and f2 over the (l̄, ḡ, p̄, L̄, Ḡ, P̄ )-space. We expand the
generating function W as
W = ǫW1 +
W2 + · · · , (11)
and specify the averaged, target Hamiltonian K =
K(ḡ, p̄, L̄, Ḡ, P̄ ) as the series [9]
K = K0 + ǫK1 +
K3 + · · · , (12)
K0 = ωL̄, (13a)
H1dl̄, (13b)
[2(H1;W1) + ((H0;W1);W1)]dl̄, (13c)
3(H1;W2) + 3((H1;W1);W1)
+ 2((H0;W2);W1) + ((H0;W1);W2)
+ (((H0;W1);W1);W1)
dl̄. (13d)
W1 and W2 are determined through solving the following
differential equations
= H1(l̄, ḡ, p̄, L̄, Ḡ, P̄ )−K1(ḡ, p̄, L̄, Ḡ, P̄ ),(14a)
= 2(H1;W1) + ((H0;W1);W1)−K2. (14b)
By substituting from (14) into (13) and evaluating the
integrals, one finds the explicit form of the new Hamilto-
nian K, which has been given in Appendix A up to the
third-order terms.
Once l̄ is removed from the Hamiltonian, L̄ becomes
an integral of motion. The slow dynamics of the system
is thus governed by the flows in the (ḡ, Ḡ)-space. We
introduce the slow time τ = p̄/Ω, ignore the fourth-order
terms in ǫ, and obtain the following differential equations
for the dynamics of (ḡ, Ḡ)
= f1(ḡ, Ḡ) + ǫh1(ḡ, Ḡ, ǫ, τ), (15a)
= −∂K
= f2(ḡ, Ḡ) + ǫh2(ḡ, Ḡ, ǫ, τ), (15b)
where
f1 = d6 + d1 cos(2ḡ), f2 = e1 sin(2ḡ),
h1 = d7 + d2 cos(2ḡ)
+ ǫ[d8 + d3 cos(2ḡ) + d4 cos(4ḡ) + d5 cos(2ḡ − 2Ωτ)],
h2 = e2 sin(2ḡ)
+ ǫ[e3 sin(2ḡ) + e4 sin(4ḡ) + e5 sin(2ḡ − 2Ωτ)]. (16)
In these equations, di (i = 1, · · · , 8) and ej (j = 1, · · · , 5)
are functions of L̄ and Ḡ (see Appendix A). It is re-
marked that the action P̄ appears only in K1 via the
term ΩP̄ . It then disappears in the normalized equa-
tions (15) after taking the partial derivatives of K with
respect to ḡ and Ḡ. The partial derivative of K with
respect to P̄ determines the evolution of p̄, which is in
accordance with the simple linear law p̄(τ) = Ωτ + p̄(0).
The dynamics of P̄ itself is governed by
= −∂K
= − 1
K(ḡ, Ḡ, τ). (17)
One may integrate (17) to obtain P̄ (τ) once equations
(15) are solved. The behavior of P̄ is thus inherited from
ḡ(τ) and Ḡ(τ).
IV. THE MELNIKOV FUNCTION
There are few analytical methods in the literature for
the detection of chaos in perturbed Hamiltonian systems
[10, 11]. Melnikov’s [10] method is the most powerful
technique when the governing equations take the form
= f(x) + ǫh(x, ǫ, τ), x ∈ R2, (18)
so that the unperturbed system dx/dτ = f(x) is in-
tegrable and possesses a homoclinic (heteroclinic) orbit
qh(τ) to a hyperbolic saddle point, and h(x, ǫ, τ) is T -
periodic in τ . The occurrence of chaos is examined by
the Melnikov function
M(τ0, ǫ) = ǫM1(τ0) + ǫ
2M2(τ0) + . . . , (19)
where Mk(τ0) denotes the kth-order Melnikov function.
Assume that Mi(τ0) is the first nonzero term, i.e.,
Mk(τ0) ≡ 0 for 1 ≤ k ≤ i − 1. If Mi(τ0) has simple
zeros, then, for sufficiently small ǫ, the system (18) has
transverse homoclinic (heteroclinic) orbits, which imply
chaos due to the Smale-Birkhoff homoclinic theorem [12].
The first-order term in (19) is determined by the classical
formula
M1(τ0) =
f (qh (τ)) ∧ h (qh (τ) , 0, τ + τ0) dτ,(20)
where the wedge operator ∧ is defined as f ∧ h =
f1h2 − f2h1. Although the Hamiltonian equations (15)
have a suitable form for the application of Melnikov’s
method, they are autonomous up to the first-order terms
in ǫ. Consequently, M1(τ0) vanishes identically for all
τ0 ∈ [0, T ]. We thus need to investigate the second-order
Melnikov function. For doing so, we begin with solving
the unperturbed system
= f1(ḡ, Ḡ), (21a)
= −∂K1
= f2(ḡ, Ḡ), (21b)
along homoclinic (heteroclinic) orbits. Jalali and An-
goshtari [2] showed that for L̄ > ηω3/γ, equations (21)
have hyperbolic stationary points at S0 ≡ (ḡ0, Ḡ0) =
(−π, 0), S2 ≡ (ḡ2, Ḡ2) = (0, 0), and S4 ≡ (ḡ4, Ḡ4) =
(π, 0). The implicit equation of the invariant manifolds
that terminate at the saddle points are
cos[2ḡh(τ)] =
L̄− γ
Ḡ2h(τ)
L̄2 − Ḡ2h(τ)
]−1/2
For γL̄ ≥ 2ηω3, equation (22) represents a heteroclinic
orbit which connects S0 to S2. For ηω
3 < γL̄ < 2ηω3,
the heteroclinic orbit disappears and it is replaced by a
homoclinic orbit (see Figure 1). To compute the explicit
form of the homoclinic (or heteroclinic) orbit of (21), we
use (21b) and (22), and obtain
∫ Ḡh(τ)
Ḡh(0)
1− αḠ2
βτ, (23)
16ηβω4
, β =
γL̄− ηω3
where the lower integration limit is Ḡh(0) = 1/
α. After
taking the integral (23), we arrive at
Ḡh(τ) =
sech(
β τ)√
, (24)
for the Ḡ ≥ 0 branch of the homoclinic (heteroclinic)
orbit. Having Ḡh(τ), it is straightforward to calculate
cos[2ḡh(τ)] and sin[2ḡh(τ)], and determine the explicit
form of qh(τ) =
ḡh(τ), Ḡh(τ)
For constructing M2(τ0), we use Françoise’s [13, 14]
algorithm that has been devised for dynamical systems
with polynomial nonlinearities. To express the averaged
Hamiltonian K in terms of polynomial functions of some
new dependent variables, we utilize Hopf’s variables
L̄2 − Ḡ2 cos(2ḡ), (25a)
L̄2 − Ḡ2 sin(2ḡ), (25b)
and obtain the following differential 1-form for the evo-
lution of the averaged system
dQ1 +
+ ǫ [(s1 + ǫs2) dQ2 − (z1 + ǫz2) dQ1] = 0. (26)
Here, the first-order Hamiltonian is
K1(Q1, Q2) =
(Q21 +Q
Q1 + C, (27)
s1 = m1Q2, (28a)
z1 = n1Q1 + n6, (28b)
s2 = m2(Q
2)Q2 +m3Q1Q2 +m4Q2
+ m5 sin(2Ωτ), (28c)
z2 = n2(Q
2)Q1 + n3(3Q
2) + n4Q1
+ n5 cos(2Ωτ) + n7. (28d)
The constant coefficients C, mi (i = 1, · · · , 5), and nj
(j = 1, · · · , 7) have been given in Appendix B. A prereq-
uisite for the application of Françoise’s [13] algorithm is
that for all polynomial 1-forms D that satisfy the condi-
D ≡ 0, (29)
there must exist polynomials A(Q1, Q2) and r(Q1, Q2)
such that D = dA+ rdK1. We call this the condition (∗)
and prove in Appendix C that K1 satisfies the condition
Françoise’s algorithm states that if M1(τ0) = · · · =
Mk−1(τ0) ≡ 0 for some integer k ≥ 2, it follows that
Mk(τ0) =
Dk, (30a)
D1 = δ1, Dm = δm +
i+j=m
riδj , (30b)
δj = zjdQ1 − sjdQ2, (30c)
for 2 ≤ m ≤ k. The functions ri are then determined
successively from the formulas Di = dAi + ridK1 for
i = 1, · · · , k − 1. We have already found that
M1(τ0) =
δ1 = 0, (31)
δ1 = (n1Q1 + n6)dQ1 −m1Q2dQ2. (32)
From (30) and (C2) it can be shown that
M2(τ0) =
δ2. (33)
Substituting (30c) and (28) into (33), and carrying out
the integration along qh(τ), result in
M2(τ0) =
3γF 2I
sin(2Ωτ0), (34)
with I being a constant (see Appendix D). Equation (34)
shows that τn = nπ/(2Ω) (n = 1, 2, · · · ) are simple zeros
of M2(τ0) so that
M2(τn) = 0,
∂M2(τ0)
τ0=τn
6= 0. (35)
FIG. 1: Possible topologies of the phase space flows of the averaged system for F = 0. In all panels the horizontal axis indicates
ḡ/π and the vertical axis indicates Ḡ/L̄. In the first topology (left panel), all stationary points are of center type. In the second
topology (middle panel) two new centers, with non-zero Ḡ-coordinates, have emerged for ḡ = ±nπ (n = 0, 1) and symmetrical
homoclinic loops (thick lines) connect saddle points to themselves. In the third topology (right panel) the off-axis centers (and
their surrounding tori) are still present but the separatrix curves (thick lines) are of heteroclinic type.
FIG. 2: Phase space structure of the normalized equations (15) with L̄ = 2, ω = η = F = 1, and Ω = ǫ = 0.1. Left panel:
γ = 0.8. Right panel: γ = 2. In both panels the horizontal and vertical axes indicate ḡ/π and Ḡ/L̄, respectively.
Thus, we conclude that the global stable and unstable
manifolds of the saddle point Sτ0n , W
s(Sτ0n ) andW
u(Sτ0n ),
always intersect transversely. Transversal intersections
cause a sensitive dependence on initial conditions due to
the Smale-Birkhoff homoclinic theorem. This is a route
to chaos. On the other hand this means that the reduced
equations (15) are non-integrable for F 6= 0.
V. CLASSIFICATION OF CIRCUMFERENTIAL
WAVES
For F = 0, h does not depend on τ and the normalized
equations (15) are integrable. In such a circumstance, the
phase space structure can take three general topologies
(depending on the values of the system parameters and
L̄) as shown in Figure 1. In the first topology all sta-
tionary points with the coordinates (ḡs, Ḡs) calculated
f(ḡs, Ḡs) + ǫh(ḡs, Ḡs, ǫ) = 0, (36)
are centers and they lie on the Ḡ = 0 axis with ḡs =
−π + nπ/2 (n = 0, · · · , 4). In the second and third
topologies, two off-axis centers (with Ḡs 6= 0) come to
existence for ḡs = ±nπ (n = 0, 1) and the on-axis sta-
tionary points with the same ḡs = ±nπ become saddles.
In the second topology, each saddle point is connected to
itself by a homoclinic orbit, and in the third topology, a
heteroclinic orbit connects two neighboring saddle points.
The system with heteroclinic orbits allows for rotational
ḡ(τ) while in the system with homoclinic orbits ḡ(τ) is
always librating. Beware that this classification of phase
space flows is valid as long as ǫ is sufficiently small.
For F = 0, the phase space flows of (15) are
structurally stable (with no unbounded branches) and
the whole (ḡ, Ḡ)-space is occupied by periodic orbits
of period T (K). At the stationary points, one has
T (K(ḡs, Ḡs, L̄)) = 0. Given the invariance of L̄, and
the periodic solutions ḡ(τ) = ḡ(τ + T (K)) and Ḡ(τ) =
Ḡ(τ + T (K)), the anomaly l̄ is determined through solv-
= ω + ǫ
, (37)
which results in l̄ = ωt + ǫR(τ) with R(τ) = R(τ +
T (K)). According to (9), the functions g(t), G(t) and
L(t) are also periodic in t and we conclude that l(t) =
ωt + RW (t) with RW (t) = EW (l̄) − ωt being a small-
amplitude periodic function of t. The explicit from of
the circumferential wave will then become
w(R, φ, t)
Um(R)
L(t)+G(t)
cos [nφ−ωt−RW (t)−g(t)]
L(t)−G(t)
cos [nφ+ωt+RW (t)−g(t)] , (38)
which is composed of a forward and a backward traveling
wave. Due to the periodic nature of L(t) and G(t), when
the amplitude of the forward traveling wave is maximum,
that of the backward wave is minimum and vice versa.
As our results of §IV shows, the regular nature of trav-
eling waves is destroyed for F 6= 0 and a chaotic layer oc-
curs through the destruction of the homoclinic and hete-
roclinic orbits of (21). This happens over the time scale
τ ∼ O(ǫ−2) or t ∼ O(ǫ−3) (because p̄ is present only
in K3). Figure 2 shows Poincaré maps of the system
(15) for F 6= 0. The sampling time step in generating
the Poincaré maps has been 2π/Ω. It is seen that most
tori around elliptic fixed points are preserved. They cor-
respond to regular periodic and quasi-periodic solutions
of the normalized system. For chaotic flows, the func-
tions ḡ(τ) and Ḡ(τ) randomly change within the invari-
ant measure of the chaotic set. Consequently, the orig-
inal Lissajous variables g(t), L(t), G(t) and also RW (t)
become chaotic too.
For Ḡ(τ) > 0 and Ḡ(τ) < 0 the forward and the back-
ward traveling waves are the dominant components of
the circumferential wave, respectively. When the chaotic
layer emerges from the destroyed homoclinic orbits (left
panel in Figure 2), the sign of Ḡ(τ) is randomly switched
along a chaotic trajectory. This means a random trans-
fer of kinetic/potential energy between the forward and
backward traveling wave components. For chaotic trajec-
tories of this kind the angle ḡ(τ) randomly fluctuates near
ḡ ≈ ±nπ (n = 0, 1) with an almost zero average. The
evolution of circumferential waves is quite different when
the chaotic layer emerges due to destroyed heteroclinic
orbits (right panel in Figure 2). In this case chaos means
a random change between the librational and rotational
states of ḡ(τ). Such a change induces an unpredictable
phase shift of magnitude π for both forward and back-
ward traveling wave components. We note that Ḡ(τ) can
flip sign on a chaotic trajectory only when ḡ(τ) is in its
librational state.
VI. CONCLUDING REMARKS
Resonance overlapping [11, 16] is the main cause for
chaotic behavior in spinning disks with near-resonant an-
gular velocities [1, 2]. The chaos predicted in this paper,
however, happens far from fundamental resonances. Op-
tical and HHDs are usually operated below critical res-
onant speeds and the lateral force F due to magnetic
head is very small. We showed that whatever the magni-
tude of F may be, a chaotic layer fills some parts of the
phase space because the Melnikov function of the nor-
malized equations has always simple zeros. Dynamics of
rotating disks is regular only if F vanishes, which is an
unrealistic assumption for disk drives. In low-speed disks
with small F , diffusion of chaotic orbits (within their in-
variant measure) takes a long time of t ∼ O(ǫ−3). The
slow development of chaotic circumferential waves makes
them undetectable in short time scales at which most
controllers work. The Melnikov function (34) depends
not only on F , but also on the parameter η through the
constant I. The parameter η is a contribution of imper-
fections, which are likely because of limited fabrication
precision in micro/nano scales. For a perfect disk with
η = 0, the off-axis elliptic stationary points of (21), and
consequently, homoclinic and heteroclinic orbits disap-
pear. In such a condition the Melnikov function is indef-
inite, but the system admits an exact first integral and
the dynamics is governed by the Hamiltonian function
given in equation (11) of Jalali and Angoshtari [2].
One of the most important achievements of this work
was to unveil the fact that it is premature to truncate the
series of canonical perturbation theories before recording
the role of all participating variables. In systems with
non-autonomous governing ODEs (non-conservative sys-
tems), one must be cautious while removing a fast angle
through averaging schemes. The removal of the fast an-
gle may also wipe out time-dependent terms, up to some
finite orders of ǫ, and hide some essential information of
the underlying dynamical process. Strange irregular so-
lutions can indeed occur at any order and influence the
long term response of dynamical systems as we observed
for the spinning disk problem by keeping the third-order
terms.
Acknowledgments
We are indebted to the anonymous referee, who discov-
ered an error in the early version of the paper and led us
to investigate the second-order Melnikov function. MAJ
thanks the Research Vice-Presidency at Sharif University
of Technology for partial support.
APPENDIX A: THE NORMALIZED
HAMILTONIAN
By evaluating the integrals in (13), we obtain the first,
second and third order terms of the normalized Hamilto-
nian as
−γḠ2
L̄(3γL̄+ 2ηω3)
L̄2 − Ḡ2 cos(2ḡ), (A1)
K2 = −
2L̄(−9Ḡ2 + 17L̄2)γ2 + 64F 2ω3
+ 8(3L̄2 − Ḡ2)γηω3 + 8L̄η2ω6
+ 4ω3[−6L̄γη − 2η2ω3]
L̄2 − Ḡ2 cos(2ḡ)
, (A2)
512ω8
11Ḡ4γ3 − 258Ḡ2L̄2γ3 + 375L̄4γ3
+ 1024F 2L̄γω3 − 180Ḡ2L̄γ2ηω3 + 340L̄3γ2ηω3
+ 256F 2ηω6 − 48Ḡ2γη2ω6 + 176L̄2γη2ω6
+ 32L̄η3ω9 − 2ω3
17(10L̄2 − Ḡ2)γ2η
+ 96L̄γη2ω3 + 16η3ω6
L̄2 − Ḡ2 cos(2ḡ)
− 16(Ḡ2 − L̄2)γη2ω6 cos(4ḡ) + 256F 2ηω6 cos(2p̄)
− 512F 2γω3
L̄2 − Ḡ2 cos(2ḡ − 2p̄)
. (A3)
Consequently, the functions di(L̄, Ḡ) in equations (16)
are found to be
L̄2 − Ḡ2
)−1/2
d2 = −
(6L̄γη + 2η2ω3)(L̄2 − Ḡ2)−1/2,
256ω5
16η3ω6 + 51(4L̄2 − Ḡ2)γ2η
+ 96L̄γη2ω3
(L̄2 − Ḡ2)−1/2,
d4 = −
3Ḡγη2
3ḠF 2γ
(L̄2 − Ḡ2)−1/2,
d6 = −
(18L̄γ2 + 8γηω3),
256ω8
2γ3(11Ḡ2 − 129L̄2)− 180L̄γ2ηω3
− 48γη2ω6
. (A4)
Defining S = (L̄2 − Ḡ2)/Ḡ, the functions ej(L̄, Ḡ) in
equations (16) become
ej = −2djS, j = 1, . . . , 5, j 6= 3,
e3 = −
128ω5
16η3ω6 + 17(10L̄2 − Ḡ2)γ2η
+ 96L̄γη2ω3
L̄2 − Ḡ2. (A5)
APPENDIX B
The constant coefficients of equations (27) and (28) are
as follows
L̄(γL̄+ ηω3)
(9L̄γ2 + 4γηω3),
n2 = −
17γ2η
n4 = −
(59L̄2γ3 + 45L̄γ2ηω3 + 16γη2ω6),
n6 = −
3L̄γη + η2ω3
256ω4
(153L̄2γ2η + 96L̄γη2ω3 + 16η3ω6),
mi = −ni, i = 1, 2, 5,
m3 = −2n3,
(59L̄2γ3 + 45L̄γ2ηω3 + 8γη2ω6).
APPENDIX C
In this appendix we prove that K1(Q1, Q2) given in
(27), satisfies the condition (∗). To this end, we need the
following theorem.
Theorem 1. Any polynomial 1-form D of degree n in
Q1 and Q2 can be expressed as
D = dA+ rdK1 + ξ(K1)Q2dQ1, (C1)
where A(Q1, Q2) and r(Q1, Q2) are polynomials of degree
(n+1) and (n−1) respectively, and ξ(K1) is a polynomial
of degree [ 1
(n− 1)] where [x] denotes the greatest integer
in x.
Iliev [15] has proved the same theorem for H = (Q21 +
Q22)/2. Theorem 1 can thus be proved in a similar man-
ner. Here we only present a useful result.
Let D be a general polynomial 1-form of degree 1,
D = (a10Q1 + a01Q2 + a00)dQ1
+(b10Q1 + b01Q2 + b00)dQ2, (C2a)
then in (C1) we have
A(Q1, Q2) =
Q21 + b10Q1Q2 +
+ a00Q1 + b00Q2, (C2b)
r(Q1, Q2) = 0, (C2c)
ξ(K1) = a01 − b10. (C2d)
Since dK1 = 0 along any phase space orbit characterized
by K1(Q1, Q2) = k, and since the integral of an exact
differential dA around any closed curve is zero, from (C1)
we obtain
D = ξ(k)
Q2dQ1
= ξ(k)
Q2(τ)
On the other hand, from (25a) we have
= E(τ)− 2Q2,
where E(τ) is an even function of τ . Given the fact that
Q2 is an odd function of τ , we conclude that
D = −2ξ(k)
Q22dτ.
Consequently, if
D ≡ 0, it follows that ξ(k) ≡ 0 and
therefore D = dA+ rdK1, which completes the proof.
APPENDIX D
In equation (34), the constant coefficient I is
I = −
ηΩ2ω2
η(3γL̄− 4ηω3)(β + 4Ω2)
η(3β − 4Ω2)
4πη3/2ω4
[1] A. Raman and C. D. Mote Jr., Int. J. Non-Linear Mech.,
36, 261 (2001).
[2] M. A. Jalali and A. Angoshtari, Int. J. Non-Linear
Mech., 41 , 726 (2006).
[3] J. N. Reddy, Applied Functional Analysis and Variational
Methods in Engineering ( McGraw-Hill, New York, 1986).
[4] A. Raman and C. D. Mote Jr., Int. J. Non-Linear Mech.,
34, 139 (1999).
[5] J. Nowinski, ASME J. Appl. Mech., 72 (1964).
[6] N. Baddour and J. W. Zu, Appl. Math. Modeling, 25, 541
(2001).
[7] A. Deprit, Celest. Mech. Dyn. Astron., 51, 201 (1991).
[8] A. Deprit, Celest. Mech. Dyn. Astron., 1, 12 (1969).
[9] A. Deprit and A. Elipe, Celest. Mech. Dyn. Astron., 51,
227 (1991).
[10] V. K. Melnikov, Trans. Moscow Math., 12, 1 (1963).
[11] B. V. Chirikov, Physics Reports, 52, 263 (1979).
[12] J. Guckenheimer and P. Holmes, Nonlinear Oscillations,
Dynamical Systems, and Bifurcations of Vector Fields
(Springer, New York, 1983).
[13] J. P. Françoise, Ergod. Theory Dynam. Syst., 16, 87
(1996).
[14] L. Perko, Differential Equations and Dynamical Systems,
3rd edition (Springer, New York, 2001).
[15] I. D. Iliev, Math. Proc. Cambridge Phil. Soc., 127, 317
(1999).
[16] G. Contopoulos, Order and Chaos in Dynamical Astron-
omy (Springer, New York, 2002).

  We use a third-order perturbation theory and Melnikov's method to prove the
existence of chaos in spinning circular disks subject to a lateral point load.
We show that the emergence of transverse homoclinic and heteroclinic points
respectively lead to a random reversal in the traveling direction of
circumferential waves and a random phase shift of magnitude $\pi$ for both
forward and backward wave components. These long-term phenomena occur in
imperfect low-speed disks sufficiently far from fundamental resonances.

To appear in CHAOS, Volume 17, Issue 2, June 2007
On the existence of chaotic circumferential waves in spinning disks
Arzhang Angoshtari∗ and Mir Abbas Jalali†
Center of Excellence in Design, Robotics and Automation
Department of Mechanical Engineering, Sharif University of Technology
P.O.Box: 11365-9567, Azadi Avenue, Tehran, Iran
(Dated: August 4, 2021)
We use a third-order perturbation theory and Melnikov’s method to prove the existence of chaos
in spinning circular disks subject to a lateral point load. We show that the emergence of transverse
homoclinic and heteroclinic points respectively lead to a random reversal in the traveling direction
of circumferential waves and a random phase shift of magnitude π for both forward and backward
wave components. These long-term phenomena occur in imperfect low-speed disks sufficiently far
from fundamental resonances.
Keywords: Chaos, canonical perturbation theory, Melnikov’s method, spinning disks, traveling wave reversal
Transversal vibration modes of hard disk drives
(HDDs) are excited by the lateral aerodynamic
force of the magnetic head. Previous works [1, 2]
revealed that chaotic orbits are inevitable ingre-
dients of phase space flows when the lateral force
is large, or the disk is rotating near the critical
resonant speed. For low-speed disks, however, an
adiabatic invariant (a first integral) was found [2]
using a first-order averaging based on canonical
Lie transforms. According to the first-order the-
ory, regular vibrating modes of imperfect, low-
speed disks are independent of the angular veloc-
ity of the disk, Ω0. In such a circumstance, the
speed of circumferential waves is the natural fre-
quency of the lateral mode, ω, derived from linear
vibration analysis. HDDs are usually operated
with angular velocities smaller than ω (safely be-
low resonance). Moreover, the magnitude of the
lateral force F is very small. Given these condi-
tions, we show that it is impossible to continue
the Lie perturbation scheme [2] up to terms of
arbitrary order and remove the time variable t
from the Hamiltonian. In fact, due to the spe-
cial forms of nonlinearities in the dynamical equa-
tions of spinning disks, one can not remove t from
third-order terms. Subsequent application of a
second-order Melnikov theory reveals that trans-
verse homoclinic and heteroclinic points do exist
for all F,Ω0 6= 0. This implies chaos, or equiva-
lently, non-integrability of governing equations.
I. INTRODUCTION
Dynamics of continuum media, like fluids, rods, plates
and shells is usually formulated as a system of partial
differential equations (PDEs) for physical quantities in
∗Electronic address: arzhang˙a@mehr.sharif.edu
†Electronic address: mjalali@sharif.edu;
URL: http://sharif.edu/~mjalali
terms of the spatial coordinates x and the time t as
L(u) = 0. (1)
Here L is a nonlinear operator and u(x, t) is the vector
of dependent variables. When the boundary conditions
are somehow simple, approximate variational methods
based on modal decomposition and Galerkin’s projection
[3] can be used to reduce the order of governing equa-
tions. These methods begin with solving an auxiliary
eigenvalue problem, which is usually the variational field
equation δL(u,Λ) = 0, and build some complete basis
set Uk(x,Λk) for expanding u(x, t) in the spatial do-
main. Here Λk is an eigenvalue, which is characterized
by the vectorial index k. EachUk is called an eigenmode
or a shape function. Such a basis set should preferably
satisfy boundary conditions and be orthogonal.
Once a complete basis set is constructed, one may sup-
pose a solution of the form
u(x, t) =
Vk(t) ·Uk(x,Λk). (2)
Substituting from (2) into (1) and taking the inner prod-
L(u) ·Uk′(x)dx = 0, (3)
leaves us with a system of nonlinear ordinary differen-
tial equations (ODEs) for the amplitude functions Vk(t).
The evolution of the reduced ODEs shows the interaction
of different modes and their influence on the development
of spatiotemporal patterns.
In a series of papers, Raman and Mote [1, 4] used the
modal decomposition method to investigate transversal
oscillations of spinning disks whose deformation field is
described in terms of the displacement vector (u, v, w)
with u and v being the in-plane components. The most
important application of the spinning disk problem is in
design, fabrication and control of HDDs. The governing
PDEs for the evolution of displacement components were
first derived by Nowinski [5] and reformulated more re-
cently by Baddour and Zu [6]. Let us define (R, φ) as
http://arxiv.org/abs/0704.1406v1
mailto:arzhang_a@mehr.sharif.edu
mailto:mjalali@sharif.edu
http://sharif.edu/~mjalali
the usual polar coordinates. For a rotating disk with the
angular velocity Ω0, Nowinski’s theory assumes that the
in-plane inertias (u, v)Ω20, 2Ω0(u,t, v,t) and (u,tt, v,tt) are
ignorable againstRΩ20. This is a rough approximation for
high-speed disks and one needs to use the complete set
of equations as Baddour and Zu [6] suggest. Nowinski’s
theory, however, has its own advantages (like the exis-
tence of a stress function) that facilitate the study of the
most important transversal modes. In low-speed disks, or
disks with high flexural rigidity, one has Ω0 ≪ ω. Hence,
it is legitimate to apply Nowinski’s governing equations
in such systems.
In this paper we analytically prove the existence of
chaos, and therefore, non-integrability of the reduced
ODEs that govern the double-mode oscillations of im-
perfect spinning disks. We investigate low-speed disks
subject to a lateral point force exerted by the magnetic
head. The lateral force in HDDs is very small and its
origin is the aerodynamic force due to air flow in the gap
between the disk and the head. We show that chaotic
circumferential waves dominate some zones of the phase
space over the time scale t ∼ O(ǫ−3) with ǫ being a small
perturbation parameter. This indicates very slow evolu-
tion of random patterns, and the practical difficulties of
their identification.
The paper is organized as follows. In §II, we present
the Hamiltonian function in terms of Deprit’s [7] Lis-
sajous variables. In §III, we use a canonical perturba-
tion theory [8, 9] to eliminate the fast anomaly l from
the Hamiltonian. The action associated with l then be-
comes an adiabatic invariant. Transversal intersections
of destroyed invariant manifolds, and therefore, non-
integrability of the normalized equations, is proved by
a second-order Melnikov method in §IV. We present a
complete classification of circumferential waves in §V and
end up the paper with concluding remarks in §VI.
II. PROBLEM FORMULATION
Let us assume Um(R) as an orthogonal basis set that
represents the disk deformation in the radial direction.
The index m stands for the number of radial nodes
that Um(R) has. According to Raman and Mote’s [1]
treatment of imperfect disks, the following choice of the
transversal displacement field
w(R, φ, t) = Um(R)[x(t) cosnφ+ y(t) sinnφ], (4)
reduces Nowinski’s governing equations to a system of
ODEs for the amplitude functions x(t) and y(t) as
ẍ+ λ2x+ ǫγ
x2 + y2
x = ǫF cos(nΩ0t), (5a)
ÿ + ω2y + ǫγ
x2 + y2
y = ǫF sin(nΩ0t), (5b)
where λ, ω and γ are constant parameters that depend
on the geometry and material of the disk. ǫ is a small
perturbation parameter, F is the weighted integral of the
lateral point force, and Ω0 is the angular velocity of the
disk.
We suppose small deviations from perfect disks and
write the constant parameter of (5) as λ2/ω2 = 1 + ǫη.
We also define nΩ0 = ǫΩ with O(Ω) ∼ O(ω). Denoting
(px, py) as the momenta conjugate to (x, y), it can be
verified that equations (5) are derivable from the Hamil-
tonian function
p2x + p
x2 + y2
+ ΩP+
x2 + y2
)2 − F (x cos p+ y sin p)
. (6)
We have introduced the action P and its conjugate an-
gle p = ǫΩt to make our equations autonomous, which
is a preferred form for the application of canonical per-
turbation theories. The extended phase space has now
dimension six. Dynamics generated by (6) is better un-
derstood after carrying out a canonical transformation
(x, y, px, py) → (l, g, L,G) to the space of Lissajous vari-
ables [7] so that
x = s cos(g + l)− d cos(g − l), (7a)
y = s sin(g + l)− d sin(g − l), (7b)
px = −ω [s sin(g + l) + d sin(g − l)] , (7c)
py = ω [s cos(g + l) + d cos(g − l)] , (7d)
, d =
, L ≥ 0, |G| ≤ L.
In the space of Lissajous variables, the Hamiltonian de-
fined in (6) becomes
H = H0(L) + ǫH1(l, g, p, L,G, P ), (8)
H0 = ωL,
H1 = ΩP − F [s cos(g + l − p)− d cos(g − l − p)]
(s2 + d2)− 2sd cos(2l)
(s2 + d2) + s2 cos(2g + 2l)− 2sd cos(2g)
+ d2 cos(2g − 2l)− 2sd cos(2l)
From (8) we conclude that l is the fast angle, and g and
p are the slow ones. Therefore the long-term behavior
of the flows generated by (8) can be analyzed by averag-
ing H over l. After removing l, its corresponding action
L will be a constant of motion for the flows generated
by the averaged Hamiltonian 〈H〉l, and the phase space
dimension reduces from 6 to 4.
III. CANONICAL THIRD-ORDER AVERAGING
In order to average H over l, we use the normal-
ization procedure of Deprit and Elipe [9]. Denoting
X ≡ (l, g, p) and Y ≡ (L,G, P ), we define a Lie transfor-
mation (l, g, p, L,G, P ) → (l̄, ḡ, p̄, L̄, Ḡ, P̄ ) as
X = EW (X̄), Y = EW (Ȳ ), (9)
so that the Hamiltonian function in terms of the new
variables, K ≡ 〈H〉l, does not depend on l̄. EW is the Lie
transform generated by the function W and it is defined
EW (Z̄) = Z̄ + (Z̄;W ) +
((Z̄;W );W )
(((Z̄;W );W );W ) + · · · . (10)
In this equation, (f1; f2) denotes the Poisson bracket of
f1 and f2 over the (l̄, ḡ, p̄, L̄, Ḡ, P̄ )-space. We expand the
generating function W as
W = ǫW1 +
W2 + · · · , (11)
and specify the averaged, target Hamiltonian K =
K(ḡ, p̄, L̄, Ḡ, P̄ ) as the series [9]
K = K0 + ǫK1 +
K3 + · · · , (12)
K0 = ωL̄, (13a)
H1dl̄, (13b)
[2(H1;W1) + ((H0;W1);W1)]dl̄, (13c)
3(H1;W2) + 3((H1;W1);W1)
+ 2((H0;W2);W1) + ((H0;W1);W2)
+ (((H0;W1);W1);W1)
dl̄. (13d)
W1 and W2 are determined through solving the following
differential equations
= H1(l̄, ḡ, p̄, L̄, Ḡ, P̄ )−K1(ḡ, p̄, L̄, Ḡ, P̄ ),(14a)
= 2(H1;W1) + ((H0;W1);W1)−K2. (14b)
By substituting from (14) into (13) and evaluating the
integrals, one finds the explicit form of the new Hamilto-
nian K, which has been given in Appendix A up to the
third-order terms.
Once l̄ is removed from the Hamiltonian, L̄ becomes
an integral of motion. The slow dynamics of the system
is thus governed by the flows in the (ḡ, Ḡ)-space. We
introduce the slow time τ = p̄/Ω, ignore the fourth-order
terms in ǫ, and obtain the following differential equations
for the dynamics of (ḡ, Ḡ)
= f1(ḡ, Ḡ) + ǫh1(ḡ, Ḡ, ǫ, τ), (15a)
= −∂K
= f2(ḡ, Ḡ) + ǫh2(ḡ, Ḡ, ǫ, τ), (15b)
where
f1 = d6 + d1 cos(2ḡ), f2 = e1 sin(2ḡ),
h1 = d7 + d2 cos(2ḡ)
+ ǫ[d8 + d3 cos(2ḡ) + d4 cos(4ḡ) + d5 cos(2ḡ − 2Ωτ)],
h2 = e2 sin(2ḡ)
+ ǫ[e3 sin(2ḡ) + e4 sin(4ḡ) + e5 sin(2ḡ − 2Ωτ)]. (16)
In these equations, di (i = 1, · · · , 8) and ej (j = 1, · · · , 5)
are functions of L̄ and Ḡ (see Appendix A). It is re-
marked that the action P̄ appears only in K1 via the
term ΩP̄ . It then disappears in the normalized equa-
tions (15) after taking the partial derivatives of K with
respect to ḡ and Ḡ. The partial derivative of K with
respect to P̄ determines the evolution of p̄, which is in
accordance with the simple linear law p̄(τ) = Ωτ + p̄(0).
The dynamics of P̄ itself is governed by
= −∂K
= − 1
K(ḡ, Ḡ, τ). (17)
One may integrate (17) to obtain P̄ (τ) once equations
(15) are solved. The behavior of P̄ is thus inherited from
ḡ(τ) and Ḡ(τ).
IV. THE MELNIKOV FUNCTION
There are few analytical methods in the literature for
the detection of chaos in perturbed Hamiltonian systems
[10, 11]. Melnikov’s [10] method is the most powerful
technique when the governing equations take the form
= f(x) + ǫh(x, ǫ, τ), x ∈ R2, (18)
so that the unperturbed system dx/dτ = f(x) is in-
tegrable and possesses a homoclinic (heteroclinic) orbit
qh(τ) to a hyperbolic saddle point, and h(x, ǫ, τ) is T -
periodic in τ . The occurrence of chaos is examined by
the Melnikov function
M(τ0, ǫ) = ǫM1(τ0) + ǫ
2M2(τ0) + . . . , (19)
where Mk(τ0) denotes the kth-order Melnikov function.
Assume that Mi(τ0) is the first nonzero term, i.e.,
Mk(τ0) ≡ 0 for 1 ≤ k ≤ i − 1. If Mi(τ0) has simple
zeros, then, for sufficiently small ǫ, the system (18) has
transverse homoclinic (heteroclinic) orbits, which imply
chaos due to the Smale-Birkhoff homoclinic theorem [12].
The first-order term in (19) is determined by the classical
formula
M1(τ0) =
f (qh (τ)) ∧ h (qh (τ) , 0, τ + τ0) dτ,(20)
where the wedge operator ∧ is defined as f ∧ h =
f1h2 − f2h1. Although the Hamiltonian equations (15)
have a suitable form for the application of Melnikov’s
method, they are autonomous up to the first-order terms
in ǫ. Consequently, M1(τ0) vanishes identically for all
τ0 ∈ [0, T ]. We thus need to investigate the second-order
Melnikov function. For doing so, we begin with solving
the unperturbed system
= f1(ḡ, Ḡ), (21a)
= −∂K1
= f2(ḡ, Ḡ), (21b)
along homoclinic (heteroclinic) orbits. Jalali and An-
goshtari [2] showed that for L̄ > ηω3/γ, equations (21)
have hyperbolic stationary points at S0 ≡ (ḡ0, Ḡ0) =
(−π, 0), S2 ≡ (ḡ2, Ḡ2) = (0, 0), and S4 ≡ (ḡ4, Ḡ4) =
(π, 0). The implicit equation of the invariant manifolds
that terminate at the saddle points are
cos[2ḡh(τ)] =
L̄− γ
Ḡ2h(τ)
L̄2 − Ḡ2h(τ)
]−1/2
For γL̄ ≥ 2ηω3, equation (22) represents a heteroclinic
orbit which connects S0 to S2. For ηω
3 < γL̄ < 2ηω3,
the heteroclinic orbit disappears and it is replaced by a
homoclinic orbit (see Figure 1). To compute the explicit
form of the homoclinic (or heteroclinic) orbit of (21), we
use (21b) and (22), and obtain
∫ Ḡh(τ)
Ḡh(0)
1− αḠ2
βτ, (23)
16ηβω4
, β =
γL̄− ηω3
where the lower integration limit is Ḡh(0) = 1/
α. After
taking the integral (23), we arrive at
Ḡh(τ) =
sech(
β τ)√
, (24)
for the Ḡ ≥ 0 branch of the homoclinic (heteroclinic)
orbit. Having Ḡh(τ), it is straightforward to calculate
cos[2ḡh(τ)] and sin[2ḡh(τ)], and determine the explicit
form of qh(τ) =
ḡh(τ), Ḡh(τ)
For constructing M2(τ0), we use Françoise’s [13, 14]
algorithm that has been devised for dynamical systems
with polynomial nonlinearities. To express the averaged
Hamiltonian K in terms of polynomial functions of some
new dependent variables, we utilize Hopf’s variables
L̄2 − Ḡ2 cos(2ḡ), (25a)
L̄2 − Ḡ2 sin(2ḡ), (25b)
and obtain the following differential 1-form for the evo-
lution of the averaged system
dQ1 +
+ ǫ [(s1 + ǫs2) dQ2 − (z1 + ǫz2) dQ1] = 0. (26)
Here, the first-order Hamiltonian is
K1(Q1, Q2) =
(Q21 +Q
Q1 + C, (27)
s1 = m1Q2, (28a)
z1 = n1Q1 + n6, (28b)
s2 = m2(Q
2)Q2 +m3Q1Q2 +m4Q2
+ m5 sin(2Ωτ), (28c)
z2 = n2(Q
2)Q1 + n3(3Q
2) + n4Q1
+ n5 cos(2Ωτ) + n7. (28d)
The constant coefficients C, mi (i = 1, · · · , 5), and nj
(j = 1, · · · , 7) have been given in Appendix B. A prereq-
uisite for the application of Françoise’s [13] algorithm is
that for all polynomial 1-forms D that satisfy the condi-
D ≡ 0, (29)
there must exist polynomials A(Q1, Q2) and r(Q1, Q2)
such that D = dA+ rdK1. We call this the condition (∗)
and prove in Appendix C that K1 satisfies the condition
Françoise’s algorithm states that if M1(τ0) = · · · =
Mk−1(τ0) ≡ 0 for some integer k ≥ 2, it follows that
Mk(τ0) =
Dk, (30a)
D1 = δ1, Dm = δm +
i+j=m
riδj , (30b)
δj = zjdQ1 − sjdQ2, (30c)
for 2 ≤ m ≤ k. The functions ri are then determined
successively from the formulas Di = dAi + ridK1 for
i = 1, · · · , k − 1. We have already found that
M1(τ0) =
δ1 = 0, (31)
δ1 = (n1Q1 + n6)dQ1 −m1Q2dQ2. (32)
From (30) and (C2) it can be shown that
M2(τ0) =
δ2. (33)
Substituting (30c) and (28) into (33), and carrying out
the integration along qh(τ), result in
M2(τ0) =
3γF 2I
sin(2Ωτ0), (34)
with I being a constant (see Appendix D). Equation (34)
shows that τn = nπ/(2Ω) (n = 1, 2, · · · ) are simple zeros
of M2(τ0) so that
M2(τn) = 0,
∂M2(τ0)
τ0=τn
6= 0. (35)
FIG. 1: Possible topologies of the phase space flows of the averaged system for F = 0. In all panels the horizontal axis indicates
ḡ/π and the vertical axis indicates Ḡ/L̄. In the first topology (left panel), all stationary points are of center type. In the second
topology (middle panel) two new centers, with non-zero Ḡ-coordinates, have emerged for ḡ = ±nπ (n = 0, 1) and symmetrical
homoclinic loops (thick lines) connect saddle points to themselves. In the third topology (right panel) the off-axis centers (and
their surrounding tori) are still present but the separatrix curves (thick lines) are of heteroclinic type.
FIG. 2: Phase space structure of the normalized equations (15) with L̄ = 2, ω = η = F = 1, and Ω = ǫ = 0.1. Left panel:
γ = 0.8. Right panel: γ = 2. In both panels the horizontal and vertical axes indicate ḡ/π and Ḡ/L̄, respectively.
Thus, we conclude that the global stable and unstable
manifolds of the saddle point Sτ0n , W
s(Sτ0n ) andW
u(Sτ0n ),
always intersect transversely. Transversal intersections
cause a sensitive dependence on initial conditions due to
the Smale-Birkhoff homoclinic theorem. This is a route
to chaos. On the other hand this means that the reduced
equations (15) are non-integrable for F 6= 0.
V. CLASSIFICATION OF CIRCUMFERENTIAL
WAVES
For F = 0, h does not depend on τ and the normalized
equations (15) are integrable. In such a circumstance, the
phase space structure can take three general topologies
(depending on the values of the system parameters and
L̄) as shown in Figure 1. In the first topology all sta-
tionary points with the coordinates (ḡs, Ḡs) calculated
f(ḡs, Ḡs) + ǫh(ḡs, Ḡs, ǫ) = 0, (36)
are centers and they lie on the Ḡ = 0 axis with ḡs =
−π + nπ/2 (n = 0, · · · , 4). In the second and third
topologies, two off-axis centers (with Ḡs 6= 0) come to
existence for ḡs = ±nπ (n = 0, 1) and the on-axis sta-
tionary points with the same ḡs = ±nπ become saddles.
In the second topology, each saddle point is connected to
itself by a homoclinic orbit, and in the third topology, a
heteroclinic orbit connects two neighboring saddle points.
The system with heteroclinic orbits allows for rotational
ḡ(τ) while in the system with homoclinic orbits ḡ(τ) is
always librating. Beware that this classification of phase
space flows is valid as long as ǫ is sufficiently small.
For F = 0, the phase space flows of (15) are
structurally stable (with no unbounded branches) and
the whole (ḡ, Ḡ)-space is occupied by periodic orbits
of period T (K). At the stationary points, one has
T (K(ḡs, Ḡs, L̄)) = 0. Given the invariance of L̄, and
the periodic solutions ḡ(τ) = ḡ(τ + T (K)) and Ḡ(τ) =
Ḡ(τ + T (K)), the anomaly l̄ is determined through solv-
= ω + ǫ
, (37)
which results in l̄ = ωt + ǫR(τ) with R(τ) = R(τ +
T (K)). According to (9), the functions g(t), G(t) and
L(t) are also periodic in t and we conclude that l(t) =
ωt + RW (t) with RW (t) = EW (l̄) − ωt being a small-
amplitude periodic function of t. The explicit from of
the circumferential wave will then become
w(R, φ, t)
Um(R)
L(t)+G(t)
cos [nφ−ωt−RW (t)−g(t)]
L(t)−G(t)
cos [nφ+ωt+RW (t)−g(t)] , (38)
which is composed of a forward and a backward traveling
wave. Due to the periodic nature of L(t) and G(t), when
the amplitude of the forward traveling wave is maximum,
that of the backward wave is minimum and vice versa.
As our results of §IV shows, the regular nature of trav-
eling waves is destroyed for F 6= 0 and a chaotic layer oc-
curs through the destruction of the homoclinic and hete-
roclinic orbits of (21). This happens over the time scale
τ ∼ O(ǫ−2) or t ∼ O(ǫ−3) (because p̄ is present only
in K3). Figure 2 shows Poincaré maps of the system
(15) for F 6= 0. The sampling time step in generating
the Poincaré maps has been 2π/Ω. It is seen that most
tori around elliptic fixed points are preserved. They cor-
respond to regular periodic and quasi-periodic solutions
of the normalized system. For chaotic flows, the func-
tions ḡ(τ) and Ḡ(τ) randomly change within the invari-
ant measure of the chaotic set. Consequently, the orig-
inal Lissajous variables g(t), L(t), G(t) and also RW (t)
become chaotic too.
For Ḡ(τ) > 0 and Ḡ(τ) < 0 the forward and the back-
ward traveling waves are the dominant components of
the circumferential wave, respectively. When the chaotic
layer emerges from the destroyed homoclinic orbits (left
panel in Figure 2), the sign of Ḡ(τ) is randomly switched
along a chaotic trajectory. This means a random trans-
fer of kinetic/potential energy between the forward and
backward traveling wave components. For chaotic trajec-
tories of this kind the angle ḡ(τ) randomly fluctuates near
ḡ ≈ ±nπ (n = 0, 1) with an almost zero average. The
evolution of circumferential waves is quite different when
the chaotic layer emerges due to destroyed heteroclinic
orbits (right panel in Figure 2). In this case chaos means
a random change between the librational and rotational
states of ḡ(τ). Such a change induces an unpredictable
phase shift of magnitude π for both forward and back-
ward traveling wave components. We note that Ḡ(τ) can
flip sign on a chaotic trajectory only when ḡ(τ) is in its
librational state.
VI. CONCLUDING REMARKS
Resonance overlapping [11, 16] is the main cause for
chaotic behavior in spinning disks with near-resonant an-
gular velocities [1, 2]. The chaos predicted in this paper,
however, happens far from fundamental resonances. Op-
tical and HHDs are usually operated below critical res-
onant speeds and the lateral force F due to magnetic
head is very small. We showed that whatever the magni-
tude of F may be, a chaotic layer fills some parts of the
phase space because the Melnikov function of the nor-
malized equations has always simple zeros. Dynamics of
rotating disks is regular only if F vanishes, which is an
unrealistic assumption for disk drives. In low-speed disks
with small F , diffusion of chaotic orbits (within their in-
variant measure) takes a long time of t ∼ O(ǫ−3). The
slow development of chaotic circumferential waves makes
them undetectable in short time scales at which most
controllers work. The Melnikov function (34) depends
not only on F , but also on the parameter η through the
constant I. The parameter η is a contribution of imper-
fections, which are likely because of limited fabrication
precision in micro/nano scales. For a perfect disk with
η = 0, the off-axis elliptic stationary points of (21), and
consequently, homoclinic and heteroclinic orbits disap-
pear. In such a condition the Melnikov function is indef-
inite, but the system admits an exact first integral and
the dynamics is governed by the Hamiltonian function
given in equation (11) of Jalali and Angoshtari [2].
One of the most important achievements of this work
was to unveil the fact that it is premature to truncate the
series of canonical perturbation theories before recording
the role of all participating variables. In systems with
non-autonomous governing ODEs (non-conservative sys-
tems), one must be cautious while removing a fast angle
through averaging schemes. The removal of the fast an-
gle may also wipe out time-dependent terms, up to some
finite orders of ǫ, and hide some essential information of
the underlying dynamical process. Strange irregular so-
lutions can indeed occur at any order and influence the
long term response of dynamical systems as we observed
for the spinning disk problem by keeping the third-order
terms.
Acknowledgments
We are indebted to the anonymous referee, who discov-
ered an error in the early version of the paper and led us
to investigate the second-order Melnikov function. MAJ
thanks the Research Vice-Presidency at Sharif University
of Technology for partial support.
APPENDIX A: THE NORMALIZED
HAMILTONIAN
By evaluating the integrals in (13), we obtain the first,
second and third order terms of the normalized Hamilto-
nian as
−γḠ2
L̄(3γL̄+ 2ηω3)
L̄2 − Ḡ2 cos(2ḡ), (A1)
K2 = −
2L̄(−9Ḡ2 + 17L̄2)γ2 + 64F 2ω3
+ 8(3L̄2 − Ḡ2)γηω3 + 8L̄η2ω6
+ 4ω3[−6L̄γη − 2η2ω3]
L̄2 − Ḡ2 cos(2ḡ)
, (A2)
512ω8
11Ḡ4γ3 − 258Ḡ2L̄2γ3 + 375L̄4γ3
+ 1024F 2L̄γω3 − 180Ḡ2L̄γ2ηω3 + 340L̄3γ2ηω3
+ 256F 2ηω6 − 48Ḡ2γη2ω6 + 176L̄2γη2ω6
+ 32L̄η3ω9 − 2ω3
17(10L̄2 − Ḡ2)γ2η
+ 96L̄γη2ω3 + 16η3ω6
L̄2 − Ḡ2 cos(2ḡ)
− 16(Ḡ2 − L̄2)γη2ω6 cos(4ḡ) + 256F 2ηω6 cos(2p̄)
− 512F 2γω3
L̄2 − Ḡ2 cos(2ḡ − 2p̄)
. (A3)
Consequently, the functions di(L̄, Ḡ) in equations (16)
are found to be
L̄2 − Ḡ2
)−1/2
d2 = −
(6L̄γη + 2η2ω3)(L̄2 − Ḡ2)−1/2,
256ω5
16η3ω6 + 51(4L̄2 − Ḡ2)γ2η
+ 96L̄γη2ω3
(L̄2 − Ḡ2)−1/2,
d4 = −
3Ḡγη2
3ḠF 2γ
(L̄2 − Ḡ2)−1/2,
d6 = −
(18L̄γ2 + 8γηω3),
256ω8
2γ3(11Ḡ2 − 129L̄2)− 180L̄γ2ηω3
− 48γη2ω6
. (A4)
Defining S = (L̄2 − Ḡ2)/Ḡ, the functions ej(L̄, Ḡ) in
equations (16) become
ej = −2djS, j = 1, . . . , 5, j 6= 3,
e3 = −
128ω5
16η3ω6 + 17(10L̄2 − Ḡ2)γ2η
+ 96L̄γη2ω3
L̄2 − Ḡ2. (A5)
APPENDIX B
The constant coefficients of equations (27) and (28) are
as follows
L̄(γL̄+ ηω3)
(9L̄γ2 + 4γηω3),
n2 = −
17γ2η
n4 = −
(59L̄2γ3 + 45L̄γ2ηω3 + 16γη2ω6),
n6 = −
3L̄γη + η2ω3
256ω4
(153L̄2γ2η + 96L̄γη2ω3 + 16η3ω6),
mi = −ni, i = 1, 2, 5,
m3 = −2n3,
(59L̄2γ3 + 45L̄γ2ηω3 + 8γη2ω6).
APPENDIX C
In this appendix we prove that K1(Q1, Q2) given in
(27), satisfies the condition (∗). To this end, we need the
following theorem.
Theorem 1. Any polynomial 1-form D of degree n in
Q1 and Q2 can be expressed as
D = dA+ rdK1 + ξ(K1)Q2dQ1, (C1)
where A(Q1, Q2) and r(Q1, Q2) are polynomials of degree
(n+1) and (n−1) respectively, and ξ(K1) is a polynomial
of degree [ 1
(n− 1)] where [x] denotes the greatest integer
in x.
Iliev [15] has proved the same theorem for H = (Q21 +
Q22)/2. Theorem 1 can thus be proved in a similar man-
ner. Here we only present a useful result.
Let D be a general polynomial 1-form of degree 1,
D = (a10Q1 + a01Q2 + a00)dQ1
+(b10Q1 + b01Q2 + b00)dQ2, (C2a)
then in (C1) we have
A(Q1, Q2) =
Q21 + b10Q1Q2 +
+ a00Q1 + b00Q2, (C2b)
r(Q1, Q2) = 0, (C2c)
ξ(K1) = a01 − b10. (C2d)
Since dK1 = 0 along any phase space orbit characterized
by K1(Q1, Q2) = k, and since the integral of an exact
differential dA around any closed curve is zero, from (C1)
we obtain
D = ξ(k)
Q2dQ1
= ξ(k)
Q2(τ)
On the other hand, from (25a) we have
= E(τ)− 2Q2,
where E(τ) is an even function of τ . Given the fact that
Q2 is an odd function of τ , we conclude that
D = −2ξ(k)
Q22dτ.
Consequently, if
D ≡ 0, it follows that ξ(k) ≡ 0 and
therefore D = dA+ rdK1, which completes the proof.
APPENDIX D
In equation (34), the constant coefficient I is
I = −
ηΩ2ω2
η(3γL̄− 4ηω3)(β + 4Ω2)
η(3β − 4Ω2)
4πη3/2ω4
[1] A. Raman and C. D. Mote Jr., Int. J. Non-Linear Mech.,
36, 261 (2001).
[2] M. A. Jalali and A. Angoshtari, Int. J. Non-Linear
Mech., 41 , 726 (2006).
[3] J. N. Reddy, Applied Functional Analysis and Variational
Methods in Engineering ( McGraw-Hill, New York, 1986).
[4] A. Raman and C. D. Mote Jr., Int. J. Non-Linear Mech.,
34, 139 (1999).
[5] J. Nowinski, ASME J. Appl. Mech., 72 (1964).
[6] N. Baddour and J. W. Zu, Appl. Math. Modeling, 25, 541
(2001).
[7] A. Deprit, Celest. Mech. Dyn. Astron., 51, 201 (1991).
[8] A. Deprit, Celest. Mech. Dyn. Astron., 1, 12 (1969).
[9] A. Deprit and A. Elipe, Celest. Mech. Dyn. Astron., 51,
227 (1991).
[10] V. K. Melnikov, Trans. Moscow Math., 12, 1 (1963).
[11] B. V. Chirikov, Physics Reports, 52, 263 (1979).
[12] J. Guckenheimer and P. Holmes, Nonlinear Oscillations,
Dynamical Systems, and Bifurcations of Vector Fields
(Springer, New York, 1983).
[13] J. P. Françoise, Ergod. Theory Dynam. Syst., 16, 87
(1996).
[14] L. Perko, Differential Equations and Dynamical Systems,
3rd edition (Springer, New York, 2001).
[15] I. D. Iliev, Math. Proc. Cambridge Phil. Soc., 127, 317
(1999).
[16] G. Contopoulos, Order and Chaos in Dynamical Astron-
omy (Springer, New York, 2002).
