{'fraction_non_alphanumeric': 0.040961371527777776, 'fraction_numerical': 0.02318431712962963, 'mean_word_length': 4.346838135757107, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 42, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 30, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Paradoxes of dissipation-induced destabilization or who opened Whitney's umbrella?

The paradox of destabilization of a conservative or non-conservative system by small dissipation, or Ziegler's paradox (1952), has stimulated an ever growing interest in the sensitivity of reversible and Hamiltonian systems with respect to dissipative perturbations. Since the last decade it has been widely accepted that dissipation-induced instabilities are closely related to singularities arising on the stability boundary. What is less known is that the first complete explanation of Ziegler's paradox by means of the Whitney umbrella singularity dates back to 1956. We revisit this undeservedly forgotten pioneering result by Oene Bottema that outstripped later findings for about half a century. We discuss subsequent developments of the perturbation analysis of dissipation-induced instabilities and applications over this period, involving structural stability of matrices, Krein collision, Hamilton-Hopf bifurcation and related bifurcations.


Introduction
'Il n'y a de nouveau que ce qui est oublié'-this paraphrase of the Ecclesiastes 1:10, attributed to Marie-Antoinette, perfectly summarizes the story of the mathematical description of the destabilizing effect of vanishing dissipation in non-conservative systems.
There is a fascinating category of mechanical and physical systems which exhibit the following paradoxical behavior: when modeled as systems without damping they possess stable equilibria or stable steady motions, but when small damping is introduced, some of these equilibria or steady motions become unstable.
The paradoxical effect of damping on dynamic instability was noticed first for rotor systems which have stable steady motions for a certain range of speed, but which become unstable when the speed is changed to a value outside the range.
In 1879 Thomson and Tait [79] showed that a statically unstable conservative system which has been stabilized by gyroscopic forces could be destabilized again by the introduction of small damping forces. More generally, they consider conservative and nonconservative linear two degrees of freedom systems in remarkable detail. The destabilization by damping, using Routh's theorems, is implicit in their calculations, it is not formulated as a paradox.
In 1924, to explain the destabilization of a flexible rotor in stable rotation at a speed above the critical speed for resonance, Kimball [36] introduced a damping of the rotation, which has lead to non-conservative positional (circulatory) forces in the equations of motion of a gyroscopic system. In 1933 Smith [76] found that this non-conservative rotor system loses stability when the speed of rotation Ω > ω 1 + δ ν , where ω is the undamped natural whirling frequency (the critical speed for resonance) and δ and ν are the viscous damping constants for the stationary and rotating damping mechanisms. In Smith's model, the destabilizing effect of the damping of rotation (ν), observed also by Kapitsa [34], was compensated by the stationary damping (δ).
This was a first demonstration of a strong influence of the spatial distribution of damping (or equivalently the modal distribution) on the borderlines between stability and instability domains in multi-modal non-conservative systems [19,89].
In the 1950s and 1960s the publications of Ziegler [90,91], Bolotin [10,11], and Herrmann [25,26], motivated by aerodynamics applications, initiated a considerable activity in the investigation of dynamic instability of equilibrium configurations of structures under non-conservative loads. The canonical problem was the flutter of a vertical flexible cantilever column under a compressive non-conservative or follower load which remains tangent to the bending column. In the flutter mode the tip of the column is preponderantly slanted towards the left during the half-cycle in which the tip is moving towards the right and vice versa in the following half-cycle.
This snake-like oscillation permits the follower force to do positive work on each cycle [19].
The strong influence of the spatial or modal distribution of damping within the structure on its stability under non-conservative loading, observed in these publications, should not have been surprising in the light of earlier findings of rotor dynamists. However, they revealed explicitly the most dramatic and paradoxical aspect of the sensitivity of the stability of the nonconservative structures to small damping forces. It turned out that the critical load for a structure with small damping may be considerably smaller than that for the same structure without damping. In other words, there is a wide range of loads for which the undamped structure is stable, but which produce instability as soon as a tiny bit of damping is added to the structure.
These phenomena were actively studied in the 1960's to provide more basic understanding and they have continued to be studied with more sophisticated tools, including early attempts to employ singularity theory [81], until in the mid 1990s it was understood [30,71] that the destabilization paradox is related to the Whitney umbrella singularity of the stability boundary [84,85]. After describing in the first sections Whitney's umbrella and Ziegler's paradox, we make in section 4 a sharp turn to the 1950s to revisit an article of Oene Bottema [14], who in 1956 first made this discovery and clarified the paradox. Surprisingly, this paper surpassed the attention of most scientists during five decades.
In section 2 we will relate these results to singularity theory, in sections 5 and 7 we show in various ways their extension to finite-and infinite-dimensional systems using perturbation theory of multiple eigenvalues, in section 6 we focus on periodic systems, and in the remainder we discuss applications in physics and engineering.

Whitney's umbrella
In a remarkable paper of 1943 [84], Hassler Whitney described singularities of maps of a differential n-manifold into E m with m = 2n − 1. It turns out that in this case a special kind of singularity plays a prominent role. Later, the local geometric structure of the manifold near the singularity has been aptly called 'Whitney's umbrella'. In Fig. 1 we reproduce the original sketch of the singular surface from the companion article [85].
The paper contains two main theorems. Consider the C k map f : E n → E m with m = 2n − 1.
1. The map f can be altered slightly, forming f * , for which the singular points are isolated. For each such an isolated singular point p, a technical regularity condition C is valid which relates to the map f * of the independent vectors near p and of the differentials, the vectors in tangent space.

2.
Consider the map f * which satisfies condition C. Then we can choose coordinates x = (x 1 , x 2 , · · · , x n ) in a neighborhood of p and coordinates y = (y 1 , y 2 , · · · , y m ) (with m = 2n − 1) in a neighborhood of y = f (p) such that in a neighborhood of f * (p) we have exactly If for instance n = 2, m = 3, the simplest interesting case, we have near the origin so that y 1 ≥ 0 and on eliminating x 1 and x 2 : Starting on the y 2 -axis for y 1 = y 3 = 0, the surface widens up for increasing values of y 1 . For each y 2 , the cross-section is a parabola; as y 2 passes through 0, the parabola degenerates to a half-ray, and opens out again  Note that because of the C k assumption for the differentiable map f , the analysis is delicate. There is a considerable simplification of the treatment if the map is analytical.
The analysis of singularities of functions and maps is a fundamental ingredient for bifurcation studies of differential equations. After the pioneering work of Hassler Whitney and Marston Morse, it has become a huge research field, both in theoretical investigations and in applications. We can not even present a summary of this field here, so we restrict ourselves to citing a number of survey texts and discussing a few key concepts and examples. In particular we mention [3], [21], [22], [4], [2] and [5]. A monograph relating bifurcation theory with normal forms and numerics is [55].
The relation between singularities of functions and critical points or equilibria of differential equations becomes relatively simple when considering Hamiltonian and gradient systems. Consider for instance the timeindependent Hamiltonian function H(p, q) with p, q ∈ R n . Singularities of the function H are found in the set These points correspond with the critical points (equilibria) of the Hamiltonian equations of motioṅ More in general, consider the dynamical system described by the autonomous ODĖ An equilibrium x 0 of the system arises if f (x 0 ) = 0. With a little smoothness of the map f we can linearize near x 0 so that we can writeẋ with A a constant n × n− matrix, g(x) contains higher-order terms only. In other words The properties of the matrix A determine in a large number of cases the local behavior of the dynamical system. In a seminal paper [3], Arnold considers families of matrices, smoothly depending on a number of parameters (denoted by vector p). So, for the constant n × n− matrix we write A p . Suppose that for p = 0, A 0 is in Jordan normal form. Choosing p in a neighborhood of p = 0 produces a deformation (or perturbation) of A 0 , assuming that near p = 0 the entries of A p can be expanded in a convergent power series in the parameters.
A deformation is versal if all other deformations near p = 0 are equivalent under smooth change of parameters.
The paper [3] uses normal forms to obtain suitable versal deformations. These are associated with the bifurcations of the linearized system (3). Note that although a matrix induces a linear map, the corresponding eigenvalue problem produces a nonlinear characteristic equation. In addition, the parameters involved, make it necessary to analyze maps of R n into R m . For instance in the following sections we meet with maps from R 2 into R 3 as studied by Whitney [84]. Nevertheless, in 1943 it was hard to imagine that this study of global analysis, a pure mathematical abstraction, would find already an industrial application in the next decade.

Ziegler's paradox
In 1952 Hans Ziegler of ETH Zurich published a paper [90] that became classical and widely known in the community of mechanical engineers; it also attracted the attention of mathematicians. Studying a simplified twodimensional model of an elastic rod, fixed at one end and compressed by a tangential end load, he unexpectedly encountered a phenomenon with a paradoxal character: the domain of stability of the Ziegler's pendulum changes in a discontinuous way when one passes from the case in which the damping is very small to that where it has vanished [90,91].
Ziegler's double pendulum presented in Fig. 2(a) consists of two rigid rods of length l each, whose inclinations with respect to the vertical are denoted as ϕ 1 and ϕ 2 . Two masses m 1 and m 2 with the weights G 1 and G 2 are concentrated at the distances a 1 and a 2 from the joints. The elastic restoring torques and the damping torques at the joints are cϕ 1 , c(ϕ 2 − ϕ 1 ) and b 1φ1 , b 2 (φ 2 −φ 1 ), respectively. With these assumptions the kinetic energy of the system is while the potential energy reads The generalized dissipative and non-conservative forces are then Writing the Lagrange's equations of motionLφ i − L ϕi = Q i , where L = T − V and a dot denotes time differentiation, and assuming G 1 = 0 and G 2 = 0 for simplicity, we find   m 1 a 2 1 + m 2 l 2 m 2 la 2 m 2 la 2 m 2 a 2 With the substitution ϕ i = A i exp(λt), equation (7) yields a 4-dimensional eigenvalue problem with respect to the spectral parameter λ.
Putting m 1 = 2m, m 2 = m, a 1 = a 2 = l, b 1 = b 2 = b and assuming that dissipation is absent (b = 0), Ziegler found from the characteristic equation that the vertical equilibrium position of the pendulum looses its stability when the magnitude of the follower force exceeds the critical value P k , where In the presence of damping (b > 0) the Routh-Hurwitz condition yields the new critical follower load that depends on the square of the damping coefficient b Ziegler found that the domain of asymptotic stability for the damped pendulum is given by the inequalities P < P k (b) and b > 0 and he plotted it against the stability interval P < P k of the undamped system, Fig. 2(b).
Surprisingly, the limit of the critical load P k (b) when b tends to zero turned out to be significantly lower than the critical load of the undamped system Note that in the original work of Ziegler, formula (9) contains a misprint which yields linear dependency of the critical follower load on the damping coefficient b. Nevertheless, the domain of asymptotic stability found in [90] and reproduced in Fig. 2(b), is correct.
Some authors considered extensions of Ziegler's model by adding a conservative load and by assuming unequal damping coefficients [11,25,37,51,80]. Fig. 3 demonstrates how the domain of instability for the undamped Ziegler's pendulum with the partially follower load (η = 1 corresponds to the pure follower load), shown in dark gray in the (η, p)-plane, extends in a discontinuous manner in the presence of dissipation when b 2 = 0.3b 1 and b 1 → 0. The portion of the stability domain that became unstable is depicted in light gray [37,80]. Therefore, the two-dimensional stability diagrams of the undamped system and the system with vanishingly small damping differ by a region of positive measure.
Ziegler drew attention both to the substantial decrease in the critical load of the damped non-conservative system with vanishingly small dissipation and to the high sensitivity of the critical follower load with respect to the variation of the damping distribution. In the mechanical engineering literature these two effects are called the Ziegler's paradox of destabilization by small damping. In the conclusion to his classical book [10], Bolotin emphasized that the discrepancy between the stability domains of undamped non-conservative systems and that of systems with infinitesimally small dissipation is a topic of the greatest theoretical interest in stability theory. Encouraging further research of the destabilization paradox, Bolotin was not aware that the crucial ideas for its explanation were formulated as early as 1956.

Bottema's solution
In 1956, in the journal 'Indagationes Mathematicae', there appeared an article by Oene Bottema (1901Bottema ( -1992 [14], then Rector Magnificus of the Technical University of Delft and an expert in classical geometry and mechanics, that outstripped later findings for decades. Bottema's work in 1955 [13] can be seen as an introduction, it was directly motivated by Ziegler's paradox. However, instead of examining the particular model of Ziegler, he studied in [14] a much more general class of non-conservative systems. Following [13,14], we consider a holonomic scleronomic linear system with two degrees of freedom, of which the coordinates x and y are chosen in such a way that the kinetic energy is T = (ẋ 2 +ẏ 2 )/2. Hence the Lagrange equations of small oscillations near the equilibrium configuration x = y = 0 are as follows The matrices A and B can both be written uniquely as the sum of symmetrical and antisymmetrical parts: with k 11 = a 11 , k 22 = a 22 , k 12 = k 21 = (a 12 + a 21 )/2, ν = (a 12 − a 21 )/2 and The system (12) has a potential energy function (disregarding damping) when ν = 0, it is purely circulatory for k 11 = k 12 = k 22 = 0, it is non-gyroscopic for Ω = 0, and has no damping when d 11 = d 12 = d 22 = 0. If damping exists, we suppose in this section that it is positive.
In order to solve the equations (12) we put x = C 1 exp(λt), y = C 2 exp(λt) and obtain the characteristic equation for the frequencies of the small oscillations around equilibrium where [37,38,45] a 1 = trD, a 2 = trK + det D + Ω 2 , a 3 = trKtrD − trKD + 2Ων, For the equilibrium to be stable all roots of the characteristic equation (13) must be semi-simple and have real parts which are non positive.
It is always possible to write, in at least one way, the left hand-side as the product of two quadratic forms with real coefficients, Q = (λ 2 + p 1 λ + q 1 )(λ 2 + p 2 λ + q 2 ). Hence For all the roots of the equation (13) to be in the left side of the complex plane (L) it is obviously necessary and sufficient that p i and q i are positive. Therefore in view of (15) we have: a necessary condition for the roots Q = 0 having negative real parts is a i > 0 (i = 1, 2, 3, 4). This system of conditions however is not sufficient, that either one root of three roots lies in L (for then a 4 ≤ 0); it is also impossible that no root is in it (for, then a 4 ≤ 0). Hence if a i > 0 at least two roots are in L; the other ones are either both in L, or both on the imaginary axis, or both in R. In order to distinguish between these cases we deduce the condition for two roots being on the imaginary axis. If µi (µ = 0 is real) is a root, then µ 4 − a 2 µ 2 + a 4 = 0 and −a 1 µ 2 + a 3 = 0. Hence H := a 2 1 a 4 + a 2 3 − a 1 a 2 a 3 = 0. Now by means of (15) we have In view of a 1 > 0, a 3 > 0 the second factor is positive; furthermore a 1 = p 1 + p 2 > 0, hence p 1 and p 2 cannot both be negative. Therefore H < 0 implies p 1 > 0, p 2 > 0, for H = 0 we have either p 1 = 0 or p 2 = 0 (and not both, because a 3 > 0), for H > 0 p 1 and p 2 have different signs. We see from the decomposition of the polynomial (13) that all its roots are in L if p 1 and p 2 are positive.
Hence: a set of necessary and sufficient conditions for all roots of (13) to be on the left hand-side of the complex plane is We now proceed to the cases where all roots have non-positive real parts, so that they lie either in L or on the imaginary axis. If three roots are in L and one on the imaginary axis, this root must be λ = 0. Reasoning along the same lines as before we find that necessary and sufficient conditions for this are a i > 0 (i = 1, 2, 3), a 4 = 0, and H < 0. If two roots are in L and two (different) roots on the imaginary axis we have p 1 > 0, q 1 > 0, p 2 = 0, q 2 > 0 and the conditions are a i > 0 (i = 1, 2, 3, 4) and H = 0. If one root is in L and three are on the imaginary axis, then p 1 > 0, q 1 = 0, p 2 = 0, q 2 > 0 and the conditions are a i > 0 (i = 1, 2, 3), a 4 = 0, and The obtained conditions are border cases of (17). This does not occur with the last type we have to consider: all roots are on the imaginary axis. We now have p 1 = 0, p 2 = 0, q 1 > 0, q 2 > 0. Hence a 2 > 0, a 4 > 0, a 1 = a 3 = 0 and therefore H = 0. This set of relations is necessary, but not sufficient, as the example Q = λ 4 + 6λ 2 + 25 = 0 (which has two roots in L and two in the righthand side of the complex plane (R)) shows. The proof given above is not valid because as seen from (17), H = 0 does not imply now p 1 p 2 = 0, the second factor being zero for a 1 a 3 = 0 and q 1 = q 2 . The condition can of course easily be given; the equation (13) is λ 4 + a 2 λ 2 + a 4 = 0 and therefore it reads a 2 > 0, a 4 > 0, a 2 2 > 4a 4 .
Summing up we have: all roots of (13) (assumed to be different) have non-positive real parts if and only if one of the two following sets of conditions is satisfied [14] A : Note that a 1 represents the damping coefficients b 11 and b 22 in the system. One could expect B to be a limit of A, so that for a 1 → 0, a 3 → 0 the set A would continuously tend to B. That is not the case.
Furthermore, if A is satisfied and we take a 1 = εb 1 , a 3 = εb 3 , where b 1 and b 3 are fixed and ε → 0, the last condition of A reads (ε = 0) Obviously we have [14]  Following Bottema [14] we substitute in (13) λ = cµ, where c is the positive fourth root of a 4 > 0. The new A and B we get the same condition as when we write b i for a i , which was to be expected, because if the roots of (13) are outside R, those of P = 0 are also outside R and inversely. We can therefore restrict ourselves to the case a 4 = 1, so that we have only three parameters a 1 , a 2 , a 3 . We take them as coordinates in an orthogonal coordinate system.
The condition H = 0 or is the equation of a surface V of the third degree, which we have to consider for a 1 ≥ 0, a 3 ≥ 0, Fig. 4. Obviously The line is parallel to the 0a 1 a 3 -plane and intersects the a 2 -axis in a 1 = a 3 = 0, a 2 = m + 1/m ≥ 2. The a 2 -axis is the double line of V , a 2 > 2 being its active part. Two generators pass through each point of it; they coincide for a 2 = 2 (m = 1), and for a 2 → ∞ their directions tend to those of the a 1 and a 3 -axis (m = 0, m = ∞). The conditions A and B express that the image point (a 1 , a 2 , a 3 ) lies on V or above V . The point (0, 2, 0) is on V , but if we go to the a 2 -axis along the line a 3 = ma 1 the coordinate a 2 has the limit m + 1/m, which is > 2 but for m = 1. Curiously enough, even half a century later, there appear papers repeating this reasoning and the result almost literally, see for instance [68].  (19) Note that we started off with 8 parameters in eq. (11), but that the surface V bounding the stability domain is described by 3 parameters. It is described by a map of E 2 into E 3 as in Whitney's papers [84,85]. Explicitly, a transformation of (19) to (2) is given by Returning to the non-conservative system (11) (ν = 0), with damping, but without gyroscopic forces, so Ω = 0, and assuming as in [13] that k 12 = 0, k 11 > 0, and k 22 > 0 (a similar setting but with d 12 = 0 and k 12 = 0 was considered later by Bolotin in [10,12]), we find that the condition for stability H ≤ 0 reads Suppose now that the damping force decreases in a uniform way, so we put 22 are constants and ε → 0. Then, for the inequality (20) we get But if there is no damping, we have to make use of condition B, which gives Obviously where the expressions written in terms of the invariants of the matrices involved [45] are valid also without the restrictions on the matrices D and K that were adopted in [10,12,13]. For the values of 2trKD−trKtrD 2trD which are small with respect to ν 0 we can approximately write [39,40] If D depends on two parameters, say δ 1 and δ 2 , then (24) has a canonical form (2) for the Whitney's umbrella in the (δ 1 , δ 2 , ν)-space. Due to discontinuity existing for 2trKD − trKtrD = 0 the equilibrium may be stable if there is no damping, but unstable if there is damping, however small it may be. We observe also that the critical non-conservative parameter, ν cr , depends on the ratio of the damping coefficients and thus is strongly sensitive to the distribution of damping similarly to how it happens in rotor dynamics. This is the results which Ziegler [90,91] found in a special case.

'Hopf meets Hamilton under Whitney's umbrella'
The title of this section derives from a nice tutorial paper by Langford [56]. As we have seen, Bottema was the first who established that the asymptotic stability domain of a real polynomial of fourth order in the space of its coefficients consists of one of the 'pockets' of the Whitney umbrella. The corresponding singularity was later identified as generic in the three parameter families of real matrices by V.I. Arnold [3,4], who named it 'deadlock of an edge'. In this respect Bottema's results in [14] can be seen as an early study of bifurcations and structural stability of polynomials and matrices, and therefore of the singularities of their stability boundaries whose systematical treatment was initiated since the beginning of the 1970s in [3,4,57,58].
Although Bottema applied his result to nonconservative systems without gyroscopic forces, there are reasons for the singularity to appear in the case when gyroscopic forces are taken into account because the stability is determined by the roots of a similar fourth order characteristic polynomial. In order to study this case we consider separately the following m-dimensional version of the non-conservative system (11) x + (ΩG + δD)ẋ where a dot stands for time differentiation, x ∈ R m , and real matrix K = K T corresponds to potential forces.
Real matrices D = D T , G = −G T , and N = −N T are related to dissipative (damping), gyroscopic, and nonconservative positional (circulatory) forces with magnitudes controlled by scaling factors δ, Ω, and ν respectively.
A circulatory system, to which the undamped Ziegler's pendulum is attributed [40,62,70], is obtained from (25) by neglecting velocity-dependent forcesẍ while a gyroscopic one has no damping and non-conservative positional forces x + ΩGẋ + Kx = 0.
Circulatory and gyroscopic systems (26) and (27) possess fundamental symmetries that are evident after transformation of equation (25) to the formẏ = Cy with where I is the identity matrix.
In the absence of damping and gyroscopic forces (δ = Ω = 0), RCR = −C with This means that the matrix C has a time reversal symmetry, and equation (26) describes a reversible dynamical system [62]. Due to this property, and the eigenvalues of circulatory system (26) appear in pairs (−λ, λ). Without damping and non-conservative positional forces (δ = ν = 0) the matrix C possesses the Hamiltonian symmetry JCJ = C T , where J is a symplectic matrix [4,60,9] with As a consequence, which implies that if λ is an eigenvalue of C then so is −λ, similar to the reversible case. Therefore, an equilibrium of a circulatory or of a gyroscopic system is either unstable or all its eigenvalues lie on the imaginary axis of the complex plane, in the last case implying marginal stability if they are semi-simple.
It is well known that in the Hamiltonian case, the transition from gyroscopic stability to flutter instability occurs through the interaction of simple purely imaginary eigenvalues with the opposite Krein signature known as the Krein collision or the Hamiltonian Hopf bifurcation [23,56,59,60,61]. The collision occurs at the border of marginal stability, say at Ω = Ω 0 for (27), and it yields a double pure imaginary eigenvalue with the Jordan chain of vectors, which splits into a a complex conjugate pair under destabilizing variation of the parameter Ω.
Let iω 0 be the double eigenvalue at Ω = Ω 0 with the Jordan chain of generalized eigenvectors u 0 , u 1 , satisfying the equations [46] (−Iω 2 0 + iω 0 Ω 0 G + K)u 0 = 0, Then, the Krein collision in the gyroscopic system (27) is described by the following expressions where the real coefficient µ is according to [46] with the bar over a symbol denoting complex conjugate.
Perturbing the system (27) by small damping and circulatory forces yields an increment to a simple pure imaginary eigenvalue [40,46] With the expressions (34), equation (36) is used for the calculation of the deviation from the imaginary axis of the eigenvalues that participated in the Krein collision in the presence of the non-Hamiltonian perturbation that makes the merging of modes an imperfect one [28].
Since D and K are real symmetric matrices and N is a real skew-symmetric one, the first-order increment to the eigenvalue iω(Ω) given by (36) With the expansions (34) the formula (37) reads where we define From (38) it follows that in the vicinity of γ := ν/δ = γ * the limit of the critical value of the gyroscopic parameter Ω cr of the near-Hamiltonian system as δ → 0 exceeds the threshold of gyroscopic stabilization determined by the Krein collision (see [46]) Substituting γ = νδ in expression (40) yields a simple estimate for the critical value of the gyroscopic parameter Ω cr (δ, ν) that has a canonical form (2) and therefore describes the Whitney's umbrella surface in the (δ, ν, Ω)-space [46] Ω cr (δ, ν) = Ω 0 + n 2 In case of two oscillators (m = 2) the approximation (41) is transformed to [44,45,46] where ω 0 = 4 √ det K and Ω 0 = −trK + 2 √ det K in the assumption that det K > 0 and trK < 0. Due to the singularity the gyroscopic stabilization in the presence of dissipative and non-conservative positional forces depends on the ratio ν δ and is thus very sensitive to non-Hamiltonian perturbations. We will discuss gyroscopic stabilization in more detail in section 7.1.
We note that the sensitivity of simple eigenvalues of Hamiltonian and gyroscopic systems to dissipative perturbations was a subject of intensive investigations, see, e.g., MacKay [60], Haller [23], and Bloch et al. [9]. MacKay pointed out the necessity to extend such a perturbation analysis to multiple eigenvalues [60].
Maddocks and Overton [61] initiated the study of multiple eigenvalues and showed that for an appropriate class of dissipatively perturbed Hamiltonian systems, the number of unstable modes of the dynamics linearized at a nondegenerate equilibrium is determined solely by the index of the equilibrium regarded as a critical point of the Hamiltonian. They analyzed the movement of the eigenvalues in the limit of vanishing dissipation without direct application, however, to the destabilization paradox and approximation of the singular stability boundary. Our calculations performed in this section use the ideas developed in [44,45,46,50] that, however, can be traced back to the works of Andreichikov and Yudovich [1] and Crandall [19].
We see that in Hamiltonian mechanics, the Hamiltonian-Hopf bifurcation in which two pairs of complex conjugate eigenvalues approach the imaginary axis symmetrically from the left and right, then merge in double purely imaginary eigenvalues and separate along the imaginary axis (or the reverse) has codimension one. In the general case of non-Hamiltonian vector fields, the occurrence of double imaginary eigenvalues has codimension three. The interface between these two cases possesses the Whitney umbrella singularity; the Hamiltonian systems lie on its handle. Quoting Langford from his introductory paper [55] linking Hopf bifurcation, Hamiltonian mechanics and Whitneys umbrella: 'Hopf meets Hamilton under Whitney's umbrella', which, we add, was opened by Bottema. 6 Parametric resonance in systems with dissipation. x =ẋ = 0 stable or unstable? Another basic question is, what happens on adding damping effects? In the theory, certain resonance relations between the frequencies ω and ν play a crucial part. See for instance [4], [10], [73], [86] or [83] and Fig. 7(a) for this classical case.
In applications with parametric excitation where usually more degrees of freedom play a part, many combination resonances are possible. For a number of interesting cases, analysis and more references see [10,73]. In what follows, the so-called sum resonance will be important.
First we will consider the general procedure for systems with this combination resonance, after which we will discuss an application.

Normalization of oscillators in sum resonance
In [30] a geometrical explanation is presented for damping induced instability in parametric systems using 'all' the parameters of the system as unfolding parameters. It will turn out that, using symmetry and normalization, four parameters are needed to give a complete description in a two degrees of freedom system, or more generally systems where three frequencies are in resonance, but three parameters suffice to visualize the situation.
Consider the following type of nonlinear differential equation with three frequencieṡ which describes for instance a system of two parametrically forced coupled oscillators. A is a 4 × 4 matrix, containing a number of parameters, with purely imaginary eigenvalues ±iω 1 and ±iω 2 . Assume that A is semi-simple, so, if necessary, we can put A into diagonal form. The vector valued function f contains both linear and nonlinear terms and is 2π-periodic in ω 0 t, f (0, ω 0 t) = 0 for all t. Eq. (43) can be resonant in many different ways, but as announced, we consider here the sum resonance where the system may exhibit instability. The parameter δ is used to control the detuning δ = (δ 1 , δ 2 ) of the frequencies (ω 1 , ω 2 ) near resonance and the parameter µ = (µ 1 , µ 2 ) derives from the damping coefficients. So we may put A = A(δ, µ). We summarize the analysis from [30].
The basic approach will be to put eq. (43) into normal form by normalization or averaging whereas the theory from [3] will play a part. In the normalized equation the time-dependence is removed from lower order and appears only in the higher order terms. It turns out that the autonomous, linear part of this equation contains already enough information to determine the stability regions of small amplitude oscillations near the origin. The linear part of the normal form can be written aṡ and Figure 5: The critical surface in (µ + , µ − , δ + ) space for eq. (43). µ + = µ 1 + µ 2 , µ − = µ 1 − µ 2 , δ + = δ 1 + δ 2 .
Only the part µ + > 0 and δ + > 0 is shown. The parameters δ 1 , δ 2 control the detuning of the frequencies, the parameters µ 1 , µ 2 the damping of the oscillators (vertical direction). The base of the umbrella lies along the δ + -axis.
Since A(δ, µ) is the complexification of a real matrix, it commutes with complex conjugation. Furthermore, according to the normalization described in [4], [32] and [69] and if ω 1 and ω 2 are independent over the integers, the normal form of eq. (43) has a continuous symmetry group. The second step is then to test the linear part A(δ, µ) of the normalized equation for structural stability i.e. to answer the question whether there exist open sets in parameter space where the dynamics is qualitatively the same. The analysis follows [3] and [4]. The family of matrices A(δ, µ) is parameterized by the detuning δ and the damping µ. The procedure is to identify the most degenerate member N of this family, which turns out to be A(δ, 0) and then show that A(δ, µ) is its versal unfolding in the sense of [4]. The family A(δ, µ) is equivalent to a versal unfolding of the degenerate member N. For details we refer again to [30,83], an explicit example is discussed in the next subsection.
We can put the conclusions in a different way: the family A(δ, µ) is structurally stable for δ, µ > 0, whereas 0) is not. This has interesting consequences in applications as small damping and zero damping may exhibit very different behavior. In parameter space, the stability regions of the trivial solution are separated by a critical surface which is the hypersurface where A(δ, µ) has at least one pair of purely imaginary complex conjugate eigenvalues. As before, this critical surface is diffeomorphic to the Whitney umbrella, see Fig. 5.
It is the singularity of the Whitney umbrella that causes the discontinuous behavior displayed in the stability diagram in the subsection 6.3. The structural stability argument guarantees that the results are 'universally valid', i.e. they qualitatively hold for generic systems in sum resonance.
Above we have described the basic normalization approach, but if we are interested only in the shape of the resonance (instability) tongues, there are faster methods. For instance using the Poincaré-Linstedt method, see [83].
The parameter α is proportional to the rotation speed Ω. System (46) constitutes a conservative system of coupled Mathieu-like equations. Abbreviating P (t) = 4η 2 cos 2ηt, the corresponding Hamiltonian is: where p x , p y are the momenta. The natural frequencies of the unperturbed system (46), ε = 0, are ω 1 = √ α 2 + 1 + α and ω 2 = √ α 2 + 1 − α. By putting z = x + iy, system (46) can be written as: Introducing the new variable: and rescaling time ηt = τ , we obtain: where the prime denotes differentiation with respect to τ. By writing down the real and imaginary parts of this equation, we have actually got two identical Mathieu equations.
Using the classical and well-known results on the Mathieu equation, we conclude that the trivial solution is stable for ε small enough, provided that √ 1 + α 2 is not close to nη , for n = 1, 2, 3, .... The first-order and most prominent interval of instability, n = 1, arises if: If condition (50) is satisfied, the trivial solution of equation (49) is unstable. Therefore, the trivial solution of system (46) is also unstable. Note that this instability arises when: i.e. when the sum of the eigenfrequencies of the unperturbed system equals the excitation frequency 2η which is the sum resonance of first order. The domain of instability is bounded by: See Fig. 7(b) where the V-shaped instability domain is presented in the case of rotor rotation (α = 0) without damping.
Higher order combination resonances can be studied in the same way; the domains of instability in parameter space continue to narrow as n increases. As noted, the parameter α is proportional to the rotation speed Ω of the disk and also to the ratio of the moments of inertia.
Because of the damping term, we can no longer reduce the complex eq. (53) to two identical second order real equations, as we did previously.
In the sum resonance of the first order, we have ω 1 + ω 2 ≈ 2η and the solution of the unperturbed (ε = 0) equation can be written as: with Applying variation of constants leads to equations for z 1 and z 2 : To calculate the instability interval around the value η 0 = 1 2 (ω 1 + ω 2 ) = √ α 2 + 1, we apply normal form or (periodic solution) perturbation theory, see [66] for details, to find for the stability boundary: It follows that, as in other examples we have seen, the domain of instability actually becomes larger when damping is introduced. See Fig. 7b.
The instability interval, shows a discontinuity at κ = 0.
If κ → 0, then the boundaries of the instability domain tend to the limits which differs from the result we found when κ = 0 : η b = √ 1 + α 2 (1 ± ε). For reasons of comparison, we display the instability tongues in Fig. 7 in the four cases with and without rotation, with and without damping.
Mathematically, the bifurcational behavior is again described by the Whitney umbrella as indicated in subsection 6.1. In mechanical terms, the broadening of the instability-domain is caused by the coupling between the two degrees of freedom of the rotor in lateral directions which arises in the presence of damping. (b) The instability tongues for the rotor system. Again, because of damping the instability tongue is lifted off from the η-axis, but the tongue broadens. The boundaries of the V -shaped tongue without damping are to first approximation described by the expression η =

Manifestation of the destabilization paradox in other applications
In this section we discuss additional applications from physics and engineering, both finite-and infinitedimensional.

Gyroscopic systems of rotor dynamics
Investigation of the stability of equilibria of the Hauger's [24] and Crandall's [19,67] gyropendulums as well as of the Tippe Top [15,53] and the Rising Egg [15] leads to the system of linear equations known as the modified Maxwell-Bloch equations [9].
The modified Maxwell-Bloch equations are the normal form for rotationally symmetric, planar dynamical systems [9,15]. They follow from the equation (25) for m = 2, D = I, and K = κI, where κ corresponds to potential forces, and thus can be written as a single differential equation with complex coefficients According to (17) the solution x = 0 of equation (57) is asymptotically stable if and only if For κ > 0 the domain of asymptotic stability is a dihedral angle with the Ω-axis serving as its edge, Fig. 8(b).
Its sections by the planes Ω = const are contained in the angle-shaped regions with the boundaries At Ω = 0 the angle is bounded by the lines ν = ±δ √ κ and thus it is less than π. The domain of asymptotic stability is twisting around the Ω-axis in such a manner that it always remains in the half-space δ > 0, Fig. 8(b).
Consequently, the system that is statically stable at Ω = 0 and δ ≥ 0 can become unstable at greater Ω in the presence of non-conservative positional forces, as shown in Fig. 8(b) by the dashed line. The larger magnitudes of circulatory forces, the lower |Ω| at the onset of instability. This is a typical example of dissipation-induced instability in the sense of [9,52,53,54] when only non-Hamiltonian perturbations can cause the destabilizing movements of eigenvalues with definite Krein signature [48]. As κ > 0 decreases, the hypersurfaces forming the dihedral angle approach each other so that, at κ = 0, they temporarily merge along the line ν = 0 and a new configuration originates for κ < 0, Fig. 8(c). The new domain of asymptotic stability consists of two disjoint parts that are pockets of two Whitney's umbrellas singled out by inequality δ > 0. The absolute values of the gyroscopic parameter Ω in the stability domain are always not less than Ω 0 = 2 √ −κ. As a consequence, the system that is statically unstable at Ω = 0 can become asymptotically stable at greater Ω in the presence of circulatory forces, as shown in Fig. 8(c) by the dashed line.
As a mechanical example we consider Hauger's gyropendulum [24], which is an axisymmetric rigid body of mass m hinged at the point O on the axis of symmetry as shown in Fig. (8) Linearization of the nonlinear equations of motion derived in [24] with the new variables x 1 = ψ and x 2 = θ and the subsequent nondimensionalization yield the Maxwell-Bloch equations (57) where the dimensionless parameters are given by The domain of asymptotic stability of the Hauger gyropendulum, given by (58), is shown in Fig. 8(b,c).
For the statically unstable gyropendulum (κ < 0) the singular points on the Ω-axis correspond to the critical values ±Ω 0 = ±2 √ −κ and the critical frequency ω 0 = √ −κ. We find approximations of the stability boundary near the Whitney umbrella singularity as derived in [46,50]: Thus, Hauger's gyropendulum, which is statically unstable at Ω = 0, can become asymptotically stable for sufficiently large |Ω| ≥ Ω 0 under a suitable distribution of dissipative and nonconservative positional forces.
For almost all combinations of δ and ν the onset of gyroscopic stabilization of the non-conservative system is greater than that of a pure gyroscopic one (destabilization paradox: Ω cr (ν, δ) ≥ Ω 0 ). The obtained results are valid also for the equilibria of Tippe Top, Rising Egg, and Crandall's gyropendulum [44,45].

Circulatory systems of rotor dynamics
In some rotor dynamics applications gyroscopic effects are neglected [34,53]. For example, in the modeling of friction-induced oscillations in disc-and drum brakes, clutches and other machinery, the speed of rotation is assumed to be small. This frequently yields the linearized equations of motion in the form of a circulatory system with or without damping. In recent models the damping is included because it is believed that high sensitivity of the squeal onset to the damping distribution might be responsible for the poor reproducibility of the laboratory experiments with the squealing machinery.
Hoffmann and Gaul [28] studied a model of a mass sliding over a conveyor belt with friction and detected that small damping in this circulatory system destroys the reversible Hopf bifurcation and makes the collision of eigenvalues imperfect, exactly as it happens with the eigenvalues of Ziegler's pendulum [40,47].
In order to study squeal vibration in drum brakes Hultén [31,75] introduced a model shown in Fig. 9(a).
This model is composed of a mass m held against a moving band; the contact between the mass and the band is modeled by two plates supported by two different springs. It is assumed that the mass and band surfaces are always in contact and that the contact can be expressed by two cubic stiffnesses. Damping is included as shown in Fig. 9(a). The friction coefficient at contact is assumed to be constant and the band moves at a constant velocity. Then it is assumed that the direction of friction force does not change because the relative velocity between the band speed andẋ 1 orẋ 2 is assumed to be positive. The tangential force F T due to friction contact is assumed to be proportional to the normal force F N as given by Coulomb's law: F T = µF N . Assuming the Figure 9: (a) A model of a drum brake [31]; (b) Its numerically calculated critical friction coefficient at the onset of flutter instability as a function of damping parameters [75]; (c) A model of a disc brake [64]; (d) Its critical friction coefficient at the onset of flutter instability as a function of damping parameters calculated both numerically and perturbatively [47].
normal force F N is linearly related to the displacement of the mass normal to the contact surface, the resulting equations of motion can be expressed as  being exactly of the form considered by Bottema. Here the relative damping coefficients are denoted by η i = c i / √ m i k i (i = 1, 2) and natural pulsations are ω 0,i = k i /m i (i = 1, 2). Fig. 9(b) shows the numerically calculated domain of asymptotic stability of the drum brake in the space of the friction coefficient µ and two damping coefficients η 1 and η 2 with the Whitney umbrella singularity [75].
In Fig. 9(c) a model of a disc brake proposed in [64] is demonstrated. Its linearized equations of motion are again that of a circulatory system with small damping. It is not surprising that the critical friction coefficient at the onset of friction-induced vibrations as a function of two damping coefficients is represented in Fig. 9(d) by a surface with the Whitney umbrella singularity [47].
In both examples a selected distribution of damping exists that yields an increase in the critical load rather than decrease that happens for all other distributions. This possibility for stabilization was first pointed out in [71] for the Ziegler's pendulum. We will discuss this effect below in more detail.

Infinite-dimensional near-reversible and near-Hamiltonian systems
Dynamic instability, or flutter, is a general phenomenon which commonly occurs in coupled fluid-structure systems including pipes conveying fluids and airfoils [10,26,27]. Typically, the models are finite dimensional or continuous reversible systems that demonstrate the destabilization paradox in the presence of damping. In a recent study [87] Ziegler's paradox was observed in a problem of a vocal fold vibration (phonation) onset.

Near-reversible case: Beck's column with external and internal damping
Beck's column loaded by a follower force is a paradigmatic model for studying dynamical instability of structures. In 1969 Bolotin and Zhinzher [11] investigated the effects of damping distribution on its stability. They considered on the interval x ∈ [0, 1] the non-selfadjoint boundary eigenvalue problem of the form where λ is an eigenvalue with the eigenfunction u(x). The operators in the differential expression depend on the magnitude of the follower load q and the parameters of external, d 2 , and internal (Kelvin-Voight), and the vector u = (u(0), x u(1)) T . Some authors considered different boundary conditions that depend both on the physical parameters and on the spectral parameter [63,88] The undamped Beck's column is stable for q < q 0 ≃ 20.05 [17]. Stability is lost at q ≥ q 0 when after the reversible Hopf bifurcation the double pure imaginary eigenvalue iω 0 ≃ 11.02 splits into a pair of complex eigenvalues. In [11] it was found that in the presence of infinitesimally small Kelvin-Voight damping the critical load is reduced to q = q cr ≃ 10.94 and the critical frequency drops to ω = ω cr ≃ 5.40.
There were numerous attempts to find an approximation of the new critical load by studying the splitting of the double eigenvalue iω 0 of the unperturbed reversible system due to dissipative perturbations [70]. Banichuk et al. [6,7] have emphasized the importance of degenerate perturbations, the linear part of which is in the tangent plane to the Whitney umbrella singularity. Nevertheless, their analysis is not complete.
Further development of the approach of [6,7] in [37,40,41,42,43,74] resulted in the approximation to the critical load in the form where the vector of the damping parameters d=(d 1 , d 2 ) and angular brackets denote the scalar product in R 2 .
The components of the vector f and the real scalar f are and the components of the vector h are defined as and the real matrix G is defined by the expression whereŵ 2 is the solution of the boundary value problem The eigenfunctions u 0 and v 0 and the associated functions u 1 and v 1 of the original and adjoint eigenvalue problems are chosen to satisfy the bi-orthogonality and normalization conditions where the adjoint boundary value problems are connected by the Lagrange formula Formula (65) can serve for the approximation of the jump in the critical load. In the finite dimensional case with two degrees of freedom the expression for the limit of the critical load reduces to (24) [40]. For the Beck column described by the equations (62) we calculate the critical load as [41,42] q cr (d 1 , d 2 ) = q 0 − 1902d 2 1 (14.34d 1 + 0.091d 2 ) 2 + 12.68d 1 d 2 + 0.053d 2 2 .
Additionally, (65) captures more information on the geometry of the stability domain. For example one can plot the cross sections of the stability domain (73) for the different levels of q and find that some combinations of internal and external damping increase the critical load for the Beck's column, as shown in the central and right pictures of Fig. 10. The form of the stability boundary with the Whitney umbrella singularity approximated by equation (73) was confirmed later by numerical computations in [33]. The limit in the critical load following from (73) agrees well with the numerical data of [1]. changes in the stability domain. As a modern application we mention a moving beam with frictional contact investigated in [77]. Below we will consider an interesting example of the occurrence of the destabilization paradox in fluid dynamics.

Near-Hamiltonian case: The instability of baroclinic zonal currents
In the 1940s the first studies appeared of instability of baroclinic zonal (west-east) currents in the Earth's atmosphere [18,20]. It is remarkable that the unexpected destabilizing effect due to the introduction of dissipation was discovered in the linear stability analyses of this hydrodynamical problem by Holopainen (1961) [29] and Romea (1977) [65] at the very same period of active research on the destabilization paradox in structural mechanics. Recently these studies were revisited by Krechetnikov and Marsden [54] with the aim to handle rigorously the treatment of dissipation-induced instability.
Romea considered an infinite channel in the periodic zonal direction x of width L in the meridional direction y that is rotating with an angular velocity Ω. Two layers of incompressible, homogeneous fluids of slightly different densities (the lighter fluid on top) are confined by the side walls and by horizontal planes, a distance D apart. For simplicity, it is assumed that, in the absence of motion, the interface is located halfway between the horizontal planes, and is flat so that centrifugal effects may be ignored. Each layer moves downstream with a constant velocity and the slope of the interface is related to these velocities through the thermal wind relation.
It is implicitly assumed that this basic state is maintained against dissipation by an external energy source which is unimportant with respect to the rest of the problem [65].
The linearized equations for each layer near the basic state, characterized by the geostrophic streamfunctions −U 1 y and −U 2 y, are according to [65,54]: where F is the internal rotational Froude number, r ≥ 0 is the measure of the effect of Ekman suction (Ekman layer dissipation), and β is the planetary vorticity factor introduced to take into account the variation of the Coriolis parameter with latitude (β-effect).
Assuming the wave solutions ϕ 1,2 ∼ e iα(x−ct) sin(mπy), where real α > 0 is the x wavenumber, Romea obtained a dispersion relation for the complex phase speed c = c r + ic i in the form of the second-order complex polynomial. The real part of c is the speed of propagation of the perturbation, while αc i is the growth rate of the wave. If c i > 0, the wave grows, and the system is unstable.
In the inviscid case when the Ekman layer dissipation is set to zero, the transition to instability occurs through the Krein collision that occurs at U c := U 1 − U 2 = U cI , where [65,54] U cI = 2βF with a 2 = α 2 + m 2 π 2 . The critical shear U cI as a function of the wavenumber is plotted in Fig. 11(left). This curve bounds the region of marginal stability of the system without dissipation.
In the limit of vanishing viscosity (r → 0), the stability boundary differs from (75) The discrepancy between the stability domains of viscous and inviscid systems is clearly seen in Fig. 11(left).
Therefore, Romea demonstrated that an introduction of infinitesimally small dissipation destabilizes the system, lowering the curve of marginal stability by an O(1) amount. This is the appearance of the destabilization paradox in a continuous near-Hamiltonian system, which is similar to that found in near-reversible systems like Ziegler's pendulum (cf. Fig. 3) and Beck's column with dissipation [11,80]. Fig. 11(right) reproduces the original drawing from [65] showing the typical imperfect merging of modes [28] that substitutes the 'perfect' Krein collision in near-Hamiltonian and near-reversible systems. Approximation to the eigenvalue branches in imperfect merging can be efficiently calculated by means of the perturbation theory of multiple eigenvalues for a wide class of non-conservative systems [37,39,40,42,43,49].

Conclusion
We have revisited the pioneering result of Oene Bottema who in 1956 resolved the paradox of destabilization by small damping and interpreted it by means of what is now called the Whitney umbrella singularity. We have Figure 11: Original drawings from the 1977 work of Romea [65]: (left) Critical shear as a function of wavenumber demonstrates a discontinuous transition from the case when the Ekman layer dissipation r = 0 initially (U cI ) to the case when r → 0 (U cR ); (right) a typical imperfect merging of modes (growth rates) that substitutes the 'perfect' Krein collision in near-Hamiltonian and near-reversible systems and is characteristic for the destabilization paradox [11], [12].
shown that this phenomenon frequently occurs in near-Hamiltonian and near-reversible systems originating in very different areas of mechanics and physics ranging from hydrodynamics to contact mechanics and we have presented a unified treatment of it. There are a few related topics upon we did not touch. We mention interesting connections of this effect to structured pseudospectra [35] and to eigenvalue optimization problems [16]. We did not even consider the effect of nonlinearites. We mention the closely related effect of discontinuous change of the critical flutter frequency due to small dissipation [11,42] and its connection to the Whitney umbrella singularity at the exceptional points on the eigenvalue surfaces [72]. Another related topic is the role of the spectral exceptional points in modern non-Hermitian physics including crystal optics, open quantum systems, and PT -symmetric quantum mechanics [8]. All this shows that modern non-conservative and non-Hermitian problems are a perfect field of applied mathematics with a big potential for new discoveries.