From rotating fluid masses and Ziegler's paradox to Pontryagin- and Krein spaces and bifurcation theory

TL;DR

This paper explores the historical development and applications of Pontryagin and Krein spaces, linking classical fluid and mechanical systems to advanced mathematical theories in stability, bifurcation, and singularity analysis across various scientific fields.

Contribution

It provides a comprehensive overview of how classical mechanical systems inspired the development of indefinite metric spaces and their role in modern stability and bifurcation theory.

Findings

Connection between classical systems and Pontryagin/Krein space theory

Application of spectral singularities in astrophysics and engineering

Insight into dissipation-induced instabilities and eigenvalue collisions

Abstract

Three classical systems, the Kelvin gyrostat, the Maclaurin spheroids, and the Ziegler pendulum have directly inspired development of the theory of Pontryagin and Krein spaces with indefinite metric and singularity theory as independent mathematical topics, not to mention stability theory and nonlinear dynamics. From industrial applications in shipbuilding, turbomachinery, and artillery to fundamental problems of astrophysics, such as asteroseismology and gravitational radiation --- that is the range of phenomena involving the Krein collision of eigenvalues, dissipation-induced instabilities, and spectral and geometric singularities on the neutral stability surfaces, such as the famous Whitney's umbrella.

Equations

$$L=C_1+C_2-A_1-A_2-\frac{K}{{\Omega }^2},K=g(l_1{M}_1+l_2{M}_2).$$

$$\frac{dZ}{dt}$$

$$(A_1+\rho {\kappa }^2)\frac{dW}{dt}$$

$${\partial }_t{v}_x$$

$${\partial }_t{v}_y$$

$${\partial }_t{v}_z$$

$${\nabla }^2\chi =0,{{\partial }_n\chi |}_{\partial V}=z(\cos nx+i\cos ny)-(x+iy)\cos nz,$$

$$Q(R_1,R_2)=L\Omega Z_1{\overline{Z}}_2+(A_1+\rho {\kappa }^2)W_1{\overline{W}}_2+\frac{\rho }{2{\Omega }^2}{\int }_V{\overline{\text{v}}}_2^T{\text{v}}_{1}dV$$

$$C_2=\frac{8\pi \rho }{15}{a}^{4}c,A_2=l_2^2{M}_2+\frac{4\pi \rho }{15}{a}^{2}c(a^2+c^2),$$

$$v_x=(z-l_2)a^{2}m\xi ,v_y=-i(z-l_2)a^{2}m\xi ,v_z=-(x-iy)c^{2}m\xi ,$$

$$\frac{d\text{x}}{dt}=i\Omega {\text{A}}^{-1}\text{C}\text{x}=i\Omega \text{B}\text{x},$$

$$=$$

$$\textbf{G}=\left(\begin{array}[]{ccc}L&0&0\\ 0&A_{1}{+}l_{2}^{2}M_{2}{+}\frac{4\pi\rho}{15}a^{2}c\frac{(c^{2}-a^{2})^{2}}{c^{2}+a^{2}}&0\\ 0&0&\frac{4\pi\rho}{15}a^{4}c^{3}\frac{(c^{2}-a^{2})^{2}}{c^{2}+a^{2}}\\ \end{array}\right)$$

$$[\text{u},\text{u}]:=(\text{G}\text{u},\text{u})={\overline{\text{u}}}^T\text{G}\text{u},\text{u}\in {\mathbb{C}}^3,$$

$$(\lambda -\overline{\lambda }){\overline{\text{u}}}^T\text{G}\text{u}=0,$$

$${\lambda }_2^{\pm }=-\frac{1}{2}\pm \frac{1}{2}\sqrt{1+\frac{32\pi \rho }{15}\frac{c{a}^4}{L}},L=\frac{4\pi \rho }{15}{a}^{2}c(a^2-c^2),$$

$$\frac{x^2}{a_1^2}+\frac{y^2}{a_2^2}+\frac{z^2}{a_3^2}=1,a_3<a_2=a_1$$

$${\Omega }^2(e)=2e^{-3}(3-2e^2){\sin }^{-1}(e)\sqrt{1-e^2}-6e^{-2}(1-e^2).$$

$$\frac{x^2}{a_1^2}+\frac{y^2}{a_2^2}+\frac{z^2}{a_3^2}=1,a_3<a_2<a_1.$$

$${\mathbf{L}}_{i}(\lambda)=\lambda^{2}\left(\begin{array}[]{cc}1&0\\ 0&1\\ \end{array}\right)+\lambda\left(\begin{array}[]{cc}0&-4\Omega\\ \Omega&0\\ \end{array}\right)+\left(\begin{array}[]{cc}4b-2\Omega^{2}&0\\ 0&4b-2\Omega^{2}\\ \end{array}\right),$$

$$b=\frac{\sqrt{1-e^2}}{4e^5}\left\{e(3-2e^2)\sqrt{1-e^2}+(4e^2-3){\sin }^{-1}(e)\right\}.$$

$$\lambda =\pm \left(i\Omega \pm i\sqrt{4b-{\Omega }^2}\right).$$

$$4b(e)=2{\Omega }^2(e).$$

$$4b(e)={\Omega }^2(e),$$

$$e=\sin \left(\frac{e(3+4e^2)\sqrt{1-e^2}}{3+2e^2-4e^4}\right)$$

$${\Omega }^2(e)<4b(e)<2{\Omega }^2(e)$$

$$\left(\begin{array}[]{cc}4b-2\Omega^{2}&0\\ 0&4b-2\Omega^{2}\\ \end{array}\right)$$

$${\mathbf{L}}_{v}(\lambda)=\lambda^{2}\left(\begin{array}[]{cc}1&0\\ 0&1\\ \end{array}\right)+\lambda\left(\begin{array}[]{cc}10\mu&-4\Omega\\ \Omega&10\mu\\ \end{array}\right)+\left(\begin{array}[]{cc}4b-2\Omega^{2}&0\\ 0&4b-2\Omega^{2}\\ \end{array}\right),$$

$$25{\Omega }^2{\mu }^2+(\mathrm{R}e\lambda +5\mu {)}^2({\Omega }^2-\mathrm{R}e{\lambda }^2-10\mathrm{R}e\lambda \mu -4b)=0.$$

$${\text{L}}_g(\lambda )={\lambda }^2+\lambda (\text{G}+\text{D})+\text{K}+\text{N}$$

$$\textbf{G}=\frac{5}{2}\left(\begin{array}[]{rr}0&-\Omega\\ \Omega&0\\ \end{array}\right),\quad\textbf{D}=\left(\begin{array}[]{cc}\delta 16\Omega^{2}(6b-\Omega^{2})&-3\Omega/2\\ -3\Omega/2&\delta 16\Omega^{2}(6b-\Omega^{2})\\ \end{array}\right)$$

$$\textbf{K}=\left(\begin{array}[]{cc}4b-\Omega^{2}&0\\ 0&4b-\Omega^{2}\\ \end{array}\right),\quad\textbf{N}=\delta\left(\begin{array}[]{cc}2q_{1}&2q_{2}\\ -q_{2}/2&2q_{1}\\ \end{array}\right),$$

Full text

1 Historical background

The purpose of this paper is to show how a curious phenomenon observed in the natural sciences, destabilization by dissipation, was solved by mathematical analysis. After the completion of the analysis, the eigenvalue calculus of matrices which it involved was (together with other applications) an inspiration for the mathematical theory of structural stability of matrices. But the story is even more intriguing. Later it was shown that the bifurcation picture in parameter space is related to a seemingly pure mathematical object in singularity theory, Whitney’s umbrella.

In many problems in physics and engineering, scientists found what looked like a counter-intuitive phenomenon: certain systems that without dissipation show stable behaviour became unstable when any form of dissipation was introduced. This can be viscosity in fluids and solids, magnetic diffusivity, or losses due to radiation of waves of different nature, including recently detected gravitational waves Abb17 ; A2003 , to name a few. For instance, dissipation-induced modulation instabilities BD2007 are widely discussed in the context of modern nonlinear optics PTS2018 .

The destabilization by dissipation is especially sophisticated when several dissipative mechanisms are acting simultaneously. In this case, “no simple rule for the effect of introducing small viscosity or diffusivity on flows that are neutral in their absence appears to hold” TSL2013 and “the ideal limit with zero dissipation coefficients has essentially nothing to do with the case of small but finite dissipation coefficients” M1993 .

In hydrodynamics, a classical example is given by secular instability of the Maclaurin spheroids due to both fluid viscosity (Kelvin and Tait, 1879) and gravitational radiation reaction (Chandrasekhar, 1970), where the critical eccentricity of the meridional section of the spheroid depends on the ratio of the two dissipative mechanisms and reaches its maximum, corresponding to the onset of dynamical instability in the ideal system, when this ratio equals 1 LD1977 . In meteorology this phenomenon is known as the ‘Holopäinen instability mechanism’ (Holopäinen, 1961) for a baroclinic flow when waves that are linearly stable in the absence of Ekman friction become dissipatively destabilized in its presence, with the result that the location of the curve of marginal stability is displaced by an order one distance in the parameter space, even if the Ekman number is infinitesimally small Ho61 ; WE2012 . For a baroclinic circular vortex with thermal and viscous diffusivities this phenomenon was studied by McIntyre in 1970 MI1970 .

