Gradient catastrophe and flutter in vortex filament dynamics
aa r X i v : . [ n li n . S I] J un Gradient catastrophe and flutter in vortex filament dynamics
B.G.Konopelchenko and G.Ortenzi ∗ October 30, 2018
Abstract
Gradient catastrophe and flutter instability in the motion of vortex filament within the localizedinduction approximation are analyzed. It is shown that the origin of this phenomenon is in thegradient catastrophe for the dispersionless Da Rios system which describes motion of filament withslow varying curvature and torsion. Geometrically this catastrophe manifests as a rapid oscillation ofa filament curve in a point that resembles the flutter of airfoils. Analytically it is the elliptic umbilicsingularity in the terminology of the catastrophe theory. It is demonstrated that its double scalingregularization is governed by the Painlev´e-I equation.
PACS:47.32.C, 02.30.Ik, 47.35.JkMSC: 76B47, 58K35, 35Q05Keywords: Vortex Filament Dynamics, Gradient catastrophe, Flutter
Motion of a thin vortex filament in an incompressible inviscid fluid in an infinite three-dimensional domainwithin the localized induction approximation (LIA) is governed by the simple equation for the inducedvelocity ~v ~X t ( s, t ) ≡ ~v ( s, t ) = ~X s ∧ ~X ss = K ( s, t ) ~b (1)which implies the following intrinsic equations for the curvature K and torsion τK t = − K s τ − Kτ s ,τ t = KK s − τ τ s + (cid:18) K ss K (cid:19) s . (2)Here ~X ( s, t ) is a position vector for a point on the curve representing the filament, t is time, s is thearclenght parameter, subscripts indicate the differentiation with respect to the indicated variables, and ~b is the binormal.Equations (1) and (2) have been derived by L. S. Da Rios in 1906 [1] and have been rediscovered sixtyyears later in [2, 3]. The detailed history of Da Rios (DR) system is presented in [4] For the localizedinduction approximation see e.g. [5, 6].In 1972 Hasimoto has observed [7] that equation (1) is equivalent to the focusing nonlinear Schr¨odingerequation (NLS) iψ t + ψ ss + 12 | ψ | ψ = 0 (3)via the transformation ψ ( s, t ) = K ( s, t ) exp (cid:18) i Z s τ ( s ′ , t ) ds ′ (cid:19) . (4)Only one year before Zakharov and Shabat [8] discovered that NLS equation (3) is integrable by theinverse scattering transform (IST) method. Hasimoto’s result has demonstrated that the whole powerful ∗ Corresponding author. E-mail: [email protected], Phone +39(0)264485765 Fax:+39(0)264485705. β = − τ + iK and slow variables x = ǫs , y = ǫt , ǫ ≪ β y = 12 (cid:0) β + β (cid:1) β x . (5)This system well approximates DR system (2) until the derivatives β x , β y are not large. Hodographequations for equation (5) describe critical points of the function W which obey the Euler-Poisson-Darboux equation E (cid:0) , (cid:1) .Gradient catastrophe for the dDR system happens in a point ( x , y ) for given initial data. At thispoint β x and β y (or K x , K y , τ x , τ y ) explode. The acceleration ~a = ~v t explodes too. So the filamentbecomes fast oscillating near the point x at the “time ” y . Numerical analysis (borrowed from [33])shows that such oscillations begin to expand along the filament. Intrinsically, the filament begins tooscillate around the rectifying plane near the point x . So, the gradient catastrophe for the system (5)gives rise to fast oscillation which can be referred as filament flutter by analogy with airfoil flutter (seee.g. [39]).Analysis of behavior near the point x , y of gradient catastrophe shows that β = β + ǫβ ∗ , x = x + ǫx ∗ , y = y + ǫy ∗ , ǫ ≪
1, and W = W + ǫ (cid:18) x ∗ ( f U + gV ) + 13 U − U V (cid:19) (6)where U + iV ∝ β and f , g are some constants. Thus the function W exhibits an elliptic umbilicsingularity behavior in the terminology of Thom [40]. Regularization of this singularity is one of theissues of this paper.At the point of gradient catastrophe the approximation (5) becomes invalid and should be substitutedby the full DR system (2). First order corrections can be obtained by the double-scaling technique togetherwith appropriate modification of the function W ( β, β ) to a functional W ∗ such that the equations ofcritical points for W are substituted by the Euler-Lagrange equation for W ∗ . The resulting equation forsmall corrections is equivalent to the Painlev´e-I (P-I) equationΩ ξξ = 6Ω − ξ. (7)The result of the paper [36] allows us to conjecture that any generic solution of the gradient catastrophebehaves as β = β + ǫβ ∗ where the correction β ∗ is described by the tritronque´e solution of the P-Iequation.The paper is organized as follows. Dispersionless DR system, its different forms, hodograph equationand other its properties are discussed in section 2. Gradient catastrophe for the dDR system is analysedin section 3. Geometrical implications of gradient catastrophe and flutter of filament are consideredin section 4. In section it is shown that near to the singular point filament exhibits the elliptic umbiliccatastrophe behavior. Regularization of this singularity via the Painlev`e-I equation is presented in section6. 2 Dispersionless Da Rios system
In order to describe dispersionless (or quasiclassical) limit of DR system one, in a standard manner,introduces slow variables x = ǫs , y = ǫt with ǫ ≪ K = K ( x, y ) and τ = τ ( x, y ). Under these assumptions equation (2) take theform K y = − K x τ − Kτ x ,τ y = KK x − τ τ x + ǫ (cid:18) K xx K (cid:19) x . (8)Thus in the limit ǫ → K y = − K x τ − Kτ x ,τ y = KK x − τ τ x . (9)Analogously to the dispersionless NLS equation [34]–[37] the solutions of the dDR system (9) well approx-imate solutions of the DR system in points were K x and τ x are finite. As any two-component system it islinearizable by hodograph transformation ( x, y ) ↔ ( K, τ ). The characteristic velocities for dDR systemare complex λ = − τ + iK , Riemann invariant β = − τ + iK and in terms of β the dDR system is of theform (9). We note that solutions of the dNLS system discussed in [36] and those of dDR system (9) areconnected by the simple relations u = K , v = 2 τ . For the geometrical consideration the system (9) ismore convenient.The dDR system is known to be integrable similar to its hyperbolic version, i.e. the 1-layer Benneysystem [41]. It has an infinite set of symmetries and integrals of motion. One of the form of the dDRhierarchy is given by the set of equations (see e.g. [42]) p y n = (cid:18) z n p (cid:19) + p ! x , n = 1 , , , . . . (10)where p = z + p ( x, y ) z − + p ( x, y ) z − + p ( x, y ) z − + . . . is a formal Laurent series defined by theequations p = ( z − β )( z − β ) = ( z + τ ) + K , (11) y n are time variables and f + denotes the polynomial part of f . The dDR system (9) is the first flowof the hierarchy (10) at n = 1. All p k +1 ( x, y ) are densities of integrals of motion for the dDR system.They are the dispersionless limit of the densities of integrals of motion for the full DR system (2) foundin [14, 17].Solutions of the dDR system can be calculated via the standard hodograph equation. It was shownin [43, 44] that these hodograph equations are, in fact, the equations W β = 0 , W β = 0 , (12)which define the critical points of the function of the form W = x β + β ) + y β + 2 ββ + 3 β ) + f W ( β, β ) (13)where the function f W ( β, β ) = f W ( β, β ) is defined by the initial data for β and it is such that W obeysthe Euler-Poisson-Darboux equation E (cid:0) , (cid:1) , i.e.2( β − β ) W ββ = W β − W β . (14)Since the function W is real valued, the second equation of (12) is the complex conjugated to the firstone. Hodograph equation is W β = x y β + β ) + f W β = 0 . (15)In terms of K and τ the function W is W = − xτ + y (cid:18) τ − K (cid:19) + f W ( τ, K ) , (16)3he hodograph equations (15) are given by W K = − yK + f W K = 0 W τ = − x + 2 yτ + f W τ = 0 , (17)while the Euler-Poisson-Darboux equation is of the form K ( W KK + W ττ ) + W K = 0 . (18)It is the axisymmetric three-dimensional Laplace equation studied by Beltrami [45] (see also [46]). Here K and τ play the role of the radial and axial coordinate respectively in the cylindrical system of coordinates. For the dDR hierarchy (10) the function W has the form W = ∞ X l =0 y l I γ dz πi z l +1 p ( z ) (19)where γ denotes a small circle around infinity. Euler-Poisson-Darboux equations (14) is quite useful inthe study of singular sector of the 1-layer Benney and dDR hierarchies [43, 44] Hodograph equations (12) provide us with unique solution of the dDR system if the standard condition∆ ≡ det (cid:18) W ββ W ββ W ββ W ββ (cid:19) = 0 (20)is satisfied. Due to the Euler-Poisson-Darboux equations (14) on the solutions of the dDR system with K = 0 the function W obeys the equation W ββ = 0. Hence∆ = | W ββ | (21)and, consequently, the hodograph equations (15) are uniquely solvable if W ββ = 0. In this situation thederivatives β x and β x are bounded together with β . Indeed, differentiating equation (15) with respect to x , one gets the system (cid:18) W ββ W ββ W ββ W ββ (cid:19) (cid:18) β x β x (cid:19) + 12 (cid:18) (cid:19) = 0 . (22)So for the solutions of the dDR system one has β x = −
12 1 W ββ . (23)Thus regular sector of the dDR system (9) is characterized by the condition [44] W β = 0 , W ββ = 0 . (24)In this sector β and the derivatives β x , β y , i.e. curvature K , torsion τ , and their derivatives are bounded.Singular sector of dDR system is composed by solutions for which W ββ = 0. For these solutions β (i.e. curvature and torsion) remain bounded while their derivatives explode. Such situation is usuallyreferred as gradient catastrophe (see e.g. [47, 48]). Generic gradient catastrophe is characterized by theconditions [44] W β = 0 , W ββ = 0 , W βββ = 0 . (25)At K = 0 ( β = β ) one also has W ββ = 0 and the conditions (25), in terms of curvature and torsion are W K = 0 , W τ = 0 , W KK = 0 , W Kτ = 0 , W KKK = 0 , W τττ = 0 . (26)4quations W ττ = 0, W KKτ = 0, W Kττ = 0, are consequences of (26).Last conditions (26) allow us to solve the conditions W KK = W ττ = 0 with respect to K and τ .Substituting these expressions into the first two conditions (26), one gets two equations f ( x, y ) = 0 , f ( x, y ) = 0 . (27)Thus generic gradient catastrophe for the dDR system happens in a single point ( x , y ) for given W ( β, β ),i.e. initial data for K and τ [44]. For the focusing dNLS equation this fact has been first noted in [36].We would like to emphasize that in contrast to the hyperbolic case where the catastrophe locus is a curvein x, y here the catastrophe happens in a point. Gradient catastrophe for the dDR system means a special behavior of a filament at the point x at time y . First one observes, using formula (20) in [16] with U = W = 0 and V = K , that the acceleration ~a for the velocity (1) ~a = ~v t = ǫ~v y = ǫ ( K y ~b + K~b y ) = ǫ h ( − K x τ − Kτ x ) ~b − K x ~t + τ ~n i (28)explodes at x . Thus, at the moment t = ǫy sharp bump is formed at the point s = ǫx of the filament.Local consideration reveals another peculiarity of behavior of filament around the point x . It is wellknown, that the coordinates of a curve in a neighborhood of a point in the reference system formed bythe tangent vector ~t , normal ~n , and binormal ~b and origin in the point (see e.g. formula (6.7) [49]) are x = s − K s + . . . ,x = K s − K s s + . . . ,x = Kτ s + 124 (2 K s τ + Kτ s ) s + . . . . (29)Passing to the slow variable x and slow coordinates ˜ x = ǫx , ˜ x = ǫx , ˜ x = ǫx , one gets the sameseries (29) for slow variables with substitution s → x . At ordinary point where K , τ and K x , τ x arebounded the first terms in r.h.s. of (29) dominate and the curve, locally, is a twisted cubic. Its projectionon the osculating plane (spanned by ~t and ~n ), normal plane (spanned by ~b and ~n ) and rectifying plane(spanned by ~t and ~b ) are parabola ( x − K x = 0), cusp (9 Kx − τ x = 0) and cubic ( x + Kτ x = 0)respectively (fig. 1, see e.g. [49]). Characteristic feature of a curve (and a filament) in an ordinary point xx –1012345 –3 –2 –1 1 2 3 xx –4–224 –1 1 2 3 4 5 xx –4–224 –3 –2 –1 1 2 3 Figure 1: The figure show a curve with τ >
0. At τ < ~t in the last figure.(i.e., outside the points of gradient catastrophe) is that it lies always on one side of positive direction ofthe normal and on one side of osculating plane (depending on sign of torsion τ ).5 xx –4–224 –3 –2 –1 1 2 3 Figure 2: Possible different forms of the projection on the osculating plane of a vortex filament near tothe catastrophe point.At the point x of gradient catastrophe, the behaviors of a curve changes drastically. Indeed, when K x and τ x become large, the second terms in ˜ x and ˜ x become relevant. So, parabola in the osculatingplane may change sign or even convert into cubic curve (fig. 2) while in the normal plane it could be aplane (3 ,
4) curve and so on. So, around the point x of gradient catastrophe a filament oscillates fromone side of the rectifying planes to another and back. Such oscillation is quite similar to that of airfoil(see e.g. [39]) and one can refer to these oscillation of a filament in the point of gradient catastrophe asa flutter.At each point of a curve there is a sphere which has contact of a third order at this point. It is calledthe osculating sphere (see e.g. [49]). It has radius R = K τ + K s K τ . (30)and its center has coordinates ~X o = ~X + 1 K ~n − K s K τ ~b (31)where ~X is the position vector of a point on the curve.At the point of gradient catastrophe the radius of the osculating sphere and its center blow up toinfinity.Oscillation of a curve around the rectifying plane and blow up of the osculating sphere are geometricalfeatures of the gradient catastrophe and flutter of a filament. At the point x , t of gradient catastrophederivatives K x , τ x are so large that the terms of higher order derivatives in the DR equation (8) becomerelevant for finite ǫ . Solutions of the full DR system (8) at small ǫ are almost indistinguishable from thoseof the dDR system at t < t = ǫy (for dNLS equation see [31]–[37]). At t ≥ t the behavior of solutionsof dDR and DR system are completely different. Solutions of the DR system at t ≥ t develop a zoneof rapid oscillations which expand around the point x . For the equivalent dNLS equation this fact hasbeen observed both analytically and numerically (see e.g. [31]–[37]).An example of such a behaviors is given by the figure 13 of paper [33]. In the terms of K and τ the initial data from [33] are K ( x,
0) = 2 e − x and τ ( x,
0) = − tanh( x ). Gradient catastrophe happensin x = 0 and y ≃ .
25. At this point the corresponding solution of DR system becomes oscillate. Astime y growth a zone of rapid oscillations expands as √ y − y symmetrically around the point x . Forthe filament dynamics this effect is seen as the creation of rapid oscillations at the point of gradientcatastrophe and their subsequent expansion along the growing in time piece of filament (see fig. 3). Sucha behavior is an evidence of the flutter type instability of the vortex filament motion.6igure 3: Typical behavior of a flutter vortex line. In order to understand better this phenomenon one should analyse in more detail the structure of thesingularity and the regularizing mechanism for the gradient catastrophe. For the focusing dNLS equationsuch an analysis based on the ǫ -expansion of the integrals of motions of the NLS/Toda equation has beenperformed in [36]. Here we will follow a different approach discussed recently in [50].Thus we consider a neighborhood of the gradient catastrophe point x , y and denote the values of β at this point by β . Following to the double scaling limit method (see e.g. [51]–[52]) we will look forsolutions of the DR system near the point of gradient catastrophe of the form x = x + ǫ α x ∗ y = y + ǫ σ y ∗ β = β + ǫ γ β ∗ (32)where ǫ ≪ α, σ, γ should be fixed by further consideration. For the sake of simplicity werestrict ourselves to the case y ∗ = 0. We first consider the function W ( β, β ) (13). One has W ( x + ǫ α x ∗ , y , β + ǫ α β ∗ ) = W + 12 ( β + β ) x ∗ + ǫ γ (cid:20)(cid:18) x + 14 (3 β + β ) + f W β (cid:19) β ∗ + (cid:18) x + 14 ( β + 3 β ) + f W β (cid:19) β ∗ (cid:21) + 12 ǫ α + γ x ∗ ( β ∗ + 3 β ∗ )+ 12 ǫ γ (cid:20)(cid:18) y + f W ββ (cid:19) β ∗ + (cid:18) y + f W ββ (cid:19) β ∗ + (cid:18)
12 + 2 f W ββ (cid:19) β ∗ β ∗ (cid:21) + 16 ǫ γ hf W βββ β ∗ + 2 f W βββ β ∗ β ∗ + 2 f W βββ β ∗ β ∗ + f W βββ β ∗ i + . . . (33)where W = W (cid:12)(cid:12)(cid:12) β = β , f W β = ∂ f W∂β (cid:12)(cid:12)(cid:12) β = β , f W β = ∂ f W∂β (cid:12)(cid:12)(cid:12) β = β and so on.Hodograph equation (15) and conditions (25) imply that x y β + β ) + f W β = 0 , y + f W ββ = 0 . (34)At curvature K = 0 the Euler-Poisson-Darboux equations (14) and its differential consequences2 βW ββ + 2( β − β ) W βββ = W ββ − W ββ , (35)imply that W ββ = 0 , W βββ = W βββ = 0 , (36)i.e. y f W ββ = 0 , f W βββ = f W βββ = 0 . (37)7aking into account equations (34) and (37), one finally gets W ( x + ǫ α x ∗ , y , β + ǫ γ β ∗ ) = W + 12 ǫ α ( β + β ) x ∗ + 12 ǫ α + γ ( β ∗ + β ∗ ) x ∗ + 13 ǫ γ ( aβ ∗ + aβ ∗ ) + . . . (38)where a = f W βββ . Consequently near the point x , y the hodograph equations take the form W β ∗ = 12 ǫ α + γ x ∗ + ǫ γ aβ ∗ = 0 . (39)These equations readily imply that α = 2 γ and β ∗ = x ∗ / . So, near to the point x β ∗ x ∗ ∼ ( x − x ) − / . (40)Thus near to the point of gradient catastrophe the function W is of the form W = W + 12 ǫ γ ( β + β ) x ∗ + ǫ γ W ∗ (41)where W ∗ = 12 x ∗ ( β ∗ + β ∗ ) + 13 ( aβ ∗ + aβ ∗ ) (42)and β = − τ + iK . Denoting a / β ∗ = U + iV , one gets W ∗ = x ∗ ( f U + gV ) + 13 U − U V (43)where f and g are certain constants. In the Thom’s catastrophe theory a function of this form is knownas describing elliptic umbilic singularity (e.g. [40, 53]). So gradient catastrophe and flutter for vortexfilament motion enter into the general scheme of Catastrophe Theory. For the focusing NLS equation theappearance of elliptic umbilic singularity has been observed for the first time in [36]. Next step is to regularize this singularity. Within the double scaling limit method adopted also in [36]one begins with the full dispersive system. Performing appropriate double scaling limit one calculatethe required regularizing terms. An approach discussed in [50] suggests to proceed in opposite direction,namely, to begin with the original (say dispersionless system), to modify the corresponding function W adding to it the differential terms of lowest order with appropriate scaling and then to require that thecritical point equations (12) for W are substituted by the Euler Lagrange equations for the modified W q ,i.e. δW q δβ = 0 , δW q δβ = 0 . (44)Following this approach we first observe that the contributions of β ∗ and β ∗ into the perturbation W ∗ of W are separated. It is quite natural to modify the function (33) in such a way that this separation andreality property W ( β, β ) = W ( β, β ) are preserved. Thus, omitting inessential second term in (41), andchoosing without lost of generality γ = 1, one has the following modified WW q = W + ǫ (cid:20) x ∗ ( β ∗ + β ∗ ) + 13 ( aβ ∗ + aβ ∗ ) (cid:21) + 12 ǫ δ h bβ ∗ x ∗ + b β ∗ x ∗ i (45)where δ and b are appropriate constants. The corresponding Euler-Lagrange equation for β ∗ and β ∗ is ∂W q ∂β ∗ − (cid:18) ∂W q ∂β ∗ x ∗ (cid:19) x ∗ = 0 , (46)i.e. ǫ (cid:20) x ∗ β ∗ + aβ ∗ (cid:21) − ǫ δ [ bβ ∗ x ∗ x ∗ ] = 0 . (47)8o balance the first term one should choose δ = 3.Thus, the modified W q is of the form W q = W + ǫ (cid:20) x ∗ ( β ∗ + β ∗ ) + 13 ( aβ ∗ + 12 aβ ∗ ) + 12 bβ ∗ x ∗ + 12 b β ∗ x ∗ (cid:21) (48)and the corresponding equation for β ∗ is bβ ∗ x ∗ x ∗ = aβ ∗ + 12 x ∗ . (49)This equation is converted into the classical Painlev´e-I equation (7) by simple change of variables x ∗ = λξ, β ∗ = − λ b Ω (50)with aλ = − b .So, the regularized behavior at the point of the gradient catastrophe of DR system, i.e. for thevortex filament dynamics is governed by the Painlev´e-I equation. The relevance of this equation forthe NLS/Toda system near to the critical point has been observed in a different way in [36]. Thecorrespondence between the solution of equation (49) and those given by the formula (5.21) in [36] is (at b = 1) u = K , v = 2 τ and x = − iε / (cid:0) re iψ (cid:1) − / (cid:18) − a K (cid:19) / x ∗ ,β = − ibε / (cid:0) re iψ (cid:1) (cid:18) − a K (cid:19) / β ∗ (51)where x , β , r and ψ are defined in (2.19) and (2.22) of [36]. One concludes also that ε = ǫ / . In thepaper ([36]) it was conjectured that any generic solution of he NLS equation near to the critical pointis described by the tritronque´e solution of the Painlev´e-I equation. Extension of this observation to ourcase allows to formulate the Conjecture 6.1.
Generic behavior of the vortex filament near the point x , y of gradient catastrophe isgiven by β ( s, t , ξ ) ≃ β − ǫ (cid:18) − a b (cid:19) / Ω (cid:18)(cid:16) − a b (cid:17) / s − s ǫ (cid:19) (52)where Ω is the tritronque´e solution of the Painlev´e-I equation (7).Properties of the solution (52) as well as other features of flutter for vortex filaments will be discussedelsewhere. Acknowledgements
This work has been partially supported by PRIN grant no 28002K9KXZ and byFAR 2009 (
Sistemi dinamici Integrabili e Interazioni fra campi e particelle ) of the University of MilanoBicocca.
References [1] Da Rios L.S., Rend. Circ. Mat. Palermo An introduction to fluid mechanics , Cambridge Univ. Press, Cambridge (1967)[6] Saffman P.G.,
Vortex dynamics , Cambridge Univ. Press, New York (1992)97] Hasimoto H., J. Fluid. Mech.
118 (1971); Sov. Phys. JTEP
62 (1972)[9] Lamb Jr. G.L., J. Math. Phys.
408 (1983)[12] Cieslinski J., Gragert P.K.H., and Sym A., Phys. Rev. Lett.
323 (1991); II, Physica D
267 (1991)[16] Nakayma K, Segur H. and Wadati M., Phys. Rev. Lett. , 2603-2606 (1992)[17] Ricca R., Phys. Fliuds A
31 (1999)[21] Ruban V. P., Phys. Rev. E, Geometry and modulation theory for the periodic nonlinear Schroedingerequation in: IMA Volumes in Mathematics and its applications The nonlinear Schrdinger equation as both aPDE and a dynamical system
Handbook of dynamical systems Semiclassical soliton ensembles for the focusingnonlinear Schrdinger equation in: Annals of Mathematics Studies
Princeton University PressPrinceton NJ (2003)[35] Tovbis A., Venarides S. and Zhou X., Commun. Pure Appl. Math. Theoretical Hydrodynamics
McMillan & co. (1960)[40] Thom R.,
Structural stability and morphogenesis. An outline of a general theory of models
Addison-Wesley Publishing Company Redwood City CA (1989)[41] Zakharov V.E., Func. Anal. Appl.
89 (1980)[42] Kodama Y. and Konopelchenko B., J. Phys. A: Math. Gen. Systems of quasilinear equations and their applications togas dynamics
Transl. Math. Mon. AMS Providence RI (1983)[48] Whitham G.B.,
Linear and nonlinear waves
John Wiley and Sons NY (1974)[49] Eisenhart L.P.,
An introduction to differential geometry
Princeton Univ. Press princeton (1947)[50] Konopelchenko B. and Mart´ınez Alonso L.,
Catastrophes, double scaling regularization and integrablesystems (in preparation)[51] Di Francesco P., Ginsparg P. and Zinn-Justin Z., Phys. Rep. bf 254 2 (1995)[52] Mart´ınez Alonso L. and Medina E., J. Phys. A: Math. Theor. Singularities of differentiable maps. Vol. I. Theclassification of critical points, caustics and wave fronts
Monographs in Mathematics82