{'fraction_non_alphanumeric': 0.04622778910573842, 'fraction_numerical': 0.015671424410136907, 'mean_word_length': 4.219857077694997, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A catalogue of singularities

This paper is an attempt to classify finite-time singularities of PDEs. Most of the problems considered describe free-surface flows, which are easily observed experimentally. We consider problems where the singularity occurs at a point, and where 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 fixed point of a {\it infinite dimensional dynamical system} representing the original dynamics. We show that the dynamics close to the fixed point is a useful way classifying the structure of the singularity. Specifically, we consider various types of stable and unstable fixed points, centre-manifold dynamics, limit cycles, and chaotic dynamics.


Introduction
A non-linear partial differential equation (PDE), starting from smooth initial data, will in general not remain smooth for all times. Consider for example the physical case shown in Fig. 1, which we will treat in section 3 below. Shown is a snapshot of one viscous fluid dripping into another fluid, close to the point where a drop of the inner fluid pinches off. This process is driven by surface tension, which tries to minimise the surface area between the two fluids. At a particular point x 0 , t 0 in space and time, the local radius h(z, t) of the fluid neck goes to zero. This point is a singularity of the 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 very least, but also the solution's axial extend). This absence of a characteristic scale near the singularity is the basic motivation to look for self-similar solutions.
In this paper we consider equations 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 important cases where the variable h is a vector, or (1) is a system of equations. Let us suppose that (1) forms a localised singularity at x 0 , t 0 . If t ′ = t 0 − t and x ′ = x − x 0 , we are looking for local solutions of (1) which have the structure 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(ξ, τ ), (1) is turned into the "dynamical system" If (2) is indeed a solution of (1), the right hand side of (4) is independent of τ , and self-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. For the relation of singular PDE problems to the renormalisation group, developed in the context 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 in some infinite dimensional phase phase, for instance, L 2 . If the fixed point is an attractor, blow-up will be self-similar for some class of initial initial conditions. However, other types of attractors (ω-limit sets in the notation which is customary in the context of partial differential equations, see [20] and references therein) are frequently observed, furnishing a very fruitful means of classifying singularities. In this paper we discuss the following 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, but these unstable motions can be absorbed into a redefinition of the origin of space and time; this is discussed in section 2.1. A sub-classification into self-similarity of the first and of the second kind is due to Barenblatt [21]. Self-similar solutions are of the first kind if (2) only solves (1) for one set of exponents α, β; their values are fixed by either dimensional analysis or symmetry, and are thus rational. Solutions are of the second kind (in the sense of Barenblatt) if solutions (2) exist for a continuous set 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 never reached, except for non-generic initial data. Often a stable fixed point is associated with an infinite sequence of unstable fixed points.
(iv) Centre manifold. This is also known, when it leads to singularities that develop at a faster rate than the selfsimilar scaling, as type-II self-similarity [22]; it arises if one eigenvalue of the linearisation around the fixed point is zero, and there is a non-linear dependence instead. This typically leads to corrections involving logarithmic time τ . We discuss two different cases, involving quadratic and cubic non-linearities.
(v) Limit cycles This is also known as "discrete self-similarity" [12,23]. Corresponding solutions have the form h = t ′ α ψ [ξ, τ ] with ψ being a periodic function of period T in τ . Thus at the discrete sequence of times τ n = τ 0 + nT, which approaches the singular time 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 dynamical system lie on a strange attractor. At the moment, we are not aware of an equation exhibiting such behaviour which would correspond to any physical phenomenon. However, this may simply be due to the fact that the corresponding singular behaviour is more difficult to detect numerically, and to grasp analytically. To demonstrate that such behaviour is at least possible, we show that any finite dimensional dynamical system may be "embedded" in the singular dynamics. As an explicit example, we show that the phase-space trajectory may lie on the Lorenz attractor.
(vii) Exotic objects There might be other types of behaviour that have no analogue in finite-dimensional dynamical systems. In particular, blow-up may occur at several points (x 0 , t 0 at the 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 equations of 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 results we pay particular attention to various phenomenological aspects of singularities which are often crucial for their appearance in an experiment or a numerical simulation. For example, what are the implications of the type of singularity for the approach of the PDE solution onto the self-similar form (2)? In most cases, we rely on known examples from the mathematical physics literature. To find an explicit example for limit cycle behaviour as well as chaotic dynamics, we propose a new set of model equations, inspired by a problem in general relativity.  [24]. The bottom picture shows the same sample after 10 minutes of annealing at 1100 • C. The troughs have pinched off to form a row of almost spherical voids. The dynamics is driven by surface diffusion.

Stable and unstable fixed points
2.1. Self-similarity of the first kind 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, atoms migrate along the surface driven by gradients of chemical potential, see Fig.2. The resulting equations in the axisymmetric case, where the free surface is described by the local neck radius h(x, t), are [25]: where is the mean curvature. In (5), (6), all lengths have been made dimensionless using an outer length scale R (such as the initial neck radius), and the time scale R 4 /D 4 , where D 4 is a forth-order diffusion constant. At a time t ′ ≪ 1 away from breakup, dimensional analysis implies that ℓ = t ′1/4 is a local length scale. This suggests the similarity form The approach to the self-similar profile for equation (5). The dashed line is the stable similarity solution H(ξ) as found from (8). The full lines are rescaled profiles found from the original dynamics (5) at h min = 10 −1 , 10 −2 , and h min = 10 −3 , respectively. As the singularity is approached, they converge rapidly onto the similarity solution (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 where the prime denotes differentiation with respect to ξ. Solutions of (8) have been studied extensively in [27]. To ensure matching to a timeindependent outer solution, the leading order time dependence must drop out from (7), implying that It turns out that all similarity solutions are symmetric, so only one constant c needs to be determined. Exactly as in the closely related problem of surface-tension-driven fluid pinch-off, the requirement of a certain growth condition (9) is enough to fix a unique solution 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 fluids problem [29]. The series of similarity solutions is conveniently ordered by descending values of the minimum, see table 1.
Next we turn to the dynamical system that describes the dynamics away from the fixed point, by putting  Table 1. A series of similarity solutions of (8) as given in [27]. The higher-order solutions become successively thinner and flatter.
where τ = − ln(t ′ ). The similarity form of (5) becomes which reduces to (8) if the left hand side is set to zero. To assure matching of (11) to the outer solution, we require the boundary condition Next we linearise (11) around H = H(ξ), by writing H = H(ξ) + ǫP (ξ, τ ), which gives Since H satisfies (12), P (ξ, τ ) must do the same. In particular, this means that if i.e. if ν i is an eigenvalue of the linear operator, the corresponding eigenfunction must grow like If the similarity solution H(ξ) is to be stable, the eigenvalues of L(H) must be negative. However, there are always two positive eigenvalues, which are related to the invariance of the equation of motion (5) under translations in space and time. Namely, for any ǫ, the translated similarity solution is an equally good self-similar solution of (5), and thus of (11). In particular, we can expand (16) to lowest order in ǫ, and find that is a solution of (13). Thus, since LH = 0, But this means that ν x = β ≡ 1/4 is an eigenvalue of 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 is that upon perturbing the similarity solution, the singularity time as well as the position of 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 0 and t 0 , the dynamics will converge onto a stable similarity solution.
The eigenvalues of the solutions H i have been found numerically in [27]. The result is that the linearisation around the "ground state" solution H 0 only has negative eigenvalues (apart from the two trivial ones), while all the other solutions have at least one other positive eigenvalue. This means that H 0 is the only similarity solution that can be observed, all other solutions are unstable. Close to the fixed point, the approach to H 0 will be dominated by the largest negative eigenvalue ν 1 : For large arguments, the point ξ cr where the correction becomes comparable to the similarity solution is ξ ∼ ǫt ′ν 1 ξ 1−4ν 1 , and thus ξ cr ∼ t ′1/4 . This means that the region of validity of H(ξ) expands in similarity variables, and is constant in real space. This rapid convergence is reflected by the numerical results reported in Fig. 3. More formally, one can say that for any ǫ there is a δ such that if |x ′ | ≤ δ uniformly as t ′ → 0.

Self-similarity of the second kind
In the example of the previous subsection, the exponents can be determined by dimensional analysis, and therefore assume rational values. As Barenblatt [21] points out, 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 assume any value. Often, this value is fixed by some intrinsic property of the equation, resulting in an irrational answer. We will call this situation self-similarity of the second kind (in the sense of Barenblatt). A particularly simple example of this kind of singularity is the pinch-off of a very viscous thread of liquid [30,3], which we present now. Another recent 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 without inertia, for which the equation becomes: 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) has been achieved by writing the thread radius in Lagrangian variables, i.e. as function of a  [33]. Note the long neck.
particle label s. This means the particle is at position z(s, t) at time t, and z t (s, t) is the velocity at time t. The time-dependent constant of integration C(t) must be determined from the constraint that u(s, t) ≡ z s = 1/h 2 (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 are not completely determined from dimensional analysis or from balancing powers of t ′ in the equation of motion. This is a typical situation in which self-similarity of the second kind is observed [32]. Instead, the unknown exponent is determined from a non-linear eigenvalue equation, and takes an irrational value.
Since α = 1 we introduce and Hence where K is an arbitrary constant. Imposing symmetry and regularity of f , we introduce an expansion of f (ξ) of the form (24) to obtain the condition and define β = nβ. Equation (24) can easily be integrated in terms of ln ξ and y = √ f : with C an arbitrary constant. Computing the integral above we obtain where we can fix, without loss of generality, C = 1. In terms of the profile h(s, t) this corresponds to fixing the scale of the spatial variable s. Solutions are undetermined up to such a scale factor, as is clear from the invariance of (21) under a change in spatial scale. 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 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 and hence The function K n (β) may be written explicitly as with roots given by table 2. If one converts the Lagrangian variables back to the original spatial variables, one obtains Thus for t ′ → 0 the typical radial scale t ′ of the generic n = 1 solution rapidly becomes smaller than the axial scale t ′0.175 (cf. table 2). This explains the long necks seen in Fig. 4. For generic initial data u 0 (s) = h −2 0 (s) = B 0 + B 1 s 2 + B 2 s 4 + ... + B j s 2j + ... one expects that B 1 = 0, so that the self-similar solution with n = 1 will develop. Only if B i = 0 for i = 1, 2, ..., n − 1 and B n = 0 the n-th self-similar solution will be the asymptotic description of the solution. For this reason, only the n = 1 solution is stable, since a generic perturbation of the initial data with B i = 0 for i = 1, 2, ..., n − 1 will, in general, make B 1 = 0.  Table 2. A list of exponents, found from K n (β) = 0, with K n given by (28). The number 2n gives the smallest non-vanishing power in a series expansion of the corresponding similarity solution around the origin. Only the solution with n = 1 is stable. The minimum radius is found from (25).

Travelling wave
The pinching of a liquid thread in the presence of an external fluid is described by the Stokes equation [34]. For simplicity, we consider the case that the viscosity η of the fluid in the drop and that of the external fluid are the same. An experimental photograph of this situation is shown in Fig. 1. To further simplify the problem, we make the assumption (the full problem is completely analogous) that the fluid thread is slender. Then the equations given in [1] simplify to where Here we have written the velocity in units of the capillary speed v η = γ/η. The limits of integration z − and z + are for example the positions of the plates which hold a liquid bridge [35]. Dimensionally, one would once more expect a local solution of the form and H out (ξ) has to be a linear function at infinity to match to a time-independent outer solution. In similarity variables, (31) has the form We have chosen z b as a real-space variable close to the pinch-point, such that the similarity 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 similarity variable ξ: with τ = − ln(t ′ ) as usual. This means in similarity variables the solution is a travelling wave. With this modification, the mass balance (30) becomes Now we choose b such that the logarithmic singularity cancels, namely we demand that finite for t ′ → 0. This is achieved by putting Thus defining the similarity equation is finite, and ξ 0 is an arbitrary constant. It remains as an arbitrary axial shift in the similarity solution.
The numerical solution of the integro-differential equation (39) gives The slope of the solution away from the pinch-point are given by which means the solution is very asymmetric, as confirmed directly from Fig. 1.

Centre manifold
In section 2 we described the generic situation that the behaviour of a similarity solution is 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, so higher-order non-linear terms have to be taken into account. The approach to the fixed point is now much slower, and much more of the observed behaviour is determined by the approach to the fixed point. Depending on the type of non-linearity, there is great freedom of possible behaviours, of which we discuss two.  [37]. The recording takes a few minutes, the temperature is 0.32 K. 11 mK. The crystal first "drips" down, so that a crystalline "drop" forms at the bottom (a to c); then a second drop appears (d) and comes into contact with the first one (e); coalescence is observed (f) and subsequently breakup occurs (h).

Quadratic non-linearity
Axisymmetric motion by mean curvature in three spatial dimensions is described by the equation where h(x, t) is the radius of the moving free surface. A very good physical realization of (42) is the melting and freezing of a 3 He crystal, driven by surface tension [36], see Fig. 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 for the problem is that h(0, t) = h(L, t) = R, where R is some prescribed radius. For certain initial conditions h(x, 0) ≡ h 0 (x) the interface will become singular at some time t 0 , at which h(x 0 , t 0 ) = 0 and the curvature blows up. The moment of blow-up is shown in panel h of Fig. 5, for example. Inserting the self-similar solution (2) into (42), one finds a balance for α = β = 1/2. The corresponding similarity equation is One solution of (43) is the constant solution φ(ξ) = √ 2. Another potential solution is one that grows linearly at infinity, to ensure matching onto a time-independent outer solution. However, it can be shown that no solution to (43), which also grows linearly 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 of the dynamical system can be carried out to arbitrary order.
The solution relevant to our analysis is the constant solution, but which of course does not match onto a time-independent outer solution. We thus write the solution as with τ = − ln(t ′ ) as usual. The equation for g is then which we solve by expanding into eigenfunctions of the linear part of the operator It is easily confirmed that where H n is the n-th Hermite polynomial [40]: Thus the linear part of (45) becomes which means that all eigenvalues are negative except for the first, which vanishes. To investigate the approach of the cylindrical solution, one must therefore include nonlinear terms in the equation for a 1 . If we write the equation for a 1 becomes whose solution is Thus instead of the expected exponential convergence onto the fixed point, the approach is only algebraic. Since all other eigenvalues are negative, the τ -dependence of the a i is effectively determined by a 1 . Namely, as we will see below, a j = O(τ −j ), so corrections to (51) are of higher order.
If one linearises around (52), putting a 1 = a (0) This means that the coefficient A of ǫ 1 = A/τ 2 remains undetermined, and a simple expansion of a i in powers of τ −1 yields an indeterminate system. Instead, at quadratic order, a term of the form ǫ 1 = A ln τ /τ 2 is needed. Fortunately, this is the only place in the system of nonlinear equations for a i where such an indeterminacy occurs. Thus all logarithmic dependencies can be traced, leading to the general ansatz where δ i and δ lki are coefficients to be determined. The index n is the order of the truncation. Indeed, the coefficients can be found recursively by considering terms of successively higher order in τ −1 in the first equation: The next two orders will involve the next coefficient a 3 . From (55) and (56) We constructed a MAPLE program to compute all the coefficients up to arbitrarily high order (10th, say). Up to third order in τ −1 the result is: (ln (τ )) 2 √ 2 τ 3 (57) and thus h(x, t) becomes from which one finds the minimum. To second order, the result is Two remarks are in order. First, the presence of logarithms implies that there is some dependence on initial conditions built into the description. The reason is that the argument inside the logarithm needs to be non-dimensionalised using some "external" time scale. More formally, any change in time scalet = t/t 0 leads to an identical equation if also lengths are rescaled according toh = h/ √ t 0 . This leaves the prefactor in (60) invariant, but adds an arbitrary constant τ 0 to τ . This is illustrated by comparing to a numerical 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 τ 2 . As seen in Fig.6, the remainder is varying slowly over 12 decades in t ′ . If the constant τ 0 is adjusted, this small variation is seen to be consistent with the logarithmic dependence predicted by (60). The second important point is that convergence in space is no longer uniform as implied by (20) for the case of self-similarity of the first kind. Namely, to leading order the pinching solution is a cylinder. For this to be a good approximation, one has to require that the correction is small: ξ 2 /τ ≪ 1. Thus corrections become important beyond ξ cr ∼ τ , which, in view of the logarithmic growth of τ , implies convergence in a constant region in similarity variables only. As shown in [36], the slow convergence toward the self-similar behaviour has important consequences for a comparison to experimental data.

Cubic non-linearity
The next example is that of bubble breakup [41], for which a very different form of nonlinearity is observed. As shown in [41], the equation for a slender cavity or bubble is where a(z, t) ≡ h 2 (z, t). The integral runs over the fluid domain. If for the moment one disregards boundary conditions looks for solutions to (61) of cylindrical form, a(z, t) = a 0 (t), one can do the integral to find Figure 7. The pinch-off of an air bubble in water [44]. An initially smooth shape develops a localised pinch-point. It is easy to show that an an asymptotic solution of (62) is given by corresponding to a power law with a small logarithmic correction. Indeed, initial theories of bubble pinch-off [42,43] treated the case of an approximately cylindrical cavity, which leads to the radial exponent α = 1/2, with logarithmic corrections. However both experiment [44] and simulation [41] show that the cylindrical solution is unstable; rather, the pinch region is rather localised, see Fig. 7. Therefore, it is not enough to treat the width of the cavity as a constant L; the width ∆ is itself a timedependent quantity. In [41] we show that to leading order the time evolution of the integral equation (61) can be reduced to a set of ordinary differential equations for the minimum a 0 of a(z, t), as well as its curvature a ′′ 0 . Namely, the integral in (61) is dominated by a local contribution from the pinch region. To estimate this contribution, it is sufficient to expand the profile around the minimum at z = 0: a(z, t) = a 0 + a ′′ 0 /2z 2 + O(z 4 ). As in previous theories, the integral depends logarithmically on a, but the axial length scale is provided by the inverse curvature ∆ ≡ (2a 0 /a ′′ 0 ) 1/2 . Thus evaluating (61) at the minimum, one obtains [41] to leading orderä which is a coupled equation for a 0 and ∆. Thus, a second equation is needed to close the system, which is obtained by evaluating the the second derivative of (61) at the pinch point: The two coupled equations (64),(65) are most easily recast in terms of the timedependent exponents where τ ≡ − ln t ′ and β = α − δ, are a generalisation of the usual exponents α and β. The exponent δ characterises the time dependence of the aspect ratio ∆. Returning to the collapse (62) predicted for a constant solution, one finds that α = 1/2 and δ = 0.
In the spirit of the the previous subsection, this is the fixed point corresponding to the cylindrical solution. Now we expand the values of α and δ around their expected asymptotic values 1/2 and 0: To leading order in δ, the resulting equations are which describe perturbations around the leading-order similarity solution. These equations are analogous to (51), but they have a degeneracy of third order, rather than second order. Equations (68) are easily solved to yield, in an expansion for small δ [41], Thus the exponents converge toward their asymptotic values α = β = 1/2 only very slowly, as illustrated in Fig. 8. This explains why typical experimental values are found in the range α ≈ 0.54 − 0.58 [44], and why there is a weak dependence on initial conditions [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 is easy to re-analyse the solution to the mean curvature problem, and to formulate it in terms of time-dependent exponents. Let us define a 0 and a ′′ 0 for the mean curvature problem through h 2 (z, t) = a 0 + a ′′ 0 /2z 2 + O(z 4 ). From (44) it follows that a 0 = 2t ′ and a ′′ 0 = 2( √ 2 + g)g ′′ . This gives a ′′ 0 = 4 √ 2a 1 , and thus instead of (68). This means mean curvature flow, formulated in terms of exponents, once more gives a quadratic non-linearity, rather than the cubic term (68) found for bubble collapse. From (70) one finds to leading order δ mc = 1/(2τ ) and α mc = 1/2 + O(1/τ 2 ) for the mean curvature flow.

Limit cycles
An example for this kind of blow-up was introduced into the literature in [12] in the context of cosmology. There is considerable numerical evidence [46] that discrete selfsimilarity occurs at the mass threshold for the formation of a black hole. The same type 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 cosmological problem leads to the following system: In [13], the self-similar description corresponding to the system (71)-(75) was solved using formal asymptotics and numerical shooting procedures. This leads to the solutions observed in [12]. Below we propose a very similar system, which we solve analytically, and which shows the same type of limit cycle behaviour: The simplified system (76)-(78) can be solved introducing characteristics: which leads to The system (80) has to be solved at constant ν with initial conditions 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 where f init (x) is the initial profile, which can be written in terms of ν, η as The function f can be eliminated from (80) using the transformation which results in the system for a harmonic oscillator: The solution is with C(ν) = u 2 init (ν) + v 2 init (ν) and φ 0 (ν) = arcsin(u init (ν)/C(ν)). According to (82), a singularity of the PDE system (76)-(78) first occurs at the maximum f 0 of f init (x), which we assume to occur at x = 0 without loss of generality. Thus locally we can write f (x, 0) ≈ f 0 − ax 2 , and for small t ′ = t 0 − t we have for any finite ξ. Using this explicit form of f , (84) can be integrated to find ζ. At constant ξ we have so taking the limit t ′ → 0, the leading order result for (84) is Thus as the singularity is reached, t ′ → 0, the variable ζ goes to infinity. This means the singularity corresponds to the long-time limit of the dynamical system (85), which will perform a harmonic motion according to (86). Namely, for t ′ → 0 the result is Thus the singular solution is of the general form where ψ is periodic in τ . This is a particularly simple version of discretely self-similar behaviour. Note that the character of the solution is different from the travelling waves described in section 3.

Strange attractors and exotic behaviour
In connection to limit cycles and in the context of singularities in relativity, a few interesting situations have been found numerically quite recently. One of them is the existence of Hopf bifurcations where a self-similar solution (a stable fixed point) is transformed into a discrete self-similar solution (limit cycle) as a certain parameter varies (see [49]). Other kinds of bifurcations, for example of the Shilnikov type, are found 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 differential equations for the coefficients of the expansion of an arbitrary perturbation to an explicit solution. Such high-dimensional systems in principle allow for a rich variety of dynamical behaviours, including those found in classical finite dimensional dynamical systems, such as chaos. Consider for instance an equation for the perturbation g (the analogue of (45)) of the form where Lg is a linear operator. Assuming an appropriate non-linear structure for the function 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 be generalised to produce any low-dimensional dynamics near the singularity. Namely, let us generalise the system (76)-(78) to Using the transformation (79), the first n equations are turned into: which is an ODE system for constant ν. The system (95) has to be solved with initial conditions init (ν), i = 1, . . . , n.

Outlook
The singularities described in this paper are point-like in the sense that they occur at a single point x 0 at a given time t 0 . There are two important situations we have not discussed and for which the dynamical systems point of view and analytical approach does not apply. First, the case of singularities not developing at single points, but on sets of finite measure. This is the case in a few simple examples of reaction-diffusion equations of the family where depending on values of p > 1 and q > 1 singularities in the form of blowing-up u may be regional (u blows up in subsets of Ω of finite measure) or even global (the solution blows-up in the whole domain). See for instance [53] and references therein. Another interesting possibility is that the stable fixed point which is approached depends on the initial conditions. Thus exponents could vary either discretely or continuously with the choice 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 fixed points 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 of the singular set at the time of first blow-up is at most 1 (cf. [55]). This connects with the second issue we did not address here. It is the nature of the singular sets both in space and time. In many instances, existence of global in time (for all 0 ≤ t < ∞) solutions to nonlinear problems can be established in a weak sense, that is, allowing certain kind of singularities 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 the singular set is well-known [55]. Hence, provided a certain singularity does not persist in time, the question is how to continue the solutions after a singularity has developed.