In rotor dynamics, the generic character of the discontinuity of the instability threshold in the zero dissipation limit was noticed by Smith already in 1933 S33 . In mechanical engineering such a phenomenon is called Ziegler’s paradox, it was found in the analysis of a double pendulum with a nonconservative positional force with and without damping in 1952 Z52 ; Z53 . The importance of solving the Ziegler paradox for mechanics was emphasized by Bolotin Bo63 : “The greatest theoretical interest is evidently centered in the unique effect of damping in the presence of pseudo-gyroscopic forces, and in particular, in the differences in the results for systems with slight damping which then becomes zero and systems in which damping is absent from the start.” Encouraging further research of the destabilization paradox, Bolotin was not aware that the crucial ideas for its explanation were formulated as early as 1956.

Ziegler’s paradox was solved in 1956 by an expert in classical geometry and mechanics, Oene Bottema B56 . He formulated the problem of the stability of an equilibrium in two degrees-of-freedom (4 dimensions), allowing for gyroscopic and nonconservative positional forces. The solution by Bottema in the form of concrete analysis was hardly noticed at that time. Google Scholar gives no citation of Bottema’s paper in the first 30 years after 1956.

As mentioned above, a new twist to the treatment of the problem came from identifying geometric considerations independently introduced by Whitney in singularity theory with the bifurcation analysis of Bottema. Interestingly Whitney’s “umbrella singularity” predated Bottema’s analysis, it gives the right geometric picture but it ignores the stability questions of the dynamical context which is the essential question of Ziegler’s paradox. Later, in 1971-72, V.I. Arnold showed that the umbrella singularity is generic in parameter families of real matrices. This result links to stability by linearisation of the vector field near equilibrium.

The phenomenon of dissipation-induced instability and the Ziegler paradox raised important questions in mechanics and mathematics. For instance, what is the connection between the conservative and deterministic Hamiltonian systems and real systems involving both dissipation of energy and stochastic effects La03 ? Also it added a new bifurcation in the analysis of dynamical systems. An important consequence of the results is that in a large number of problem-fields one can now predict and characterise precisely this type of instability HR95 ; KV10 ; K2013dg ; KI2017 ; KM07 ; La03 ; LFA2016 ; SJ07 .

1.1 Stability of Kelvin’s gyrostat and spinning artillery shells filled with liquid

Already von Laue (1905) VL1905 , Lamb (1908) Lamb1908 and Heisenberg (1924) H1924 realized that dissipation easily destabilizes waves and modes of negative energy of an ideal system supported by rotating and translating continua S1958 ; S1960 . Williamson (1936) proposed normal forms for Hamiltonian systems allowing sorting stable modes of negative and positive energy according to their symplectic sign MK86 ; W1936 ; W1937 . However, further generalization — the theory of spaces with indefinite metric, or Krein K1983 and Pontryagin P44 spaces, — was directly inspired by the problem of stability of the Kelvin gyrostat G1880 ; WT1880 , having both astrophysical and industrial (and even military) applications.

Kelvin WT1880 experimentally demonstrated in 1880 that a thin-walled and slightly oblate spheroid completely filled with liquid remains stable if rotated fast enough about a fixed point, which does not happen if the spheroid is slightly prolate, Figure 1. In the same year this observation was confirmed theoretically by Greenhill G1880 , who found that rotation around the center of gravity of the top in the form of a weightless ellipsoidal shell completely filled with an ideal and incompressible fluid, is unstable when $a<c<3a$, where $c$ is the length of the semiaxis of the ellipsoid along the axis of rotation and the lengths of the two other semiaxes are equal to $a$ G1880 .

Quite similarly, bullets and projectiles fired from the rifled weapons can relatively easily be stabilized by rotation, if they are solid inside. In contrast, the shells, containing a liquid substance inside, have a tendency to turn over despite seemingly revolved fast enough to be gyroscopically stabilized. Motivated by such artillery applications, in 1942 Sobolev, then director of the Steklov Mathematical Institute in Moscow, considered stability of a rotating heavy top with a cavity entirely filled with an ideal incompressible fluid S2006 —a problem that is directly connected to the classical XIXth century models of astronomical bodies with a crust surrounding a molten core S1959 .

For simplicity, the solid shell of the top and the domain $V$ occupied by the cavity inside it, can be assumed to have a shape of a solid of revolution. They have a common symmetry axis where the fixed point of the top is located. The velocity profile of the stationary unperturbed motion of the fluid is that of a solid body rotating with the same angular velocity $\Omega$ as the shell around the symmetry axis.

Following Sobolev, we denote by $M_{1}$ the mass of the shell, $M_{2}$ the mass of the fluid, $\rho$ and $p$ the density and the pressure of the fluid, $g$ the gravity acceleration, and $l_{1}$ and $l_{2}$ the distances from the fixed point to the centers of mass of the shell and the fluid, respectively. The moments of inertia of the shell and the ‘frozen’ fluid with respect to the symmetry axis are $C_{1}$ and $C_{2}$, respectively; $A_{1}$ ($A_{2}$) stands for the moment of inertia of the shell (fluid) with respect to any axis that is orthogonal to the symmetry axis and passes through the fixed point. Let, additionally,

[TABLE]

The solenoidal $(\textrm{div}\,\textbf{v}=0)$ velocity field v of the fluid is assumed to satisfy the no-flow condition on the boundary of the cavity: $\left.\textbf{ v}_{n}\right|_{\partial V}=0$.

Stability of the stationary rotation of the top around its vertically oriented symmetry axis is determined by the system of linear equations derived by Sobolev in the frame $(x,y,z)$ that has its origin at the fixed point of the top and rotates with respect to an inertial frame around the vertical $z$-axis with the angular velocity of the unperturbed top, $\Omega$. If the real and imaginary part of the complex number $Z$ describe the deviation of the unit vector of the symmetry axis of the top in the coordinates $x$, $y$, and $z$, then these equations are, see e.g. S2006 ; Y97 :

[TABLE]

where $2\kappa^{2}=\int_{V}|\nabla\chi|^{2}dV$, $E=i\int_{V}\left(\partial_{x}\overline{\chi}\partial_{y}\chi-\partial_{y}\overline{\chi}\partial_{x}\chi\right)dV$, and the function $\chi$ is determined by the conditions

[TABLE]

with $n$ the absolute value of a vector n, normal to the boundary of the cavity.

Sobolev realized that some qualitative conclusions on the stability of the top can be drawn with the use of the bilinear form

[TABLE]

on the elements $R_{1}$ and $R_{2}$ of the space $\{R\}=\{Z,W,\textbf{v}\}$. The linear operator $B$ defined by Eqs. (1.1) that can be written as $\frac{dR}{dt}=iBR$ has all its eigenvalues real when $L>0$, which yields Lyapunov stability of the top. The number of pairs of complex-conjugate eigenvalues of $B$ (counting multiplicities) does not exceed the number of negative squares of the quadratic form $Q(R,R)$, which can be equal only to one when $L<0$. Hence, for $L<0$ an unstable solution $R=e^{i\lambda_{0}t}R_{0}$ can exist with $\textrm{ Im}\lambda_{0}<0$; all real eigenvalues are simple except for maybe one S2006 ; Y97 ; Kopachevskii2001 .

In the particular case when the cavity is an ellipsoid of rotation with the semi-axes $a$, $a$, and $c$, the space of the velocity fields of the fluid can be decomposed into a direct sum of subspaces, one of which is finite-dimensional. Only the movements from this subspace interact with the movements of the rigid shell, which yields a finite-dimensional system of ordinary differential equations that describes coupling between the shell and the fluid.

Calculating the moments of inertia of the fluid in the ellipsoidal container

[TABLE]

denoting $m=\frac{c^{2}-a^{2}}{c^{2}+a^{2}}$, and assuming the field $\textbf{ v}=(v_{x},v_{y},v_{z})^{T}$ in the form

[TABLE]

one can eliminate the pressure in Eqs. (1.1) and obtain the reduced model

[TABLE]

where $\textbf{x}=(Z,W,\xi)^{T}\in\mathbb{C}^{3}$ and

[TABLE]

The matrix $\textbf{B}\neq\textbf{B}^{T}$ in Eq. (5) after multiplication by a symmetric matrix

[TABLE]

yields a Hermitian matrix $\textbf{G}\textbf{B}=\overline{(\textbf{G}\textbf{B})}^{T}$, i.e. B is a self-adjoint operator in the space $\mathbb{C}^{3}$ endowed with the metric

[TABLE]

