aa r X i v : . [ m a t h . A P ] N ov A catalogue of singularities
Jens Eggers ∗ and Marco A. Fontelos † ∗ School of Mathematics, University of Bristol, University Walk,Bristol BS8 1TW, United Kingdom † Instituto de Matem´aticas y F´ısica Fundamental,Consejo Superior de Investigaciones Cient´ıficasC/ Serrano 121, 28006 Madrid, Spain.
Abstract.
This paper is an attempt to classify finite-time singularities of PDEs.Most of the problems considered describe free-surface flows, which are easily observedexperimentally. We consider problems where the singularity occurs at a point, andwhere typical scales of the solution shrink to zero as the singularity is approached.Upon a similarity transformation, exact self-similar behaviour is mapped to the fixedpoint of a infinite dimensional dynamical system representing the original dynamics.We show that the dynamics close to the fixed point is a useful way classifying thestructure of the singularity. Specifically, we consider various types of stable andunstable fixed points, centre-manifold dynamics, limit cycles, and chaotic dynamics. ingularities Figure 1.
A drop of Glycerin dripping through PolyDimethylSiloxane (PDMS) nearsnap-off [1]. The nozzle diameter is 0 .
48 cm.
1. Introduction
A non-linear partial differential equation (PDE), starting from smooth initial data, willin general not remain smooth for all times. Consider for example the physical caseshown in Fig. 1, which we will treat in section 3 below. Shown is a snapshot of oneviscous fluid dripping into another fluid, close to the point where a drop of the innerfluid pinches off. This process is driven by surface tension, which tries to minimise thesurface area between the two fluids. At a particular point x , t in space and time,the local radius h ( z, t ) of the fluid neck goes to zero. This point is a singularity ofthe underlying equation of motion for the one-dimensional profile h ( z, t ). In particular,typical length scales near the pinch point go to zero (the minimum radius, at the veryleast, but also the solution’s axial extend). This absence of a characteristic scale nearthe singularity is the basic motivation to look for self-similar solutions.A fascinating aspect of the study of singularities is that they describe a greatvariety of phenomena which appear in the natural sciences and beyond [2]. Forexample, such singular events occur in free-surface flows [3], turbulence and Eulerdynamics (singularities of vortex tubes [4, 5] and sheets [6]), elasticity [7], Bose-Einstein condensates [8], non-linear wave physics [9], bacterial growth [10, 11], black-holecosmology [12, 13], and financial markets [14].In this paper we consider equations h t = F [ h ] , (1)where F [ h ] represents some (nonlinear) differential or integral operator. For simplicity,we will for the most part discuss the case of scalar h , but make reference to the many ingularities h is a vector, or (1) is a system of equations. Let ussuppose that (1) forms a localised singularity at x , t . If t ′ = t − t and x ′ = x − x , weare looking for local solutions of (1) which have the structure h ( x, t ) = t ′ α φ ( x ′ /t ′ β ) , (2)with appropriately chosen values of the exponents α, β .Giga and Kohn [15, 16] proposed to introduce self-similar variables τ = − ln( t ′ )and ξ = x ′ /t ′ β to study the asymptotics of blow up. Namely, putting h ( x, t ) = t ′ α H ( ξ, τ ) , (3)(1) is turned into the “dynamical system” H τ = G [ H ] ≡ αH + t ′ − α F [ t ′ α H ] . (4)If (2) is indeed a solution of (1), the right hand side of (4) is independent of τ , andself-similar solutions of the form (2) are fixed points of (4). By studying the “long-time”( τ → ∞ ) behaviour of solutions of (4) one can study the behaviour near blow-up. Forthe relation of singular PDE problems to the renormalisation group, developed in thecontext of critical phenomena, see [17, 18], for a more computational perspective, see[19]. Solutions to the original PDE (1) for given initial data can be viewed as orbits insome infinite dimensional phase phase, for instance, L . If the fixed point is an attractor,blow-up will be self-similar for some class of initial initial conditions. However, othertypes of attractors ( ω -limit sets in the notation which is customary in the context ofpartial differential equations, see [20] and references therein) are frequently observed,furnishing a very fruitful means of classifying singularities. In this paper we discuss thefollowing cases:(i) Stable fixed points
Initial conditions in some neighbourhood of the fixed point become attracted to it,so the solution converges exponentially to the self-similar solution (2). Formally,two eigenvalues of the linearisation around the fixed point are always positive, butthese unstable motions can be absorbed into a redefinition of the origin of space andtime; this is discussed in section 2.1. A sub-classification into self-similarity of thefirst and of the second kind is due to Barenblatt [21]. Self-similar solutions are ofthe first kind if (2) only solves (1) for one set of exponents α, β ; their values are fixedby either dimensional analysis or symmetry, and are thus rational. Solutions are ofthe second kind (in the sense of Barenblatt) if solutions (2) exist for a continuousset of exponents α, β ; the exponents are fixed by a non-linear eigenvalue problem,and take irrational values in general.(ii)
Unstable fixed points
The linearisation around the fixed point possesses positive eigenvalues, so it is neverreached, except for non-generic initial data. Often a stable fixed point is associatedwith an infinite sequence of unstable fixed points. ingularities
Travelling waves
Solutions of (1) converge to h = t ′ α φ ( ξ + cτ ), which is a travelling wave solution of(4) with propagation velocity c .(iv) Centre manifold .This is also known, when it leads to singularities that develop at a faster rate thanthe selfsimilar scaling, as type-II self-similarity [22]; it arises if one eigenvalue of thelinearisation around the fixed point is zero, and there is a non-linear dependenceinstead. This typically leads to corrections involving logarithmic time τ . We discusstwo different cases, involving quadratic and cubic non-linearities.(v) Limit cycles
This is also known as “discrete self-similarity” [12, 23]. Corresponding solutionshave the form h = t ′ α ψ [ ξ, τ ] with ψ being a periodic function of period T in τ .Thus at the discrete sequence of times τ n = τ + n T , which approaches the singulartime for n → ∞ , the solution looks like a simple self-similar one.(vi) Strange attractors
In principle, more complex behaviour is possible, where the orbits of the dynamicalsystem lie on a strange attractor. At the moment, we are not aware of an equationexhibiting such behaviour which would correspond to any physical phenomenon.However, this may simply be due to the fact that the corresponding singularbehaviour is more difficult to detect numerically, and to grasp analytically. Todemonstrate that such behaviour is at least possible, we show that any finitedimensional dynamical system may be “embedded” in the singular dynamics. Asan explicit example, we show that the phase-space trajectory may lie on the Lorenzattractor.(vii)
Exotic objects
There might be other types of behaviour that have no analogue in finite-dimensionaldynamical systems. In particular, blow-up may occur at several points ( x , t atthe same time, in which case the description (4) is not so useful.This paper’s aim is to assemble the body of knowledge on singularities of equationsof the type (1) that is available in both the mathematical and the applied community,and to categorise it according to the types given above. In addition to rigorous resultswe pay particular attention to various phenomenological aspects of singularities whichare often crucial for their appearance in an experiment or a numerical simulation. Forexample, what are the implications of the type of singularity for the approach of thePDE solution onto the self-similar form (2)? In most cases, we rely on known examplesfrom the mathematical physics literature. To find an explicit example for limit cyclebehaviour as well as chaotic dynamics, we propose a new set of model equations, inspiredby a problem in general relativity. ingularities Figure 2.
Scanning electron microscopy images illustrating the pinch-off of a row ofrectangular troughs in silicone (top) [24]. The bottom picture shows the same sampleafter 10 minutes of annealing at 1100 ◦ C. The troughs have pinched off to form a rowof almost spherical voids. The dynamics is driven by surface diffusion.
2. Stable and unstable fixed points
Our example, exhibiting self-similarity of the first kind (in the sense of Barenblatt) [21],is that of a solid surface evolving under the action of surface diffusion. Namely, atomsmigrate along the surface driven by gradients of chemical potential, see Fig.2. Theresulting equations in the axisymmetric case, where the free surface is described by thelocal neck radius h ( x, t ), are [25]: h t = 1 h (cid:20) h (1 + h x ) / κ x (cid:21) x , (5)where κ = 1 h (1 + h x ) / − h xx (1 + h x ) / (6)is the mean curvature. In (5),(6), all lengths have been made dimensionless using anouter length scale R (such as the initial neck radius), and the time scale R /D , where D is a forth-order diffusion constant.At a time t ′ ≪ ℓ = t ′ / isa local length scale. This suggests the similarity form h ( x, t ) = t ′ / H ( x ′ /t ′ / ) , (7) ingularities Figure 3.
The approach to the self-similar profile for equation (5). The dashed lineis the stable similarity solution H ( ξ ) as found from (8). The full lines are rescaledprofiles found from the original dynamics (5) at h min = 10 − , − , and h min = 10 − ,respectively. As the singularity is approached, they converge rapidly onto the similaritysolution (7). implying α = β = 1 /
4. This is the classical situation for self-similarity of the first kind,more examples are found in [21, 26]. The similarity form of the PDE becomes −
14 ( H − ξH ′ ) = 1 H " H (1 + H ′ ) / κ ′ ′ , ξ = x ′ t ′ / (8)where the prime denotes differentiation with respect to ξ .Solutions of (8) have been studied extensively in [27]. To ensure matching to a time-independent outer solution, the leading order time dependence must drop out from (7),implying that H ( ξ ) ∼ c | ξ | , ξ → ±∞ . (9)It turns out that all similarity solutions are symmetric, so only one constant c needs tobe determined. Exactly as in the closely related problem of surface-tension-driven fluidpinch-off, the requirement of a certain growth condition (9) is enough to fix a uniquesolution of (8) [28]. The value of c comes out as part of the solution. Solutions of (8)with the growth condition (9) form a discretely infinite set [27], again like the fluidsproblem [29]. The series of similarity solutions is conveniently ordered by descendingvalues of the minimum, see table 1.Next we turn to the dynamical system that describes the dynamics away from thefixed point, by putting h ( x, t ) = t ′ / H ( ξ, τ ) , (10) ingularities H i (0) c i Table 1.
A series of similarity solutions of (8) as given in [27]. The higher-ordersolutions become successively thinner and flatter. where τ = − ln( t ′ ). The similarity form of (5) becomes H τ = 14 ( H − ξH ξ ) + 1 H " H (1 + H ξ ) / κ ξ ξ , (11)which reduces to (8) if the left hand side is set to zero. To assure matching of (11) tothe outer solution, we require the boundary condition H τ − ( H − ξH ξ ) / → | ξ | → ∞ . (12)Next we linearise (11) around H = H ( ξ ), by writing H = H ( ξ ) + ǫP ( ξ, τ ), whichgives P τ = L ( H ) P. (13)Since H satisfies (12), P ( ξ, τ ) must do the same. In particular, this means that if L ( H ) P i = ν i P i , (14)i.e. if ν i is an eigenvalue of the linear operator, the corresponding eigenfunction mustgrow like P i ( ξ ) ∝ ξ − ν i . (15)If the similarity solution H ( ξ ) is to be stable, the eigenvalues of L ( H ) must benegative. However, there are always two positive eigenvalues, which are related to theinvariance of the equation of motion (5) under translations in space and time. Namely,for any ǫ , the translated similarity solution h ǫ ( x, t ) = H ( x ′ + ǫℓ ) (16)is an equally good self-similar solution of (5), and thus of (11). In particular, we canexpand (16) to lowest order in ǫ , and find that H ǫ ( ξ, τ ) = H ( ξ ) + ǫe βτ H ′ ( ξ ) + O ( ǫ ) (17)is a solution of (13).Thus, since L H = 0, ǫe βτ βH ′ = ∂H ǫ ( ξ, τ ) ∂τ = L H ǫ = ǫe βτ L H ′ . (18) ingularities ν x = β ≡ / L with eigenfunction H ′ ( ξ ).Similarly, considering the transformation t → t + ǫ , one finds a second positive eigenvalue ν t = 1, with eigenfunction ξH ′ . To reiterate, the physical meaning of these eigenvalues isthat upon perturbing the similarity solution, the singularity time as well as the positionof the singularity will change. Thus if the coordinate system is not adjusted accordingly,it looks as if the solution would flow away from the fixed point. If, on the other hand,the solution is represented relative to the perturbed values of x and t , the dynamicswill converge onto a stable similarity solution.The eigenvalues of the solutions H i have been found numerically in [27]. Theresult is that the linearisation around the “ground state” solution H only has negativeeigenvalues (apart from the two trivial ones), while all the other solutions have at leastone other positive eigenvalue. This means that H is the only similarity solution thatcan be observed, all other solutions are unstable. Close to the fixed point, the approachto H will be dominated by the largest negative eigenvalue ν : h ( x, t ) = t ′ / (cid:2) H ( ξ ) + ǫt ′ ν P ( ξ ) (cid:3) . (19)For large arguments, the point ξ cr where the correction becomes comparable to thesimilarity solution is ξ ∼ ǫt ′ ν ξ − ν , and thus ξ cr ∼ t ′ / . This means that the region ofvalidity of H ( ξ ) expands in similarity variables, and is constant in real space. This rapidconvergence is reflected by the numerical results reported in Fig. 3. More formally, onecan say that for any ǫ there is a δ such that (cid:12)(cid:12) h ( x, t ) − t ′ / φ ( ξ ) (cid:12)(cid:12) ≤ ǫ (20)if | x ′ | ≤ δ uniformly as t ′ → In the example of the previous subsection, the exponents can be determined bydimensional analysis, and therefore assume rational values. As Barenblatt [21] pointsout, there are problems where the scaling behaviour depends on external parameters,set for example by the initial conditions. In that case, the scaling exponent can assumeany value. Often, this value is fixed by some intrinsic property of the equation, resultingin an irrational answer. We will call this situation self-similarity of the second kind (inthe sense of Barenblatt). A particularly simple example of this kind of singularity isthe pinch-off of a very viscous thread of liquid [30, 3], which we present now. Anotherrecent example is the pinch-off of a two-dimensional inviscid sheet [31].For simplicity, we confine ourselves to the case of a slender viscous filament withoutinertia, for which the equation becomes: h t ( s, t ) = 16 (cid:18) C ( t ) h ( s, t ) (cid:19) . (21)The typical velocity scale γ/η , where γ is the surface tension and η is the viscosity,has been absorbed into the time variable. The particularly simple form of (21) hasbeen achieved by writing the thread radius in Lagrangian variables , i.e. as function of a ingularities Figure 4.
A drop of viscous fluid falling from a pipette 1 mm in diameter [33]. Notethe long neck. particle label s . This means the particle is at position z ( s, t ) at time t , and z t ( s, t ) is thevelocity at time t . The time-dependent constant of integration C ( t ) must be determinedfrom the constraint that u ( s, t ) ≡ z s = 1 /h ( s, t ).Note that the self-similar form (2) is a solution of (21) for α = 1, and any value of β . Thus contrary to the example described in the previous section, the exponents arenot completely determined from dimensional analysis or from balancing powers of t ′ inthe equation of motion. This is a typical situation in which self-similarity of the secondkind is observed [32]. Instead, the unknown exponent is determined from a non-lineareigenvalue equation, and takes an irrational value.Since α = 1 we introduce u = t ′− f ( ξ ) , with ξ = s/t ′ β (22)and C ( t ) = Kt ′ . (23)Hence 1 √ f + 3 (cid:18) f + βξf ′ f (cid:19) = K, (24)where K is an arbitrary constant. Imposing symmetry and regularity of f , we introducean expansion of f ( ξ ) of the form f ( ξ ) = R − + Cξ n + O ( ξ n ) , n = 1 , , ... into (24) to obtain the condition R = 112( nβ −
1) (25)and define β = nβ. Equation (24) can easily be integrated in terms of ln ξ and y = √ f : Z dy ((1 + 6 R ) y − y − R y ) = 16 R β ln ξ + e C = 16 R β ln ξ n + e C, with e C an arbitrary constant. Computing the integral above we obtain y − β (cid:0)(cid:0) β − (cid:1) y + 1 (cid:1) β − (1 − y ) = Cξ n , (26) ingularities C = 1. In terms of the profile h ( s, t ) thiscorresponds to fixing the scale of the spatial variable s . Solutions are undetermined upto such a scale factor, as is clear from the invariance of (21) under a change in spatialscale. As a result, the axial scale is fixed by the initial conditions.The value of the velocity C n at infinity is therefore given by C n = Z ∞ z st ds = Z ∞ u t ds = 13 Z ∞ (cid:18)(cid:18)
124 2 β − β − (cid:19) f − f (cid:19) dξ. From the condition that the velocity at infinity must vanish we thus obtain C n = 0,which will be the equation that determines the exponent β . Taking the derivative of(26) we obtain nξ n − dξdy = ddy (cid:16) y − β (cid:0)(cid:0) β − (cid:1) y + 1 (cid:1) β − (1 − y ) (cid:17) == − y − β − (cid:0) yβ − y + 1 (cid:1) β − β p (1 − y )and hence K n ( β ) ≡ β − C n = βn Z (cid:18)(cid:18)
12 2 β − β − (cid:19) y − y (cid:19)(cid:18) y − n + βn (cid:0)(cid:0) β − (cid:1) y + 1 (cid:1) −
12 2 n − β +1 n (1 − y ) −
12 2 n − n (cid:19) dy. (27)The function K n ( β ) may be written explicitly as K n ( β ) = β Γ (4 − β ) Γ (cid:0) n (cid:1) Γ (cid:0) − β + n (cid:1) nβ − nβ − F (cid:18) n + 12 n − β, − β ; 4 − β + 12 n ; 1 − nβ (cid:19) −− β Γ (3 − β ) Γ (cid:0) n (cid:1) Γ (cid:0) − β + n (cid:1) F (cid:18) n + 12 n − β, − β ; 3 − β + 12 n ; 1 − nβ (cid:19) (28)with roots given by table 2. If one converts the Lagrangian variables back to the originalspatial variables, one obtains h ( x, t ) = t ′ φ St (cid:0) x ′ /t ′ β − (cid:1) . (29)Thus for t ′ → t ′ of the generic n = 1 solution rapidly becomessmaller than the axial scale t ′ . (cf. table 2). This explains the long necks seen inFig. 4.For generic initial data u ( s ) = h − ( s ) = B + B s + B s + ... + B j s j + ... oneexpects that B = 0, so that the self-similar solution with n = 1 will develop. Onlyif B i = 0 for i = 1 , , ..., n − B n = 0 the n -th self-similar solution will be theasymptotic description of the solution. For this reason, only the n = 1 solution isstable, since a generic perturbation of the initial data with B i = 0 for i = 1 , , ..., n − B = 0. ingularities β Table 2.
A list of exponents, found from K n ( β ) = 0, with K n given by (28).The number 2 n gives the smallest non-vanishing power in a series expansion of thecorresponding similarity solution around the origin. Only the solution with n = 1 isstable. The minimum radius is found from (25).
3. Travelling wave
The pinching of a liquid thread in the presence of an external fluid is described by theStokes equation [34]. For simplicity, we consider the case that the viscosity η of the fluidin the drop and that of the external fluid are the same. An experimental photographof this situation is shown in Fig. 1. To further simplify the problem, we make theassumption (the full problem is completely analogous) that the fluid thread is slender.Then the equations given in [1] simplify to h t = − v z h/ − vh z , (30)where v = 14 Z z + z − h z ′ ( z ′ )( h ( z ′ ) + ( z − z ′ ) ) / dz ′ . (31)Here we have written the velocity in units of the capillary speed v η = γ/η . The limitsof integration z − and z + are for example the positions of the plates which hold a liquidbridge [35].Dimensionally, one would once more expect a local solution of the form h ( z, t ) = t ′ H out (cid:18) z ′ t ′ (cid:19) , (32)and H out ( ξ ) has to be a linear function at infinity to match to a time-independent outersolution. In similarity variables, (31) has the form V out ( ξ ) = 14 Z z b/t ′ − z b /t ′ H ′ out ( ξ ′ ) p H + ( ξ − ξ ′ ) dξ ′ . (33)We have chosen z b as a real-space variable close to the pinch-point, such that thesimilarity description is valid in [ − z b , z b ]. But if H out is linear, the integral in (33)diverges, which means that a simple “fixed point” solution (32) is impossible.However, the integral can be made convergent by introducing a shift in the similarityvariable ξ : h = t ′ H out ( ξ − bτ ) , (34) ingularities τ = − ln( t ′ ) as usual. This means in similarity variables the solution is a travellingwave. With this modification, the mass balance (30) becomes − H out + H ′ out ( ξ + V out + bτ ) = H out V ′ out / . (35)Now we choose b such that the logarithmic singularity cancels, namely we demand that14 Z z b /t ′ − z b /t ′ H ′ out ( ξ ′ ) p H out + ( ξ − ξ ′ ) dξ ′ − b ln( t ′ ) (36)finite for t ′ →
0. This is achieved by putting b = 14 " − H + p H + 1 + H − p H − + 1 . (37)Thus defining V fin ( ξ ) = lim Λ →∞ Z Λ − Λ H ′ out p H + ( ξ − ξ ′ ) dξ ′ + b ln Λ (38)the similarity equation − H out + H ′ out ( ξ + V fin − ξ ) = H out V ′ fin / ξ is an arbitrary constant. It remains as an arbitrary axial shift in thesimilarity solution.The numerical solution of the integro-differential equation (39) gives h min = a out v η t ′ , where a out = 0 . . (40)The slope of the solution away from the pinch-point are given by H + = 4 .
81 and H − = − . , (41)which means the solution is very asymmetric, as confirmed directly from Fig. 1.
4. Centre manifold
In section 2 we described the generic situation that the behaviour of a similarity solutionis determined by the linearisation around it. In the case of a stable fixed point,convergence is fast, and the observed behaviour is essentially that of the fixed point.In this section, we describe two different cases where the largest eigenvalue vanishes, sohigher-order non-linear terms have to be taken into account. The approach to the fixedpoint is now much slower, and much more of the observed behaviour is determined bythe approach to the fixed point. Depending on the type of non-linearity, there is greatfreedom of possible behaviours, of which we discuss two. ingularities Figure 5.
Nine images (of width 3.5 mm) showing how a He crystal “flows” downfrom the upper part of a cryogenic cell into its lower part [37]. The recording takes afew minutes, the temperature is 0.32 K. 11 mK. The crystal first “drips” down, so thata crystalline “drop” forms at the bottom (a to c); then a second drop appears (d) andcomes into contact with the first one (e); coalescence is observed (f) and subsequentlybreakup occurs (h).
Axisymmetric motion by mean curvature in three spatial dimensions is described by theequation h t = (cid:18) h xx h x − h (cid:19) , (42)where h ( x, t ) is the radius of the moving free surface. A very good physical realizationof (42) is the melting and freezing of a He crystal, driven by surface tension [36], seeFig. 5. As before, the time scale t has been chosen such that the diffusion constant,which sets the rate of motion, is normalised to one. A possible boundary condition forthe problem is that h (0 , t ) = h ( L, t ) = R , where R is some prescribed radius. For certaininitial conditions h ( x, ≡ h ( x ) the interface will become singular at some time t , atwhich h ( x , t ) = 0 and the curvature blows up. The moment of blow-up is shown inpanel h of Fig. 5, for example.Inserting the self-similar solution (2) into (42), one finds a balance for α = β = 1 / − φ ξ φ ′ (cid:18) φ ′′ φ ′ − φ (cid:19) , ξ = x ′ t ′ / . (43)One solution of (43) is the constant solution φ ( ξ ) = √
2. Another potential solutionis one that grows linearly at infinity, to ensure matching onto a time-independentouter solution. However, it can be shown that no solution to (43), which also growslinearly at infinity, exists [38, 39]. Our analysis below follows the rigorous work in [22],demonstrating type-II self-similarity. In addition, we now show how the description ofthe dynamical system can be carried out to arbitrary order. ingularities h ( x, t ) = t ′ / h √ g ( ξ, τ ) i , (44)with τ = − ln( t ′ ) as usual. The equation for g is then g τ = g − ξg ′ g ′′ g ′ − g / p g ′ , (45)which we solve by expanding into eigenfunctions of the linear part of the operator L g = g − ξg ′ / g ′′ . (46)It is easily confirmed that L H i ( ξ/
2) = (1 − i ) H i ( ξ/ , (47)where H n is the n-th Hermite polynomial [40]: H n = ( − n e x d n dx n e − x . (48)Thus the linear part of (45) becomes ∂a i ∂τ = (1 − i ) a i , (49)which means that all eigenvalues are negative except for the first, which vanishes. Toinvestigate the approach of the cylindrical solution, one must therefore include nonlinearterms in the equation for a .If we write g ( ξ, τ ) = ∞ X i =1 a i ( τ ) H i ( ξ/ , (50)the equation for a becomes da dτ = − / a + O ( a a j ) , (51)whose solution is a = 1 / (2 / τ ) . (52)Thus instead of the expected exponential convergence onto the fixed point, theapproach is only algebraic. Since all other eigenvalues are negative, the τ -dependenceof the a i is effectively determined by a . Namely, as we will see below, a j = O ( τ − j ), socorrections to (51) are of higher order.If one linearises around (52), putting a = a (0)1 + ǫ , one finds dǫ dτ = − τ ǫ + other terms . (53)This means that the coefficient A of ǫ = A/τ remains undetermined, and a simpleexpansion of a i in powers of τ − yields an indeterminate system. Instead, at quadraticorder, a term of the form ǫ = A ln τ /τ is needed. Fortunately, this is the only place in ingularities a i where such an indeterminacy occurs. Thus alllogarithmic dependencies can be traced, leading to the general ansatz a ( n ) i = δ i τ i + n X k = i +1 k − i X l =0 (ln τ ) l τ k δ lki , (54)where δ i and δ lki are coefficients to be determined. The index n is the order of thetruncation.Indeed, the coefficients can be found recursively by considering terms of successivelyhigher order in τ − in the first equation: da dτ = − / a − √ a a + 22 a − √ a − √ a + 192 a a (55) da dτ = − a − √ / a + 6 a − √ a a (56)The next two orders will involve the next coefficient a . From (55) and (56), onefirst finds δ and δ , by considering O ( τ − ) and O ( τ − ), respectively. Then, at order O ( τ − ( n +1) ) in the first equation, where n = 3, one finds all remaining coefficients δ lki in the expansion (54) up to k = n . At each order in τ − , there is of course a seriesexpansion in ln τ which determines all the coefficients.We constructed a MAPLE program to compute all the coefficients up to arbitrarilyhigh order (10th, say). Up to third order in τ − the result is: a = 1 / √ τ + 1716 ln ( τ ) √ τ − √ τ +867128 ln ( τ ) √ τ − τ )) √ τ (57) a = − / √ τ + 516 √ τ − τ ) √ τ , (58)and thus h ( x, t ) becomes h ( x, t ) = t ′ / h √ a ( τ ) (cid:0) − ξ (cid:1) + a ( τ ) (cid:0) − ξ + ξ (cid:1)i , (59)from which one finds the minimum. To second order, the result is h min = (2 t ′ ) / (cid:20) − τ − τ τ (cid:21) . (60)Two remarks are in order. First, the presence of logarithms implies that there issome dependence on initial conditions built into the description. The reason is that theargument inside the logarithm needs to be non-dimensionalised using some “external”time scale. More formally, any change in time scale ˜ t = t/t leads to an identical equationif also lengths are rescaled according to ˜ h = h/ √ t . This leaves the prefactor in (60)invariant, but adds an arbitrary constant τ to τ . This is illustrated by comparing to anumerical simulation of the mean curvature equation (42) close to the point of breakup,see Fig. 6. Namely, we subtract the analytical result (60) from the numerical solution h min / (2 √ t ′ ) and multiply by τ . As seen in Fig.6, the remainder is varying slowly over 12 ingularities τ -14 -12 -10 -8 -6 -4-6-5.5-5-4.5-4-3.5 Figure 6.
A plot of h h min / √ t ′ − / (2 τ ) i τ (dashed line) and τ / − (3+17 ln( τ + τ ) /
8) (full line) with τ = 4 . decades in t ′ . If the constant τ is adjusted, this small variation is seen to be consistentwith the logarithmic dependence predicted by (60).The second important point is that convergence in space is no longer uniform asimplied by (20) for the case of self-similarity of the first kind. Namely, to leadingorder the pinching solution is a cylinder. For this to be a good approximation, onehas to require that the correction is small: ξ /τ ≪
1. Thus corrections becomeimportant beyond ξ cr ∼ τ , which, in view of the logarithmic growth of τ , impliesconvergence in a constant region in similarity variables only . As shown in [36], theslow convergence toward the self-similar behaviour has important consequences for acomparison to experimental data. The next example is that of bubble breakup [41], for which a very different form ofnonlinearity is observed. As shown in [41], the equation for a slender cavity or bubbleis Z L − L ¨ a ( ξ, t ) dξ p ( z − ξ ) + a ( z, t ) = ˙ a a , (61)where a ( z, t ) ≡ h ( z, t ). The integral runs over the fluid domain. If for the momentone disregards boundary conditions looks for solutions to (61) of cylindrical form, a ( z, t ) = a ( t ), one can do the integral to find¨ a ln (cid:18) L a (cid:19) = ˙ a a . (62) ingularities Figure 7.
The pinch-off of an air bubble in water [44]. An initially smooth shapedevelops a localised pinch-point. −14 −12 −10 −8 −6 −4 −2 t ⁄ α α PSfrag replacements t ′ Figure 8.
A comparison of the exponent α between full numerical simulations ofbubble pinch-off (solid line) and the leading order asymptotic theory (69) (dashedline). It is easy to show that an an asymptotic solution of (62) is given by a ∝ ∆ t ln(∆ t ) / , (63)corresponding to a power law with a small logarithmic correction. Indeed, initial theoriesof bubble pinch-off [42, 43] treated the case of an approximately cylindrical cavity, whichleads to the radial exponent α = 1 /
2, with logarithmic corrections.However both experiment [44] and simulation [41] show that the cylindrical solutionis unstable; rather, the pinch region is rather localised, see Fig. 7. Therefore, it is notenough to treat the width of the cavity as a constant L ; the width ∆ is itself a time-dependent quantity. In [41] we show that to leading order the time evolution of theintegral equation (61) can be reduced to a set of ordinary differential equations for theminimum a of a ( z, t ), as well as its curvature a ′′ .Namely, the integral in (61) is dominated by a local contribution from the pinchregion. To estimate this contribution, it is sufficient to expand the profile around theminimum at z = 0: a ( z, t ) = a + a ′′ / z + O ( z ). As in previous theories, the integraldepends logarithmically on a , but the axial length scale is provided by the inversecurvature ∆ ≡ (2 a /a ′′ ) / . Thus evaluating (61) at the minimum, one obtains [41] toleading order ¨ a ln(4∆ /a ) = ˙ a / (2 a ) , (64) ingularities a and ∆. Thus, a second equation is needed to close thesystem, which is obtained by evaluating the the second derivative of (61) at the pinchpoint: ¨ a ′′ ln (cid:18) e a ′′ (cid:19) − a a ′′ a = ˙ a ˙ a ′′ a − ˙ a a ′′ a . (65)The two coupled equations (64),(65) are most easily recast in terms of the time-dependent exponents2 α ≡ − ∂ τ a /a , δ ≡ − ∂ τ a ′′ /a ′′ , (66)where τ ≡ − ln t ′ and β = α − δ , are a generalisation of the usual exponents α and β .The exponent δ characterises the time dependence of the aspect ratio ∆. Returning tothe collapse (62) predicted for a constant solution, one finds that α = 1 / δ = 0.In the spirit of the the previous subsection, this is the fixed point corresponding tothe cylindrical solution. Now we expand the values of α and δ around their expectedasymptotic values 1 / α = 1 / u ( τ ) , δ = v ( τ ) . (67)To leading order in δ , the resulting equations are ∂ τ u = − vu , ∂ τ v = − v , (68)which describe perturbations around the leading-order similarity solution. Theseequations are analogous to (51), but they have a degeneracy of third order, ratherthan second order. Equations (68) are easily solved to yield, in an expansion for small δ [41], α = 1 / √ τ + O ( τ ) , δ = 14 √ τ + O ( τ − / ) . (69)Thus the exponents converge toward their asymptotic values α = β = 1 / α ≈ . − .
58 [44], and why there is a weak dependence on initialconditions [45].The cubic equation (67) applies to the exponents α , δ , rather than the solution itself,as in the previous subsection, where we dealt with mean curvature flow. In fact, it iseasy to re-analyse the solution to the mean curvature problem, and to formulate it interms of time-dependent exponents. Let us define a and a ′′ for the mean curvatureproblem through h ( z, t ) = a + a ′′ / z + O ( z ). From (44) it follows that a = 2 t ′ and a ′′ = 2( √ g ) g ′′ . This gives a ′′ = 4 √ a , and thus ∂ τ δ = − δ , (70)instead of (68). This means mean curvature flow, formulated in terms of exponents, oncemore gives a quadratic non-linearity, rather than the cubic term (68) found for bubblecollapse. From (70) one finds to leading order δ mc = 1 / (2 τ ) and α mc = 1 / O (1 /τ )for the mean curvature flow. ingularities
5. Limit cycles
An example for this kind of blow-up was introduced into the literature in [12] in thecontext of cosmology. There is considerable numerical evidence [46] that discrete self-similarity occurs at the mass threshold for the formation of a black hole. The sametype of self-similarity has also been proposed for singularities of the Euler equation [47]and for a variety of other phenomena [48]. A reformulation of the original cosmologicalproblem leads to the following system: f x = ( a − fx , (71)( a − ) x = 1 − (1 + U + V ) /a x , (72)( a − ) t = (cid:20) ( f + x ) U − ( f − x ) V x + 1 (cid:21) /a − , (73) U x = f [(1 − a ) U + V ] − xU t x ( f + x ) , (74) V x = f [(1 − a ) U + V ] + xV t x ( f − x ) . (75)In [13], the self-similar description corresponding to the system (71)-(75) was solvedusing formal asymptotics and numerical shooting procedures. This leads to the solutionsobserved in [12]. Below we propose a very similar system, which we solve analytically,and which shows the same type of limit cycle behaviour: u t + u x = 2 f v, (76) v t + v x = − f u, (77) f t = f . (78)The simplified system (76)-(78) can be solved introducing characteristics: η = x + t, ν = x − t, (79)which leads to u η = f v, v η = − f u. (80)The system (80) has to be solved at constant ν with initial conditions u ( ν, ν ) = u init ( ν ) , v ( ν, ν ) = v init ( ν ) , (81)where u init ( x ) ≡ u ( x,
0) and v init ( x ) ≡ v ( x,
0) are the profiles at t = 0.The solution of (78) is f ( x, t ) = 11 /f init ( x ) − t , (82) ingularities f init ( x ) is the initial profile, which can be written in terms of ν, η as f ( η, ν ) = 11 /f init (( η + ν ) / − ( η − ν ) / . (83)The function f can be eliminated from (80) using the transformation ζ = Z ην f ( η ′ , ν ) dη ′ , (84)which results in the system for a harmonic oscillator: u ζ = v, v ζ = − u. (85)The solution is u = C ( ν ) sin( ζ + φ ( ν )) , v = C ( ν ) cos( ζ + φ ( ν )) , (86)with C ( ν ) = p u init ( ν ) + v init ( ν ) and φ ( ν ) = arcsin( u init ( ν ) /C ( ν )).According to (82), a singularity of the PDE system (76)-(78) first occurs at the maximum f of f init ( x ), which we assume to occur at x = 0 without loss of generality.Thus locally we can write f ( x, ≈ f − ax , and for small t ′ = t − t we have f = t ′− at ξ , ξ = xt ′ / , (87)for any finite ξ . Using this explicit form of f , (84) can be integrated to find ζ . Atconstant ξ we have ν = − t + ξt ′ / + t ′ , η = t + ξt ′ / − t ′ , (88)so taking the limit t ′ →
0, the leading order result for (84) is ζ = − ln (cid:20) t ′ (1 + at ξ )(1 + at ) t (cid:21) ≡ τ + φ ( ξ ) . (89)Thus as the singularity is reached, t ′ →
0, the variable ζ goes to infinity. This meansthe singularity corresponds to the long-time limit of the dynamical system (85), whichwill perform a harmonic motion according to (86). Namely, for t ′ → u = q u init ( − t ) + v init ( − t ) sin ( arcsin u init ( − t ) p u init ( − t ) + v init ( − t ) ! − ln (cid:20) t ′ (1 + at ξ )(1 + at ) t (cid:21)(cid:27) . (90)Thus the singular solution is of the general form u = ψ ( φ ( ξ ) + τ ) , (91)where ψ is periodic in τ . This is a particularly simple version of discretely self-similarbehaviour. Note that the character of the solution is different from the travelling wavesdescribed in section 3. ingularities
6. Strange attractors and exotic behaviour
In connection to limit cycles and in the context of singularities in relativity, a fewinteresting situations have been found numerically quite recently. One of them is theexistence of Hopf bifurcations where a self-similar solution (a stable fixed point) istransformed into a discrete self-similar solution (limit cycle) as a certain parametervaries (see [49]). Other kinds of bifurcations, for example of the Shilnikov type, arefound as well [50]. Now we demonstrate that chaotic behaviour is also possible.In section 4.1 we treated a system of an infinite number of ordinary differentialequations for the coefficients of the expansion of an arbitrary perturbation to an explicitsolution. Such high-dimensional systems in principle allow for a rich variety of dynamicalbehaviours, including those found in classical finite dimensional dynamical systems, suchas chaos. Consider for instance an equation for the perturbation g (the analogue of (45))of the form g τ = L g + F ( g, g ) , (92)where L g is a linear operator. Assuming an appropriate non-linear structure for thefunction F , an arbitrary nonlinear (chaotic) dynamics can be added.To give an explicit example of a system of PDEs exhibiting chaotic dynamics,consider the structure of the example given in the previous section 5. It can begeneralised to produce any low-dimensional dynamics near the singularity. Namely,let us generalise the system (76)-(78) to u ( i ) t + u ( i ) x = 2 f F i ( { u ( i ) } ) , i = 1 , . . . , n, (93) f t = f . (94)Using the transformation (79), the first n equations are turned into: u ( i ) η /f = F i , i = 1 , . . . , n, (95)which is an ODE system for constant ν . The system (95) has to be solved with initialconditions u ( i ) ( ν, ν ) = u ( i ) init ( ν ) , i = 1 , . . . , n. (96)As before, the function f can be eliminated using (84), resulting in the generalnon-linear dynamical system u ( i ) ζ = F i , i = 1 , . . . , n. (97)Now if one chooses n = 3 and F = σ ( u (2) − u (1) ) , F = ρu (1) − u (2) − u (1) u (3) , F = u (1) u (2) − βu (3) , (98)(97) is the Lorenz system [51].As before, for t ′ →
0, the variable ζ goes to infinity, and near the singularity oneis exploring the long-time behaviour of the dynamical system (97). In the case of (98), ingularities ρ , the resulting dynamics will be chaotic. Specifically, taking σ = 10, ρ = 28, and β = 8 /
3, as done by Lorenz [52], the maximal Lyapunov exponent is0 . U (1) ( ν, ζ ) , U (2) ( ν, ζ ) , U (3) ( ν, ζ ) is a solution of (97) with initial conditions(96), the final form of the similarity solution is u ( i ) ( x, t ) = U ( i ) ( − t + ξt ′ / , τ + φ ( ξ )) , i = 1 , , . (99)In the limit t ′ → − t + ξt ′ / ≈ − t . Thereason is that in the limit of large τ two trajectories diverge like exp( λ max τ ) = t ′− λ max ,where the largest Lyapunov exponent is larger than 1 /
7. Outlook
The singularities described in this paper are point-like in the sense that they occur ata single point x at a given time t . There are two important situations we have notdiscussed and for which the dynamical systems point of view and analytical approachdoes not apply. First, the case of singularities not developing at single points, but onsets of finite measure. This is the case in a few simple examples of reaction-diffusionequations of the family u t − ∆ u = u p − b |∇ u | q for x ∈ Ω (100)where depending on values of p > q > u may be regional ( u blows up in subsets of Ω of finite measure) or even global (the solutionblows-up in the whole domain). See for instance [53] and references therein. Anotherinteresting possibility is that the stable fixed point which is approached depends on theinitial conditions. Thus exponents could vary either discretely or continuously with thechoice of initial conditions. The infinite sequence of similarity solutions found for (21)is not an example for such behaviour, since only one solution is stable. All other fixedpoints will not be visible in practise.Singularities may even happen in sets of fractional Hausdorff dimension, i.e.,fractals. This is the case of the inviscid one-dimensional system for jet breakup (cf. [54])and might be case of Navier-Stokes system in three dimensions, where the dimension ofthe singular set at the time of first blow-up is at most 1 (cf. [55]). This connects with thesecond issue we did not address here. It is the nature of the singular sets both in spaceand time. In many instances, existence of global in time (for all 0 ≤ t < ∞ ) solutions tononlinear problems can be established in a weak sense , that is, allowing certain kind ofsingularities to develop both in space and time. In the case of 3-D Navier-Stokes system,the impossibility of singularities ”moving” in time, that is of curves x = ϕ ( t ) in thesingular set is well-known [55]. Hence, provided a certain singularity does not persist intime, the question is how to continue the solutions after a singularity has developed. ingularities Acknowledgments
This paper is an outgrowth of discussions between the authors and R. Deegan, preparinga workshop on singularities at the Isaac Newton Institute, Cambridge. We thank J. M.Martin-Garcia and J. J. L. Velazquez for fruitful discussions and for providing us withvaluable references.
References [1] Cohen I, Brenner M P, Eggers J and Nagel S R 1999
Phys. Rev. Lett. Phys. Today
Rev. Mod. Phys. J. Fluid Mech.
Phys. Rev. Lett. PNAS
Phys. Rev. Letters Phys. Lett. A
Phys. Rev. Lett. J. Math. Biol. Nonlinearity Phys. Rev. Lett. Phys. Rev. D Phys. Rep.
Comm. Pure Appl. Math. Indiana University Math. J. Lectures on phase transitions and the renormalization group (Addison-Wesley)[18] Bricmont J, Kupiainen A and Lin G 1994
Comm. Pure Appl. Math. J. Non-Newtonian Fluid Mech.
A Stability Technique for Evolution Partial DifferentialEquations: A Dynamical Systems Approach (Birkhauser)[21] Barenblatt G I 1996
Similarity Self-Similarity and Intermedeate Asymptotics (Cambridge)[22] Angenent S B and Vel´azquez J J L 1997
J. reine angew. Math.
Living Rev. Rel., to be published [24] Mizushima I, Sato T, Taniguchi S and Tsunashima Y 2000
Appl. Phys. Lett. J. Appl. Phys. ZAMM J. Stat. Phys. Phys. Rev. Lett. Phys. Fluids Phys. Fluids Phys. Fluids ?? ???[32] Bensimon D, Kadanoff L P, Liang S, Shraiman B I and Tang C 1986 Rev Mod Phys New J. Phys. art. no. 59[34] Lister J R and Stone H A 1998 Phys. Fluids Acad. Sci. Bruxelles Mem. Phys. Rev. E Phys. Rev. Lett. J. Geom. Anal. ingularities [39] Huisken G 1993 Proc. of Symposia in Pure Math. Handbook of Mathematical Functions (Dover)[41] Eggers J, Fontelos M A, Leppinen D and Snoeijer J H 2007
Phys. Rev. Lett. J. Fluid Mech.
J. Fluid Mech.
Phys. Fluids Phys.Rev. Lett. Phys. Rep.
Physica D
Phys. Rep.
Phys. Rev. D Class. Quant. Grav. S299–S306[51] Strogatz S H 1994
Nonlinear Systems and Chaos (Perseus publishing)[52] Lorenz E N 1963
J. Atmos. Sci. Electron. J. Diff. Eqns.
European J. Appl. Math. Comm. Pur. Appl. Math.35