which is definite when $L>0$ and indefinite with one negative square when $L<0$. If $\lambda$ is an eigenvalue of the matrix B, i.e. $\textbf{B}\textbf{u}=\lambda\textbf{u}$, then $\overline{\textbf{u}}^{T}\textbf{G}\textbf{B}\textbf{u}=\lambda\overline{\textbf{u}}^{T}\textbf{G}\textbf{u}$. On the other hand, $\overline{\textbf{u}}^{T}(\overline{\textbf{G}\textbf{B}})^{T}\textbf{u}=\overline{\lambda}\,\overline{\overline{\textbf{u}}^{T}\textbf{G}\textbf{u}}=\overline{\lambda}\,{\overline{\textbf{u}}^{T}\textbf{G}\textbf{u}}$. Hence,

[TABLE]

implying $\overline{\textbf{u}}^{T}\textbf{G}\textbf{u}=0$ on the eigenvector u of the complex $\lambda\neq\overline{\lambda}$. For real eigenvalues $\lambda=\overline{\lambda}$ and $\overline{\textbf{u}}^{T}\textbf{G}\textbf{u}\neq 0$. The sign of the quantity $\overline{\textbf{u}}^{T}\textbf{G}\textbf{u}$ (or Krein sign) can be different for different real eigenvalues.

For example, when the ellipsoidal shell is massless and the supporting point is at the center of mass of the system, then $A_{1}=0$, $C_{1}=0$, $M_{1}=0$, $l_{2}=0$. The matrix B has thus one real eigenvalue $(\lambda_{1}^{+}=-1,~{}{\overline{\textbf{u}}_{1}^{+}}^{T}\textbf{G}\textbf{u}_{1}^{+}>0)$ and the pair of eigenvalues

[TABLE]

which are real if $L>0$ and can be complex if $L<0$. The latter condition together with the requirement that the radicand in Eq. (16) is negative, reproduces the Greenhill’s instability zone: $a<c<3a$ G1880 . With the change of $c$, the real eigenvalue $\lambda_{2}^{+}$ with ${\overline{\textbf{u}}_{2}^{+}}^{T}\textbf{G}\textbf{u}_{2}^{+}>0$ collides at $c=3a$ with the real eigenvalue $\lambda_{2}^{-}$ with ${\overline{\textbf{u}}_{2}^{-}}^{T}\textbf{G}\textbf{u}_{2}^{-}<0$ into a real double defective eigenvalue $\lambda_{d}$ with the algebraic multiplicity two and geometric multiplicity one. This Krein collision is illustrated in Figure 2. Note that ${\overline{\textbf{u}}_{d}}^{T}\textbf{G}\textbf{u}_{d}=0$, where $\textbf{u}_{d}$ is the eigenvector at $\lambda_{d}$.

Therefore, in the case of the ellipsoidal shapes of the shell and the cavity, the Hilbert space $\{R\}=\{Z,W,\textbf{v}\}$ of the Sobolev’s problem endowed with the indefinite metric $(L<0)$ decomposes into the three-dimensional space of the reduced model (5), where the self-adjoint operator $B$ can have complex eigenvalues and real defective eigenvalues, and a complementary infinite-dimensional space, which is free of these complications. The very idea that the signature of the indefinite metric can serve for counting unstable eigenvalues of an operator that is self-adjoint in a functional space equipped with such a metric, turned out to be a concept of a rather universal character possessing powerful generalizations that were initiated by Pontryagin in 1944 P44 and developed into a general theory of indefinite inner product spaces or Krein spaces Bognar1974 ; KM2014 ; K1983 . Relation of the Krein sign to the sign of energy or action has made it a popular tool for predicting instabilities in physics N1976 ; Z2016 .

1.2 Secular instability of the Maclaurin spheroids by viscous and radiative losses

It is hard to find a physical application that would stimulate development of mathematics to such an extent as the problem of stability of equilibria of rotating and self-gravitating masses of fluids. Rooted in the Newton and Cassini thoughts on the actual shape of the Earth, the rigorous analysis of this question attracted the best minds of the XVIII-th and XIX-th centuries, from Maclaurin to Riemann, Poincaré and Lyapunov. In fact, modern nonlinear dynamics AA88 ; IA92 ; Kuz04 and Lyapunov stability theory L1992 are by-products of the efforts invested in solution of this question of the astrophysical fluid dynamics BKM2009 , which experiences a revival nowadays A2003 inspired by the recent detection of gravitational waves Abb17 .

We recall that in 1742 Maclaurin established that an oblate spheroid

[TABLE]

is a shape of relative equilibrium of a self-gravitating mass of inviscid fluid in a solid-body rotation about the $z$-axis, provided that the rate of rotation, $\Omega$, is related to the eccentricity $e=\sqrt{1-\frac{a_{3}^{2}}{a_{1}^{2}}}$ through the formula M1742

[TABLE]

A century later, Jacobi (1834) has discovered less symmetric shapes of relative equilibria in this problem that are tri-axial ellipsoids

[TABLE]

Later on a student of Jacobi, Meyer (1842) M1842 , and then Liouville (1851) L1851 have shown that the family of Jacobi’s ellipsoids has one member in common with the family of Maclaurin’s spheroids at $e\approx 0.8127$. The equilibrium with the Meyer-Liouville eccentricity is neutrally stable.

In 1860 Riemann R1860 established neutral stability of inviscid Maclaurin’s spheroids on the interval of eccentricities $(0<e<0.952..)$. At the Riemann point with the critical eccentricity $e\approx 0.9529$ the Hamilton-Hopf bifurcation sets in and causes dynamical instability with respect to ellipsoidal perturbations beyond this point DS1979 ; S1980-2 ; S1980-3 .

A century later Chandrasekhar C1969 used a virial theorem to reduce the problem to a finite-dimensional system, which stability is governed by the eigenvalues of the matrix polynomial

[TABLE]

where $\Omega(e)$ is given by the Maclaurin law (17) and $b(e)$ is as follows

[TABLE]

The eigenvalues of the matrix polynomial (18) are

[TABLE]

Requiring $\lambda=0$ we can determine the critical Meyer-Liouville eccentricity by solving with respect to $e$ the equation C1969

[TABLE]

The critical eccentricity at the Riemann point follows from requiring the radicand in (20) to vanish:

[TABLE]

which is equivalent to the equation

[TABLE]

that has a root $e\approx 0.9529$.

Remarkably, when

[TABLE]

both eigenvalues of the stiffness matrix

[TABLE]

are negative. The number of negative eigenvalues of the matrix of potential forces is known as the Poincaré instability degree. The Poincaré instability degree of the equilibria with the eccentricities (21) is even and equal to 2. Hence, the interval (21) corresponding to $0.812..<e<0.952..$, which is stable according to Riemann, is, in fact, the interval of gyroscopic stabilization L2013 of the Maclaurin spheroids, Figure 3.

In 1879 Kelvin and Tait TT realized that viscosity of the fluid can destroy the gyroscopic stabilization of the Maclaurin spheroids: “If there be any viscosity, however slight, in the liquid, the equilibrium [beyond $e\approx 0.8127$] in any case of energy either a minimax or a maximum cannot be secularly stable”.

The prediction made by Kelvin and Tait TT has been rigorously verified only in the XX-th century by Roberts and Stewartson RS1963 . Using the virial approach by Chandrasekhar, the linear stability problem can be reduced to the study of eigenvalues of the matrix polynomial

[TABLE]

where $\mu=\frac{\nu}{a_{1}^{2}\sqrt{\pi G\rho}}$, $\nu$ is the viscosity of the fluid, $G$ is the universal gravitation constant, and $\rho$ the density of the fluid C1969 . The $\lambda$ and $\Omega$ are measured in units of $\sqrt{\pi G\rho}$. The operator ${\mathbf{L}}_{v}(\lambda)$ differs from the operator of the ideal system, ${\mathbf{L}}_{i}(\lambda)$, by the matrix of dissipative forces $10\lambda\mu\textbf{I}$, where I is the $2\times 2$ unit matrix.

The characteristic polynomial written for ${\mathbf{L}}_{v}(\lambda)$ yields the equation governing the growth rates of the ellipsoidal perturbations in the presence of viscosity:

[TABLE]

The right panel of Figure 4 shows that the growth rates (23) become positive beyond the Meyer-Liouville point. Indeed, assuming ${\mathrm{R}e}\lambda=0$ in (23), we reduce it to $50\mu^{2}(\Omega^{2}-2b)=0,$ meaning that the growth rate vanishes when $\Omega^{2}=2b$ no matter how small the viscosity coefficient $\mu$ is. But, as we already know, the equation $\Omega^{2}(e)=2b(e)$ determines exactly the Meyer-Liouville point, $e\approx 0.8127$.

It turns out, that the critical eccentricity of the viscous Maclaurin spheroid is equal to the Meyer-Liouville value, $e\approx 0.8127$, even in the limit of vanishing viscosity, $\mu\rightarrow 0$, and thus does not converge to the inviscid Riemann value $e\approx 0.9529$.

Viscous dissipation destroys the interaction of two modes at the Riemann critical point and destabilizes one of them beyond the Meyer-Liouville point, showing an avoided crossing in the complex plane, Figure 5(left).

Kelvin and Tait TT hypothesised that the instability, which is stimulated by the presence of viscosity in the fluid, will result in a slow, or secular, departure of the system from the unperturbed equilibrium of the Maclaurin family at the Meyer-Liouville point and subsequent evolution along the Jacobi family, as long as the latter is stable C1969 .

Therefore, a rotating, self-gravitating fluid mass, initially symmetric about the axis of rotation, can undergo an axisymmetric evolution in which it first loses stability to a nonaxisymmetric disturbance, and continues evolving along a non-axisymmetric family toward greater departure from axial symmetry; then it undergoes a further loss of stability to a disturbance tending toward splitting into two parts. Rigorous mathematical treatment of the validity of the fission theory of binary stars by Lyapunov and Poincaré has laid a foundation to modern nonlinear analysis. In particular, it has led Lyapunov to the formulation of a general theory of stability of motion BKM2009 .

In 1970 Chandrasekhar C1970 demonstrated that there exists another mechanism making Maclaurin spheroids unstable beyond the Meyer-Liouville point of bifurcation, namely, the radiative losses due to emission of gravitational waves. However, the mode that is made unstable by the radiation reaction is not the same one that is made unstable by the viscosity, Figure 5(right).

In the case of the radiative damping mechanism stability is determined by the spectrum of the following matrix polynomial C1970

[TABLE]

that contains the matrices of gyroscopic, G, damping, D, potential, K, and nonconservative positional, N, forces

[TABLE]

[TABLE]

where $\Omega(e)$ and $b(e)$ are given by equations (17) and (19). Explicit expressions for $q_{1}$ and $q_{2}$ can be found in C1970 . The coefficient $\delta=\frac{GMa_{1}^{2}(\pi G\rho)^{3/2}}{5c^{5}}$ is related to gravitational radiation reaction, $G$ is the universal gravitation constant, $\rho$ the density of the fluid, $M$ the mass of the ellipsoid, and $c$ the velocity of light in vacuum.

In 1977 Lindblom and Detweiler LD1977 studied the combined effects of the gravitational radiation reaction and viscosity on the stability of the Maclaurin spheroids. As we know, each of these dissipative effects induces a secular instability in the Maclaurin sequence past the Meyer-Liouville point of bifurcation. However, when both effects are considered together, the sequence of stable Maclaurin spheroids reaches past the bifurcation point to a new point determined by the ratio, $X=\frac{25}{2\Omega_{0}^{4}}\frac{\mu}{\delta}$, of the strengths of the viscous and radiation reaction forces, where $\Omega_{0}=0.663490...$.

Figure 6 shows the critical eccentricity as a function of the damping ratio in the limit of vanishing dissipation. This limit coincides with the inviscid Riemann point only at a particular damping ratio. At any other ratio, the critical value is below the Riemann one and tends to the Meyer-Liouville value as this ratio tends either to zero or infinity. Lindblom and Detweiler LD1977 correctly attributed the cancellation of the secular instabilities to the fact that viscous dissipation and radiation reaction cause different modes to become unstable, see Figure 5. In fact, the mode destabilized by the fluid viscosity is the prograde moving spherical harmonic that appears to be retrograde in the frame rotating with the fluid mass and the mode destabilized by the radiative losses is the retrograde moving spherical harmonic when it appears to be prograde in the inertial frame A2003 . It is known Lamb1908 ; N1976 that to excite the positive energy mode one must provide additional energy to the mode, while to excite the negative energy mode one must extract energy from the mode C95 ; Sa08 . The latter can be done by dissipation and the former by the non-conservative positional forces K07b ; K09 ; K2013dg . Both are present in the model by Lindblom and Detweiler LD1977 .

1.3 Brouwer’s rotating vessel

In 1918, Brouwer B18 published a simple model for the motion of a point mass near the bottom of a rotating vessel. For an English translation see B76 pp. 665-675. The shape of the vessel is described by a surface $S$, rotating around a vertical axis with constant angular velocity $\omega$. With the equilibrium chosen on the vertical axis at $(x,y)=(0,0)$ on $S$, the linearized equations of motion without dissipation are

[TABLE]

The constants $k_{1}$ and $k_{2}$ are the $x,y$-curvatures of $S$ at $(0,0)$, $g$ is the gravitational constant. Suppose that $k_{1}\geq k_{2}$.

**Stability and instability without friction

**Assuming there is no damping, there are the following three cases:

• •

$0<k_{2}<k_{1}$ (single-well at equilibrium).

Stability iff $0<\omega^{2}<gk_{2}$ (slow rotation) or $\omega^{2}>gk_{1}$ (fast rotation).

• •

$k_{2}<0$ and $k_{1}>0,\,k_{1}>-k_{2}$ (saddle at equilibrium).

Stability iff $\omega^{2}>gk_{1}$.

• •

$k_{2}<0$ and $k_{1}>0,\,k_{1}<-k_{2}$ (saddle).

If $3k_{1}+k_{2}<0$: instability.

If $3k_{1}+k_{2}>0$: stability if

[TABLE]

• •

$k_{1}<0,\,k_{2}<0$: instability.

**Stability of triangular libration points $L_{4}$ and $L_{5}$

**Brouwer’s rotating vessel model includes both a well with two positive curvatures $k_{1}$ and $k_{2}$ and a saddle with the curvatures of opposite signs. Remarkably, the latter case describes stability of triangular libration points $L_{4}$ and $L_{5}$ (discovered by Lagrange in 1772) in the restricted circular three-body problem of celestial mechanics. Indeed, the linearized equation for this problem is (1.3) where $\omega=1$,

[TABLE]

$\mu=\frac{m_{1}}{m_{1}+m_{2}}$, and $m_{1}$ and $m_{2}$ are the masses of the two most massive bodies (in comparison with each of them the mass of the third body is assumed to be negligible) A1970 ; G1843 ; R2002 . Since $k_{1}$ and $k_{2}$ are of opposite signs for all $0\leq\mu\leq 1$, the linear stability of the triangular Lagrange equilibriums is determined by the stability conditions (25) for the rotating saddle BB1994 . Note that the coefficients (1.3) satisfy the constraint $gk_{1}+gk_{2}=-1$. Intervals of intersection of this line with the narrow corners of the curvilinear triangle in Figure 7b correspond to the stable Lagrange points. After substitution of (1.3) into (7), we reproduce the classical condition for their linear stability first established by Gascheau in 1843 G1843

[TABLE]

Stability in a rotating saddle potential is a subject of current active discussion in respect with the particle trapping BB1994 ; BC2018 ; K2011b ; JMT2017 . An effect of slow precession of trajectories of the trapped particles in a rotating saddle potential KL2017 ; S1998 has recently inspired new works leading to the improvement of traditional averaging methods CL2018 ; YH2018 .

**Destabilization by friction

**Adding constant (Coulomb) damping, Brouwer B18 finds a number of instability cases. The equations of motion become in this case:

[TABLE]

The friction constants $c_{1},c_{2}$ are positive. The characteristic (eigenvalue) equation takes the form:

[TABLE]

with $a_{1}=c_{1}+c_{2}$,

$a_{2}=g(k_{1}+k_{2})+2\omega^{2}+c_{1}c_{2}$

$a_{3}=c_{1}(gk_{2}-\omega^{2})+c_{2}(gk_{1}-\omega^{2})$,

$a_{4}=(gk_{1}-\omega^{2})(gk_{2}-\omega^{2})$.

There are two cases that drastically change the stability (see also Bottema Bo76 ) :

• •

$0<k_{2}<k_{1}$ (single-well).

Stability iff $0<\omega^{2}<gk_{2}$.

The fast rotation branch $\omega^{2}>gk_{1}$ has vanished.

• •

$k_{2}<0$ and $k_{1}>0,\,k_{1}<-k_{2}$ (saddle).

The requirement $a_{4}>0$ produces $0<gk_{1}<\omega\texttwosuperior$. This is not compatible with $a_{3}<0$ so a saddle is always unstable with any size of positive damping.

Brouwer studied this model probably to use in his lectures. In a correspondence with O. Blumenthal and G. Hamel he asked whether the results of the calculations were known; see B76 pp. 677-686. Hamel confirmed that the results were correct and surprising, but there is no reference to older literature in this correspondence. See also Bo76 .

**Indefinite damping and PT-symmetry

**Brouwer’s model for the case of a rotating well with two positive curvatures has a direct relation to rotordynamics as it contains the Jeffcott rotor model, see for instance K2011b . Dissipation induces instabilities in this model at high speeds of rotation $\omega$ Bo63 ; C95 ; S33 , of course, under the assumption that the damping coefficients are both positive. But what happens if we relax this constraint and extend the space of damping parameters to negative values? It turns out that at low speeds $\omega$ the domain of asymptotic stability spreads to the area of negative damping, Figure 8. Even more, a part of the neutral stability curve belongs to the line $c_{1}+c_{2}=0$ where one of the damping coefficients is positive and the other one is negative. On this line the system is invariant under time and parity reversion transformation and is therefore $PT$-symmetric. Its eigenvalues are imaginary in spite of the presence of the loss (positive damping) and gain (negative damping) in the system. $PT$-symmetric systems with the indefinite damping can easily be realized in the laboratory experiments Kottos2011 . Recent study Z2018 shows their connection to complex G-Hamiltonian systems KI2017 ; YS75 .

2 Ziegler’s paradox

We already know that Greenhill’s analysis of stability of the Kelvin gyrostat G1880 inspired the works of Sobolev S2006 and Pontryagin P44 which have led to the development of the theory of spaces with indefinite inner product Bognar1974 ; K1983 . Another work of Greenhill G1883 ultimately brought about the famous Ziegler paradox Z52 . As Gladwell remarked in his historical account of the genesis of the field of nonconservative stability G1990 , “It was Greenhill who started the trouble though he never knew it.”

Motivated by the problem of buckling of propeller-shafts of steamers Greenhill (1883) analyzed stability of an elastic shaft of a circular cross-section under the action of a compressive force and an axial torque G1883 . He managed to find the critical torque that causes buckling of the shaft for a number of boundary conditions. For the clamped-free and the clamped-hinged shaft loaded by an axial torque the question remained open until Nicolai in 1928 reconsidered a variant of the clamped-hinged boundary conditions, in which the axial torque is replaced with the follower torque N1928 . The vector of the latter is directed along the tangent to the deformed axis of the shaft at the end point G1990 .

Nicolai had established that no nontrivial equilibrium configuration of the shaft exists different from the rectilinear one, meaning stability for all magnitudes of the follower torque. Being unsatisfied with this overoptimistic result, he assumed that the equilibrium method does not work properly in the case of the follower torque. He decided to study small oscillations of the shaft about its rectilinear configuration using what is now known as the Lyapunov stability theory L1992 that, in particular, can predict instability via eigenvalues of the linearized problem.

Surprisingly, it turned out that there exist eigenvalues with positive real parts (instability) for all magnitudes of the torque, meaning that the critical value of the follower torque for an elastic shaft of a circular cross-section is actually zero. Because of its unusual behavior, this instability phenomenon received a name “Nicolai’s paradox” G1990 ; LFA2016 ; N1928 .

In 1951-56 Hans Ziegler of the ETH Zürich re-considered the five original Greenhill problems with the Lyapunov approach and found that the shaft is unstable in cases of the clamped-free and the clamped-hinged boundary conditions for all values of the axial torque, just as in Nicolai’s problem with the follower torque Z53 . Moreover, the non-self-adjoint boundary eigenvalue problem for the Greenhill’s shaft with the axial torque turned out to be a Hermitian adjoint of the non-self-adjoint boundary eigenvalue problem for the Greenhill’s shaft with the follower torque Bo63 .

In 1952, inspired by the paradoxes of Greenhill-Nicolai follower torque problems, Ziegler introduced the notion of the follower force and published a paper Z52 that became widely known in the engineering community, in particular among those interested in theoretical mechanics. It was followed by a second paper Z53 that added more details.

Ziegler considered a double pendulum consisting of two rigid rods of length $l$ each. The pendulum is attached at one of the endpoints and can swing freely in a vertical plane; see Figure 9. The angular deflections with respect to the vertical are denoted by $\phi_{1},\phi_{2}$, two masses $m_{1}$ and $m_{2}$ resulting in the external forces $G_{1}$ and $G_{2}$ are concentrated at the distances $a_{1}$ and $a_{2}$ from the joints. At the joints we have elastic restoring forces of the form $c\phi_{1}$, $c(\phi_{2}-\phi_{1})$ and internal damping torques

[TABLE]

So if $b_{1}=b_{2}=0$ we have no dissipation. We impose a follower force $P$ on the lowest hanging rod, see Figure 9. We consider only the quadratic terms of kinetic and potential energy. With these assumptions the kinetic energy $T$ of the system is:

[TABLE]

A dot denotes differentiation with respect to time $t$. The potential energy $V$ reads:

[TABLE]

This leads to the generalised dissipative and non-conservative forces $Q_{1},Q_{2}$:

[TABLE]

Writing the Lagrange’s equations of motion $\dot{L}_{\dot{\varphi}_{i}}-L_{\varphi_{i}}=Q_{i}$, where $L=T-V$ and assuming $G_{1}=0$ and $G_{2}=0$ for simplicity, we find

[TABLE]

The stability analysis of equilibrium follows the standard procedure. With the substitution $\varphi_{i}=A_{i}\exp(\lambda t)$, equation (40) yields a $4$-dimensional eigenvalue problem with respect to the spectral parameter $\lambda$.

Putting $m_{1}=2m$, $m_{2}=m$, $a_{1}=a_{2}=l$, $b_{1}=b_{2}=b$ and assuming that internal damping 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

[TABLE]

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$

[TABLE]

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, Figure 9(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

[TABLE]

Note that in the original work of Ziegler, formula (47) 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 Z52 and reproduced in Figure 9(b), is correct.

Ziegler limited his original calculation to a particular mass distribution, $m_{1}=2m_{2}$, and took into account only internal damping in the joints, neglecting, e.g, the air drag (an external damping). Later studies confirmed that the Ziegler paradox is a generic phenomenon and exists at all mass distributions, both for internal and external damping BK2018 ; B2018 ; SLK2019 ; T2016 , see Figure 10.

After invention of robust methods of practical realization of follower forces BN2011 the Ziegler paradox was immediately observed in the recent laboratory experiments BK2018 ; B2018 . Nowadays follower forces find new applications in cytosceletal dynamics BD2016 ; DC2017 and acoustics of friction HG03 . In general, the interest to mathematical models involving nonconservative positional forces (known also as circulatory Z53 or curl BS2016 forces) is growing both in traditional areas such as energy harvesting and fluid-structure interactions HD92 ; PMB2017 and in rapidly emerging new research fields of optomechanics SD2017 and light robotics LR2017 .

3 Bottema’s analysis of Ziegler’s paradox

In 1956, in the journal ‘Indagationes Mathematicae’, there appeared an article by Oene Bottema (1901-1992) B56 , 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 on stability in 1955 B55 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 B56 a much more general class of non-conservative systems.

Following B55 ; B56 , 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=({\dot{x}}^{2}+{\dot{y}}^{2})/2.$ Hence the Lagrange equations of small oscillations near the equilibrium configuration $x=y=0$ are as follows

[TABLE]

where $a_{ij}$ and $b_{ij}$ are constants, ${\bf A}:=(a_{ij})$ is the matrix of the forces depending on the coordinates, ${\bf B}:=(b_{ij})$ of those depending on the velocities. If $\bf A$ is symmetrical and while disregarding the damping associated with the matrix $\bf B$, there exists a potential energy function $V=(a_{11}x^{2}+2a_{12}xy+a_{22}y^{2})/2$, if it is antisymmetrical, the forces are circulatory. When the matrix $\bf B$ is symmetrical, we have a non-gyroscopic damping force, which is positive when the dissipative function $(b_{11}x^{2}+2b_{12}xy+b_{22}y^{2})/2$ is positive definite. If $\bf B$ is antisymmetrical the forces depending on the velocities are purely gyroscopic.

The matrices $\bf A$ and $\bf B$ can both be written uniquely as the sum of symmetrical and antisymmetrical parts: ${\bf A}={\bf K}+{\bf N}$ and ${\bf B}={\bf D}+{\bf G}$, where

[TABLE]

with $k_{11}=a_{11}$, $k_{22}=a_{22}$, $k_{12}=k_{21}=(a_{12}+a_{21})/2$, $\nu=(a_{12}-a_{21})/2$ and $d_{11}=b_{11}$, $d_{22}=b_{22}$, $d_{12}=d_{21}=(b_{12}+b_{21})/2$, $\Omega=(b_{12}-b_{21})/2$.

Disregarding damping, the system (3) has a potential energy function when $\nu=0$, it is purely circulatory for $k_{11}=k_{12}=k_{22}=0$, it is non-gyroscopic for $\Omega=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 linear stability problem for equations (3) we put $x=C_{1}\exp(\lambda t)$, $y=C_{2}\exp(\lambda t)$ and obtain the characteristic equation for the frequencies of the small oscillations around equilibrium

[TABLE]

where K03a ; K03b ; K07a

[TABLE]

For the equilibrium to be stable all roots of the characteristic equation (51) 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=(\lambda^{2}+p_{1}\lambda+q_{1})(\lambda^{2}+p_{2}\lambda+q_{2}).$ Hence

[TABLE]

For all the roots of the equation (51) 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 (53) we have: a necessary condition for the roots of $Q=0$ having negative real parts is $a_{i}>0$ $(i=1,2,3,4)$. This system of conditions however is not sufficient, as the example $(\lambda^{2}-\lambda+2)(\lambda^{2}+2\lambda+3)=\lambda^{4}+\lambda^{3}+3\lambda^{2}+\lambda+6$ shows. But if $a_{i}>0$ it is not possible that either one root of three roots lies in $L$ (for then $a_{4}\leq 0$); it is also impossible that no root is in it (for, then $a_{4}\leq 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 $\mu i$ ($\mu\neq 0$ is real) is a root, then $\mu^{4}-a_{2}\mu^{2}+a_{4}=0$ and $-a_{1}\mu^{2}+a_{3}=0$. Hence $H:=a_{1}^{2}a_{4}+a_{3}^{2}-a_{1}a_{2}a_{3}=0$. Now by means of (53) we have

[TABLE]

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 (51) 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 (51) to be on the left hand-side of the complex plane is

[TABLE]

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 $\lambda=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 $H=0$.

The obtained conditions are border cases of (55). 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=\lambda^{4}+6\lambda^{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 (55), $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 (51) is $\lambda^{4}+a_{2}\lambda^{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 (51) (assumed to be different) have non-positive real parts if and only if one of the two following sets of conditions is satisfied B56

[TABLE]

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}\rightarrow 0$, $a_{3}\rightarrow 0$ the set $A$ would continuously tend to $B$. That is not the case.

Remark first of all that the roots of (51) never lie outside $R$ if $a_{1}=0$, $a_{3}\neq 0$ (or $a_{1}\neq 0$, $a_{3}=0$). Furthermore, if $A$ is satisfied and we take $a_{1}=\varepsilon b_{1}$, $a_{3}=\varepsilon b_{3}$, where $b_{1}$ and $b_{3}$ are fixed and $\varepsilon\rightarrow 0$, the last condition of $A$ reads $(\varepsilon\neq 0)$

[TABLE]

while for $\varepsilon=0$ we have

[TABLE]

Obviously we have B56

[TABLE]

so that $(g_{1}>g_{2})$ but for $b_{3}=b_{1}\sqrt{a_{4}}$. That means that in all cases where $b_{3}\neq b_{1}\sqrt{a_{4}}$ we have a discontinuity in our stability condition.

In 1987, Leipholz remarked in his monograph on stability theory L1987 that “Independent works of Bottema and Bolotin for second-order systems have shown that in the non-conservative case and for different damping coefficients the stability condition is discontinuous with respect to the undamped case.” However, Leipholz did not mention that, in contrast to Bolotin, Bottema illustrated the phenomenon of the discontinuity in a remarkable geometric diagram, first published in B56 and reproduced in Figure 11.

Following Bottema B56 we substitute in (51) $\lambda=c\mu$, where $c$ is the positive fourth root of $a_{4}>0$. The new equation reads $P:=\mu^{4}+b_{1}\mu^{3}+b_{2}\mu^{2}+b_{3}\mu+1=0$, where $b_{i}=a_{i}/c^{i}$ $(i=1,2,3,4)$. If we substitute $a_{i}=c^{i}b_{i}$ in $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 (51) 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

[TABLE]

is the equation of a surface $V$ of the third degree, which we have to consider for $a_{1}\geq 0$, $a_{3}\geq 0$, Figure 11. Obviously $V$ is a ruled surface, the line $a_{3}=ma_{1}$, $a_{2}=m+1/m$ $(0<m<\infty)$ being on $V$. 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\geq 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}\rightarrow\infty$ their directions tend to those of the $a_{1}$ and $a_{3}$-axis $(m=0,m=\infty)$. 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$.

Note that we started off with 8 parameters in Eq. (3), 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 W43 ; W44 . Explicitly, a transformation of (19) to (2) is given by

[TABLE]

with $w^{2}=\frac{1}{4}y_{3}^{2}+y_{1}y_{2}$.

Returning to the non-conservative system (3) $(\nu\neq 0)$, with damping, but without gyroscopic forces, so $\Omega=0$, and assuming as in B55 that $k_{12}=0$, $k_{11}>0$, and $k_{22}>0$ (a similar setting but with $d_{12}=0$ and $k_{12}\neq 0$ was considered later by Bolotin in Bo63 ), we find that the condition for stability $H\leq 0$ reads

[TABLE]

Suppose now that the damping force decreases in a uniform way, so we put $d_{11}=\varepsilon d_{11}^{\prime}$, $d_{12}=\varepsilon d_{12}^{\prime}$, $d_{22}=\varepsilon d_{22}^{\prime}$, where $d_{11}$, $d_{12}$, $d_{22}$ are constants and $\varepsilon\rightarrow 0$. Then, for the inequality (3) we get

[TABLE]

But if there is no damping, we have to make use of condition $B$, which gives

[TABLE]

Obviously

[TABLE]

where the expressions written in terms of the invariants of the matrices involved K07a are valid also without the restrictions on the matrices $\bf D$ and $\bf K$ that were adopted in Bo63 ; B55 . For the values of $\frac{2{\rm tr}{\bf K}{\bf D}-{\rm tr}{\bf K}{\rm tr}{\bf D}}{2{\rm tr}{\bf D}}$ which are small with respect to $\nu_{0}$ we can approximately write Ki04 ; K05

[TABLE]

If $\bf D$ depends on two parameters, say $\delta_{1}$ and $\delta_{2}$, then (62) has a canonical form (64) for the Whitney’s umbrella in the $(\delta_{1},\delta_{2},\nu)$-space. Due to discontinuity existing for $2{\rm tr}{\bf K}{\bf D}-{\rm tr}{\bf K}{\rm tr}{\bf D}\neq 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, $\nu_{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 other applications, including the viscous Chandrasekhar-Friedman-Schutz (CFS) instability of the Maclaurin spheroids LD1977 .

The analytical approximations of the type (62) for the onset of the flutter instability in the general finite-dimensional and infinite-dimensional cases were obtained for the first time in the works K03a ; K03b ; Ki04 ; K05 ; KS05b ; K07b as a result of further development of the sensitivity analysis of simple and multiple eigenvalues in multiparameter families of non-self-adjoint operators. The previous important works include AY75 ; BBM89a ; Br94 ; C95 ; Ha92 ; DS1979 ; J88 ; MK91 ; OR96 ; S1980-2 ; S1980-3 ; S96 ; SKM05 ; Sey . Recent results on the perturbation analysis of dissipation-induced instabilities and the destabilization paradox are summarized in the works K2013dg and LFA2016 .

4 An umbrella without dynamics

Part of global analysis, a topic of pure mathematics, is concerned with singularity theory, which deals with the geometric characterisation and classification of singularities (stationary points) of vector fields. In dynamics these singularities are recognised as equilibria of dynamical systems. Well-known representatives of this singularity school are René Thom Th and Christopher Zeeman. Among pure mathematicians they were exceptional as they promoted singularity theory as useful for real-life problems in biology, the social sciences and physics. Unfortunately their approach gave singularity theory a bad name as in their examples they used geometric methods without explaining a possible relation between realistic vector fields and dynamics. It makes little sense to describe equilibria and transitions (bifurcations) between equilibria without discussing actual causes that are tied in with dynamical processes and corresponding equations of motion. We want to stress here that notwithstanding the lack of dynamics the geometry of singularities as an ingredient of dynamical systems theory can be very useful.

Before Ziegler’s results a geometric result in singularity theory was obtained (1943-44) by Hassler Whitney. This result turned out to be an excellent complement to Bottema’s analytic approach. In his paper W43 , 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 Figure 12 we reproduce a sketch of the singular surface.

The paper W43 contains two main theorems. Consider the $C^{k}$ map $f:E^{n}\mapsto E^{m}$ with $m=2n-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. 2.

Consider the map $f^{*}$ which satisfies condition $C$. Then we can choose coordinates $x=(x_{1},x_{2},\cdots,x_{n})$ in a neighborhood of $p$ and coordinates $y=(y_{1},y_{2},\cdots,y_{m})$ (with $m=2n-1$) in a neighborhood of $y=f(p)$ such that in a neighborhood of $f^{*}(p)$ we have exactly

[TABLE]

If for instance $n=2,\,m=3$, the simplest interesting case, we have near the origin

[TABLE]

so that $y_{1}\geq 0$ and on eliminating $x_{1}$ and $x_{2}$:

[TABLE]

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 [math], the parabola degenerates to a half-ray, and opens out again (with sense reversed); see Figure 12.

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 AA88 ; Ar71 ; Ar83 ; AA93 ; GS85 ; GS88 ; G1997 ; F1992 . A monograph relating bifurcation theory with normal forms and numerics is Kuz04 .

5 Hopf bifurcation near 1:1 resonance and structural stability

A study of the stability of equilibria of dynamical systems will usually involve the analysis of matrices obtained by linearisation of the equations of motion in a neighbourhood of the equilibria. This triggered off the study of structural stability of matrices as an independent topic in singularity theory Ar71 ; Ar83 .

More explicitly, consider a dynamical system described by the autonomous ODE

[TABLE]

where ${\bf p}\in{\mathbb{R}}^{k}$ is a vector of parameters. An equilibrium ${\bf x}_{0}$ of the system arises if ${\bf f}({\bf x}_{0},{\bf p})={\bf 0}$. With a little smoothness of the map ${\bf f}$ we can linearise near ${\bf x}_{0}$ so that we can write

[TABLE]

with ${\bf A}_{\bf p}$ a constant $n\times n-$ matrix, ${\bf g}({\bf x},{\bf p})$ contains higher-order terms only. In other words

[TABLE]

${\bf g}({\bf x},{\bf p})$ is tangent to the linear map in ${\bf x}_{0}$. The properties of the matrix ${\bf A}_{\bf p}$ determine in a large number of cases the local behavior of the dynamical system.

Suppose that for ${\bf p}={\bf 0}$, ${\bf A}_{\bf 0}$ has two equal non-zero imaginary eigenvalues and their complex conjugates, $\pm i\omega,\omega>0$, and no other eigenvalues with zero real part. This equilibrium is called a $1:1$ resonant double Hopf point G1997 . (Similarly, in a Hopf point the matrix of linearization has a single conjugate pair of imaginary eigenvalues $\pm i\omega$ and in a double Hopf point there are two distinct such pairs: $\pm i\omega_{1}$, $\pm i\omega_{2}$ G1997 .) Then, without loss of generality, we may assume that the system (65) has been already reduced to a centre manifold of dimension $n=4$. Considering further a generic case of double non-semisimple eigenvalues with geometric multiplicity 1, after a linear change of coordinates and re-scaling time to get $\omega=1$, we can transform ${\bf A}_{\bf 0}$ to F1992

[TABLE]

Setting $z_{1}=\Delta x_{1}+i\Delta x_{2}$ and $z_{2}=\Delta x_{3}+i\Delta x_{4}$, where $i=\sqrt{-1}$ and $\Delta{\bf x}={\bf x}-{\bf x}_{0}$, and assuming ${\bf A}_{\bf 0}$ to be in the form (66) we re-write (65) at ${\bf p}={\bf 0}$ in the complex form F1992

[TABLE]

The second pair of equations governing the conjugates $\overline{z}_{1}$, $\overline{z}_{2}$ is omitted here for simplicity.

Arnold Ar71 ; Ar83 has proven that a universal unfolding of the linear vector field with the matrix

[TABLE]

is given by the three-parameter family of complex matrices

[TABLE]

where $\alpha$, $\mu_{1}$, and $\mu_{2}$ are real parameters and versality is understood with respect to the group of similarity transformations and a real positive scaling. The set of matrices with a resonant Hopf pair is a group orbit G1997 .

The universal unfolding has a pure imaginary eigenvalue if and only if there exists a real number $\delta$ such that $(i(\delta-1)-\alpha)^{2}-(\mu_{1}+i\mu_{2})=0$. Eliminating $\delta$ yields G1997

[TABLE]

Setting $y_{2}=\alpha$, $y_{3}=\mu_{2}$, and $y_{1}=−4\mu_{1}+4\alpha^{2}$ we reduce the equation to the form $y_{3}^{2}=y_{1}y_{2}^{2}$, which is nothing else but the normal form (64) for the Whitney umbrella. The double Hopf points of (68) form the half-line $\alpha=\mu_{2}=0$, $\mu_{1}=0$. Along the continuation $\mu_{1}>0$ of this half-line the eigenvalues of (68) are given by $\lambda=i\pm\sqrt{\mu_{1}}$. We see that the double Hopf points have codimension 2 and the resonant double Hopf points are of codimension 3.

If a family of matrices ${\bf A}({\bf p})={\bf A}(p_{1},p_{2},p_{3},\ldots,p_{k})$ has a $1:1$ resonant double Hopf point, the universality of the unfolding (68) means that there exist smooth functions $\alpha({\bf p})$, $\mu_{1}({\bf p})$, $\mu_{1}({\bf p})$, such that the Hopf structure of ${\bf A}({\bf p})$ near the $1:1$ resonant point is the same as the Hopf structure of the unfolding with $\alpha$, $\mu_{1}$, $\mu_{2}$ replaced by $\alpha({\bf p})$, $\mu_{1}({\bf p})$, $\mu_{2}({\bf p})$.

Therefore, the stratified set of Hopf points in the neighborhood of a non-semisimple $1:1$ resonance is a Whitney umbrella in $\bf p$-space too G1997 . The functions $\alpha({\bf p})$, $\mu_{1}({\bf p})$, $\mu_{2}({\bf p})$ can be found approximately as truncated Taylor series with respect to the components of the vector $\bf p$ of the parameters Sey .

Similar stratification of Hopf points near a semisimple $1:1$ resonance involves pairs of Whitney umbrellas forming a Plücker conoid HK10 ; HK14 , see Figure 13.

Hoveijn and Ruijgrok (1995) were the first who applied these ideas to a practical problem exhibiting the Ziegler paradox. Namely, they considered a problem of widening due to dissipation of the zones of the combination resonance YS75 in a system of two parametrically forced coupled oscillators HR95 . The system models a rotating disk with oscillating suspension point introduced in RTV93 . Its linearized equations are

[TABLE]

It is assumed that for $\varepsilon=0$ the system (5) has two pairs of imaginary eigenvalues $\pm i\omega_{1}$, $\pm i\omega_{2}$ that depend on the parameter $\Omega$ representing the speed of rotation. Of special interest is the case of the sum resonance $\omega_{0}=\omega_{1}+\omega_{2}$.

Let parameters $\delta_{1}$ and $\delta_{2}$ control the detuning of the frequencies $\omega_{1}$ and $\omega_{2}$; then $\delta_{+}=\delta_{1}+\delta_{2}$ and $\delta_{-}=\delta_{1}-\delta_{2}$. The parameter $\delta_{+}$ is small and represents the detuning of the exact sum resonance: $\omega_{0}=\omega_{1}+\omega_{2}+\delta_{+}$ where $\delta_{+}=0$. Parameters $\mu_{1}$ and $\mu_{2}$ control the detuning of the damping from $\mu$; $\mu_{+}=\mu_{1}+\mu_{2}$, $\mu_{-}=\mu_{1}-\mu_{2}$.

The original nonlinear system that has the linearization (5) at zero detuning can be reduced to the following type of equation HR95

[TABLE]

where ${\bf A}_{0}$ is a $4\times 4$ matrix with the eigenvalues $\pm i\omega_{1},\pm i\omega_{2}$. The vector of parameters ${\bf p}=(\delta_{+},\delta_{-},\mu_{+},\mu_{-})$ is used to control detuning from resonance and damping.

The vector-valued function $\bf f$ is $2\pi$-periodic in $\omega_{0}t$ and ${\bf f}(0,\omega_{0}t;{\bf p})={\bf 0}$ for all $t$ and $\bf p$. Since the origin is a stationary point of (70), one may ask how its stability depends on the parameters. Analogous to the case of a single forced oscillator, one can make a planar stability diagram by varying the strength $\varepsilon$ and the frequency $\omega_{0}$ of the forcing while fixing the other parameters. Also in this case one obtains a resonance tongue at $\omega_{0}=\omega_{1}+\omega_{2}$. However if damping is varied, the planar stability diagram does not change continuously RTV93 . For instance, applying zero damping ($\mu=0$, no damping detuning) we find instability of the trivial solution (equilibrium) if

[TABLE]

For $\mu>0$ the trivial solution is unstable if

[TABLE]

The instability interval depends discontinuously on damping coefficient $\mu$!

Hoveijn and Ruijgrok HR95 presented a geometrical explanation of this dissipation-induced instability using ‘all’ the parameters as unfolding parameters first putting the equation (70) into a normal form Ar83 ; SVM which is similar to that of the non-semisimple $1:1$ resonance studied in F1992 .

In the normalized equation the time dependence appears only in the high order terms. But the autonomous part of this equation contains enough information to determine the stability regions of the origin.

The second step was to test the linear autonomous part ${\bf A}({\bf p})$ of the normalized equation for structural stability. This family of matrices is parameterized by the detunings of the frequencies $\omega_{1}$ and $\omega_{2}$ and of the damping parameter $\mu$.

Identifying the most degenerate member of this family one can show that ${\bf A}({\bf p})$ is its versal unfolding in the sense of Arnold Ar71 ; Ar83 . Put differently, the family ${\bf A}(\delta_{+},\delta_{-},\mu_{+},\mu_{-})$ is structurally stable, whereas ${\bf A}(\delta_{+},\delta_{-},0,0)$ is not. Therefore the stability diagram actually ‘lives’ in a four dimensional space. In this space, the stability regions of the origin are separated by a critical surface which is the hypersurface where ${\bf A}({\bf p})$ has at least one pair of imaginary complex conjugate eigenvalues. This critical surface is diffeomorphic to the Whitney umbrella, see Figure 14.

It is the singularity of the Whitney umbrella that causes the discontinuous behaviour of the planar stability diagram for the combination resonance in the presence of dissipation. The structural stability argument guarantees that the results are ‘universally valid’ and qualitatively hold for every system in sum resonance. For technical details we refer again to HR95 .

6 Abscissa minimization, robust stability and heavy damping

Let us return to the work of Bottema B56 . The conditions

[TABLE]

are necessary and sufficient for the polynomial

[TABLE]

to be Hurwitz. The domain (71) was plotted by Bottema in the $(a_{1},a_{3},a_{2})$-space, Figure 15a.

A part of the plane $a_{1}=a_{3}$ that lies inside the domain of asymptotic stability constitutes a set of all directions leading from the point (0, 0, 2) to the stability region

[TABLE]

The tangent cone (73) to the domain of asymptotic stability at the Whitney umbrella singularity, which is shown in green in Figure 15a,b, is degenerate in the $(a_{1},a_{3},a_{2})$ – space because it selects a set of measure zero on a unit sphere with the center at the singular point L80 ; L82 .

The singular point $(a_{1},a_{3},a_{2})=(0,0,2)$ corresponds to a double complex-conjugate pair of roots $\lambda=\pm i$ of the polynomial (72). The fact that multiple roots of a polynomial are sensitive to perturbation of the coefficients is a phenomenon that was studied already by Isaac Newton, who introduced the so-called Newton polygon to determine the leading terms of the perturbed roots as fractional powers of a perturbation parameter. It follows that, in matrix analysis, eigenvalues are in general not locally Lipschitz at points in matrix space with non-semi-simple eigenvalues Bu06 , and, in the context of dissipatively perturbed Hamiltonian systems, MO95 . Thus, it has been well-understood for a long time that perturbation of multiple roots or multiple eigenvalues on or near the stability boundary is likely to lead to instability KS05a .

Because of the sensitivity of multiple roots and eigenvalues to perturbation, in engineering and control-theoretical applications a natural desire is to “cut the singularities off” by constructing convex inner approximations to the domain of asymptotic stability. Nevertheless, multiple roots per se are not undesirable.

Indeed, multiple roots also occur deep inside the domain of asymptotic stability. Although it might seem paradoxical at first sight, such configurations are actually obtained by minimizing the so-called polynomial abscissa in an effort to make the asymptotic stability of a linear system more robust, as we now explain.

**Abscissa minimization and multiple roots

**The abscissa of a polynomial $p(\lambda)$ is the maximal real part of its roots BGMO2012 :

[TABLE]

We restrict our attention to monic polynomials with real coefficients and fixed degree $n$: since these have $n$ free coefficients, this space is isomorphic to $\mathbb{R}^{n}$. On this space, the abscissa is a continuous but non-smooth, in fact non-Lipschitz, as well as non-convex, function whose variational properties have been extensively studied using non-smooth variational analysis BGMO2012 ; Bu06 .

Now set $n=4$, consider the set of polynomials $p(\lambda)$ defined in (72), and consider the restricted set of coefficients

[TABLE]

On this set the roots are

[TABLE]

When $0\leq a_{1}<4$ $(a_{1}>4)$, the roots $\lambda_{1,2}$ and $\lambda_{3,4}$ are complex (real) with each pair being double, that is with multiplicity two. At $a_{1}=4$ there is a quadruple real eigenvalue $-1$. So, we refer to the set (75) as a set of exceptional points K2013dg ; KI2017 ; KI2018 (abbreviated as the EP-set).

When $a_{1}>0$, the EP-set (75) (shown by the red curve in Figure 15a,b) lies within the tangent cone (73) to the domain of asymptotic stability at the Whitney umbrella singularity $(0,0,2)$. The points in the EP-set all define polynomials with two double roots (denoted $EP_{2}$) except $(a_{1},a_{3},a_{2})=(4,4,6)$, at which $p$ has a quadruple root and is denoted $EP_{4}$; see Figure 15a,b.

Let us consider how the roots move in the complex plane when $a_{1}$ and $a_{3}$ coincide and are set to specific values and $a_{2}$ increases from zero, as shown by black curves in Figure 15c. When $a_{1}=a_{3}<4$, the roots that initially have positive real parts and thus correspond to unstable solutions move along the unit circle to the left, cross the imaginary axis at $a_{2}=2$ and merge with another complex conjugate pair of roots at $a_{2}=2+\frac{a_{1}^{2}}{4}$, i.e., at the EP-set. Further increase in $a_{2}$ leads to the splitting of the double eigenvalues, with one conjugate pair of roots moving back toward the imaginary axis. By also considering the case $a_{1}=a_{3}>4$, it is clear that when $a_{1}$ and $a_{3}$ coincide, the choice $a_{2}=2+\frac{a_{1}^{2}}{4}$ minimizes the abscissa, with the polynomial $p$ on the EP-set.

Furthermore, when $a_{1}=a_{3}$ is increased toward $4$ from below, the coalescence points ($EP_{2}$) move around the unit circle to the left. This conjugate pair of coalescence points merges into the quadruple real root $\lambda=-1$ ($EP_{4}$) when $a_{1}=a_{3}=4$ and hence $a_{2}=6$, as is visible in Figure 15c. If $a_{1}=a_{3}$ is increased above 4 the quadruple point $EP_{4}$ splits again into two exceptional points $EP_{2}$, one of them inside the unit circle.

Thus, all indications are that the abscissa is minimized by the parameters corresponding to $EP_{4}$, with a quadruple root at $-1$.

In fact, application of the following theorem shows that the abscissa of (72) is globally minimized by the $EP_{4}$ parameters.

Theorem 6.1

(BGMO2012 , Theorems 7 and 14)

Let $\mathbb{F}$ denote either the real field $\mathbb{R}$ or the complex field $\mathbb{C}$. Let $b_{0}$, $b_{1}$, $\ldots$, $b_{n}\in\mathbb{F}$ be given (with $b_{1},\ldots,b_{n}$ not all zero) and consider the following family of monic polynomials of degree $n$ subject to a single affine constraint on the coefficients:

[TABLE]

Define the optimization problem

[TABLE]

Let

[TABLE]

First, suppose $\mathbb{F}=\mathbb{R}$. Then, the optimization problem (76) has the infimal value

[TABLE]

where $h^{(i)}$ is the $i$-th derivative of $h$ and $k=\max\{j:~{}~{}b_{j}\neq 0\}$. Furthermore, the optimal value $a^{*}$ is attained by a minimizing polynomial $p^{*}$ if and only if $-a^{*}$ is a root of $h$ (as opposed to one of its derivatives), and in this case we can take

[TABLE]

Second, suppose $\mathbb{F}=\mathbb{C}$. Then, the optimization problem (76) always has an optimal solution of the form

[TABLE]

with $-\gamma$ given by a root of $h$ (not its derivatives) with largest real part.

In our case, $\mathbb{F}=\mathbb{R}$, $n=4$ and the affine constraint on the coefficients of $p$ is simply $a_{4}=1$. So, the polynomial $h$ is given by $h(\lambda)=\lambda^{4}-1$.

Its real root with largest real part is 1, and its derivatives have only the zero root. So, the infimum of the abscissa $a$ over the polynomials (72) is $-1$, and this is attained by

[TABLE]

that is, with the coefficients at the exceptional point $EP_{4}$. There is nothing special about $n=4$ here; if we replace $4$ by $n$ we find that the infimum is still $-1$ and is attained by

[TABLE]

Swallowtail singularity as the global minimizer of the abscissa

It is instructive to understand the set in the $(a_{1},a_{3},a_{2})$-space where the roots of the polynomial (72) are real and negative, but not necessarily simple, which is given by the discriminant surface of the polynomial. A part of it is shown in Figure 15a,b. At the point $EP_{4}$ with the coordinates $(4,4,6)$ in the $(a_{1},a_{3},a_{2})$- space the discriminant surface has the Swallowtail singularity, which is a generic singularity of bifurcation diagrams in three-parameter families of real matrices Ar71 ; Ar83 .

Therefore, the coefficients of the globally minimizing polynomial (77) are exactly at the Swallowtail singularity of the discriminant surface of the polynomial (72).

In the region in side the “spire” formed by the discriminant surface all the roots are simple real and negative. Owing to this property, this region, belonging to the domain of asymptotic stability (see Figure 15a), plays an important role in stability theory. Physical systems with semi-simple real and negative eigenvalues are called heavily damped. The solutions of the heavily damped systems do not oscillate and monotonically decrease, which is favorable for applications in robotics and automatic control.

Now we can give the following interpretation of the Bottema stability diagram shown in Figure 15a KO2013 . The dissipative system with the characteristic polynomial (72) is asymptotically stable inside the domain (71). The boundary of the domain (72) has the Whitney umbrella singular point at $a_{1}=0$, $a_{3}=0$, and $a_{2}=2$.

The domain corresponding to heavily damped systems is confined between three hypersurfaces of the discriminant surface and has a form of a trihedral spire with the Swallowtail singularity at its cusp at $a_{1}=4$, $a_{2}=6$, and $a_{3}=4$. The Whitney umbrella and the Swallowtail singular points are connected by the EP-set given by (75). At the Swallowtail singularity of the boundary of the domain of heavily damped systems, the abscissa of the characteristic polynomial of the damped system attains its global minimum.

Therefore, by minimizing the spectral abscissa one finds points at the boundary of the domain of heavily damped systems. Furthermore, the sharpest singularity at this boundary corresponding to a quadruple real eigenvalue $\lambda=-1$ with the Jordan block of order four is the very point where all the modes of the system with two degrees of freedom are decaying to zero as rapidly as possible when $t\rightarrow\infty$.