An asymptotic study of blow up multiplicity in fourth order parabolic partial differential equations

Blow-up in second and fourth order semi-linear parabolic partial differential equations (PDEs) is considered in bounded regions of one, two and three spatial dimensions with uniform initial data. A phenomenon whereby singularities form at multiple points simultaneously is exhibited and explained by means of a singular perturbation theory. In the second order case we pre- dict that points furthest from the boundary are selected by the dynamics of the PDE for singularity. In the fourth order case, singularities can form simultaneously at multiple locations, even in one spatial dimension. In two spatial dimensions, the singular perturbation theory reveals that the set of possible singularity points depends subtly on the geometry of the domain and the equation parameters. In three spatial dimensions, preliminary numerical simulations indicate that the multiplicity of singularities can be even more complex. For the aforementioned scenarios, the analysis highlights the dichotomy of behaviors exhibited between the second and fourth order cases.

The corresponding fourth order problem (1.1b), which had previously attracted somewhat less attention, has more recently enjoyed significant interest. In [18], the critical exponent p c (n) = 1 + 4/n was established for problem (1.1b) with the nonlinearity f (u) = |u| p and Ω = R n . The existence and stability of blow-up profiles to (1.1b) with power and exponential nonlinearities was investigated in [16,17,4] BLOW UP SETS OF FOURTH ORDER PDES 3 where it was found that, in contrast to the second order problem, there exists stable self-similar singularity profiles. In the context of Micro-Electro Mechanical Systems (MEMS), where for f (u) = 1/(1−u) 2 problem (1.1b) models the deflection of a beam/plate under Coulomb forcing, quenching solutions have been studied in [24,7,8]. Related higher order problems in the context of thin film dynamics have been studied extensively ( [11,13,14,12,15]) where de-wetting processes occur via ruptures in degenerate fourth order parabolic equations.
Fourth order problems exhibit many interesting and surprising solution features when compared to their second order counterparts. As a simple example, consider the reduced setting where Ω = R n and f (u) = 0, and the evolution problem If ψ(x) ∈ C 0 ∩ L ∞ (R n ), then the unique global solution (cf. [26,27]) of (1. In contrast to the Gaussian kernels of the corresponding second order heat equation, the Bessel functions J ν in the integrand of (1.4b) generate highly oscillatory solution behavior. Consequently, many second order features such as the positivity preserving property (ψ > 0 implies u > 0) and the maximum principle do not extend to higher order problems.
As a result of these oscillatory features, the blow up dynamics for the fourth order problem (1.1b) are quite different to that of the well-studied second order problem (1.1a). To demonstrate the contrasting behaviors of these two problems, consider the example case of the 1D strip Ω = [−L, L], with nonlinearity f (u) = e u and uniform zero initial data ψ(x) = 0. In the second order case, it is well known (cf. [9]) that for L large enough, no equilibriums solutions are present and, from maximum principle considerations, that the solution blows up uniquely at the origin. In Fig. 1, numerical solutions of the equivalent fourth order problem (1.1b) are displayed for values L = 5 and L = 7 and are observed to be very different from the previously mentioned second order behavior. For L = 5, the blow-up point is observed to occur uniquely at the origin while for L = 7, we observe two singularities forming simultaneously at distinct points. This multiple singularity phenomenon of (1.1b) was recently observed in [7] for the MEMS case f (u) = 1/(1 − u) 2 and radially symmetric solutions in one and two dimensions. In the 2D radially symmetric case, the singularities form simultaneously along a ring of points for L sufficiently large.
In this particular application, the singularities indicate the contact points between two elastic surfaces and so their multiplicity and location is of practical importance.
The multiplicity of quenching singularities for the inverse square nonlinearity has also been investigated (cf. [8]) for general 2D geometries. As before, an illustration of the dichotomy between the second (1.1a) and fourth (1.1b) order cases is provided through the f (u) = e u case and the square region Ω = [−L, L] 2 . For L sufficiently large, problems (1.1) blow-up in finite time. In the second order case, solutions of (1.1a) blow-up uniquely at the origin. As indicated by the numerical simulations of (1.1b) shown in Fig. 2, the multiplicity of singularities is remarkably different for L = 1.5 and L = 1.8. In the former case, the singularity occurs uniquely at the origin, while in the later case, the singularity occurs simultaneously at four distinct points.
In the recent works of [7,8], the location and multiplicity of singularities in problem (1.5b) was studied for the nonlinearity f (u) = 1/(1 − u) 2 . In the present work, the contribution is threefold: First, we consolidate and generalize the results of [7,8] to the fourth order problem (1.1b) with general positive convex nonlinearities in one and two dimensions. The result of this analysis is a geometric framework which predicts the multiple singularity phenomenon in (1.1b) for a general class of nonlinearities f (u). The analysis demonstrates that this phenomenon is due to a combination of boundary effects and the dynamics of fourth order problems which do not admit a maximum principle.
In addition, a geometric framework is established for predicting the possible blow up set of (1.1b), for general regions Ω.
Second, we apply the same geometric framework developed herein to the second order problem (1.1a) with general nonlinearities. This analysis elucidates the underlying reasons why the multiple singularity phenomena does not occur generically in (1.1a). In addition, it provides a simple geometric framework for predicting the singularity location in such problems. Consequently, we establish an understanding of how the location of singularities in problems (1.1) is determined by the geometry of Ω and the dynamics of the PDE in one and two dimensions.
Third and finally, we present numerical simulations of (1.1b) in three dimensions to gain insight into the corresponding multiple singularities phenomenon in higher spatial dimensions. The preliminary simulations suggest that the multiplicity of singularities in (1.5b) may be greater than in 2D. For example, in the cube Ω = [−L, L] 3 , numerical simulation show the number of blow up points changes from eight to one as L increases through L c . A natural conjecture based on observations is the following: There exists an L c (n) such that (1.1b) for Ω = [−L, L] n and ψ = 0, exhibits blow-up at 2 n distinct points whenever L > L c (n).
To analyze the multiple singularity phenomenon, and thus obtain predictions of the multiplicity and location of singularities in the system (1.1), we employ formal asymptotic methods in the limit of large domain size. If the length scale of the 6

ALAN E. LINDSAY
domain Ω is L, then the problem of interest, as motivated by the previous examples, is to understand the singularity set for large L, or equivalently small ε = L −1 . When the rescaling Ω → ε −1 Ω is applied to (1.1), the second order singularly perturbed initial boundary value problem and its fourth order equivalent are arrived at. For the purposes of the present work, general positive, convex source terms f satisfying are considered. In this particular formulation, problem (1.5a) is a generic model for a slowly diffusing concentration field with local reaction kinetics f (u). The fourth order problem (1.5b) is a ubiquitous model for deflection of a plate with small flexural rigidity [39,7], undergoing forcing f (u). e t , f (t) = (1+t) p for p > 1 will be used. The premise of our approach is to construct explicit short time solutions to (1.5) in the limit ε → 0, the critical points of which act as surrogates for the singularity location(s) of (1.5) for one dimensional strips and general bounded regions of R 2 .
Note that in (1.5a-1.5b), the initial conditions have been chosen to be uniformly zero. In general, the particular form of the initial condition will play a large role in determining the singularity set of the system. In the present work, however, we restrict our attention to the uniform initial condition case to specifically study how the geometry of Ω and the dynamics of the underlying PDE alone combine to determine the singularity set.
Heuristically, the dynamics of problems (1.5) can be decomposed into two temporal regimes. The first is a short time regime whereby the solution is shaped by the initial data, the geometry of the region and the dynamics of the governing equation. The second is the blow up regime (t − T ε ) 1, in which the solution is changing very rapidly in a highly localized vicinity of the singularity points. Whichever point(s) enter the basin of attraction of the blow-up regime first, will eventually be the blow-up locations of the system. This mechanism promotes single point blow-up which is the standard behavior observed for problems (1.5) with general initial data -consequently simultaneous multiple point blow-up is not generically stable to asymmetric perturbations in the initial condition. This instability is discussed further in the disc geometry example of §4.1 and the square example of §4.2. Once a uniformly valid asymptotic expansion is established, its validity is confirmed with comparison to numerical simulations and found to be good, even for moderately large times. In the second order case (1.5a), we find that the global maximum of the uniformly valid solution is x = 0 for all ε. This indicates, as shown in [9], that this point is selected by the dynamics of the full PDE for singularity. In contrast, when this asymptotic theory is applied to the fourth order case (1.5b), the location and multiplicity of the global maxima are found to depend on the value of ε.
2.1. Laplacian Case. In this section, a small amplitude asymptotic solution to the PDE is developed. In the outer region, away from x = ±1, the solution is spatially uniform and satisfies Boundary layers are required in the vicinity of x = ±1 to enforce the conditions u(±1) = 0. The formulation for the boundary at x = 1 is established through the stretching variables Substituting variables (2.3) into (2.1) yields the equation This equation further simplifies from noticing that u 0 = O(t) as t → 0 and that f (t) = 1 + O(t) from the assumptions (1.5c) on f (t). Therefore, a short time and small amplitude approximation of (2.4) satisfies The solution of the linear problem (2.5a) and its far field asymptotic behavior are is established. The asymptotic behavior (2.5b) of the profile v(η) is increasing monotonically to the limiting value which indicates that the solution will have larger value in regions farther away from the boundary. Therefore, the regime of interest in the global approximation (2.6) is when 1 ± x φ. After applying the far field behavior (2.5b) relevant to this regime, the solution with exponential corrections, valid away from x = ±1, is obtained. As lim t→T0 φ(t; ε) = ∞, this expression is defined for 0 < t < T 0 , however, we can only expect good quantitative validity for t 1. Indeed, from , and calculate that the relative L ∞ error estimate for the asymptotic approximation is We now posit that the global maximum x = 0 of (2.7) provides a predictor of the location of the singularity of the full problem (2.1). Heuristically, this maximum of the small amplitude solution will enter the basin of attraction of a stable similarity solution regime of (2.1) before any others.
For this case, the purely reaction problem u 0t = (1 + u 0 ) 2 , u 0 (0) = 0 has the solution u 0 (t) = t/(1 − t) which blows up at T 0 = 1. Note that since u 0 (t) is a supersolution of (2.9), standard comparison principles imply u 0 (t) > u(x, t) so that For more details on upper and lower bounds for T ε , see [10]. In Fig. 3, a comparison between the numerical solution of (2.9) and the asymptotic solution Fig. 3(a) excellent agreement is observed for even moderately large values of t. To obtain accurate numerical solutions of (2.9) very close to singularity, r-adaptive moving mesh methods are employed together with computational time stepping. For more details, see [5].   To establish a boundary layer in the vicinity of x = 1, the stretching variables are introduced. The leading order term in the expansion as φ → 0, is the profile While a closed form solution to (2.12) is available, it is rather unsightly and not very useful. The key feature is the far field behavior, which can be obtained from a WKB analysis. By applying the large η anzatz v(η) at leading order. This ODE admits three non-trivial solutions however, only exp[ψ 1 ] and exp[ψ 2 ] decay as η → ∞. Therefore the complete specification of (2.12) with the obtained far field behavior is where A, θ are arbitrary constants and ω = 3 · 2 −11/3 . It is crucial to notice that the far field behavior of the equivalent second order problem (2.5b) is strictly increasing towards the limiting value, while the solution of the fourth order problem (2.15b) exhibits oscillations about it (cf. Fig. 4). This loss of monotonicity in the fourth order case has a dramatic effect on the location of the singularity selected by the dynamics of the full PDE.    .15). The important distinction between the two cases is that the second order profile is monotone while the fourth order attains a global maximum at η = η 0 .
The uniformly valid short time asymptotic solution to (2.10) is now given by where v(η) is the profile determined in (2.15). The aim now is to understand which value(s) x c (t) ∈ (−1, 1) will be global maxima of expression (2.16). In the previ- In addition, if we posit that the global maxima of (2.16) act as surrogates for the singularity points selected by the dynamics of the full PDE, then a crude approximation for such points is where T ε is the finite time singularity of the full PDE (2.10). It is crucial to note here that for φ c to be well-defined, the fact that T ε < T 0 for sufficiently small ε has been used in (2.17). In other words, the PDE problem (1.5b) can blow up faster than the ODE problem u 0t = f (u 0 ).
In the second order case, the lower bound T 0 ≤ T ε is a simple consequence of the maximum principle. In the fourth order problem (1.1b), this approach to bounding T ε from below is not valid. The challenge of obtaining lower bounds to the blowup time of (1.1b) in Ω = R N via a classical comparison argument with u 0 (t) was discussed in [18,22]. The authors' solution was to instead compare with the socalled order-preserving majorizing equation for which a comparison principle can be established. It is also known that the blow-up time of (1.5b) satisfies for constant C > 0 [33]. It would be interesting to extend these results and ideas to rigorously establish a lower bound on T ε and therefore show that for problem (1.5b) on bounded regions, T ε < T 0 for certain ranges of ε.
very close to the singularity time T ε , the r-adaptive techniques described in [5,6,7] are employed. As shown is Fig. 5(a), the leading order asymptotic prediction of (2.16) is accurate, even for moderately large values of time. The relative error of the approximation satisfies which is observed in Fig. 5(b). It is interesting to note here that T ε < T 0 = 1 for the value ε = 0.1 used here. In Fig. 6(b), the asymptotic approximation is seen to provide a rather crude quantitative prediction of the singularity point x c (T ε ). In order. However, the main goal herein is to demonstrate the underlying principle behind the phenomenon, a purpose for which the leading order theory suffices.    In panel (b), a comparison between the singularity points as predicted by the asymptotic formula (2.17) and full numerical simulations, is displayed. 3. Two Dimensional Theory. The generalization of the previous analysis to bounded two-dimensional star-shaped regions is facilitated by implementing an arclength tangent coordinate system (ρ, s), where ρ > 0 measures the distance from x ∈ Ω to ∂Ω, whereas on ∂Ω the coordinate s denotes arc-length along the boundary.
In this co-ordinate system, the Laplacian operator admits the representation where κ(s) indicates the curvature of ∂Ω as a function of arc-length s along it. In this section, we analyze the second order problem and its fourth order equivalent where Ω T and ∂Ω T are defined in (1.1c). As previously, we posit that a uniform solution u 0 (t) satisfying and expansions in terms of φ 1. After substituting (3.4) and the expansion v = v 0 + φv 1 + · · · into (3.2a), we arrive at the problems (3.5b) The correction equation (3.5b) also admits the decomposition v 1 (η, s) = κ(s)v 1 (η) The leading order problem (3.5a) depends only on the perpendicular distance from the boundary while the correction term (3.5b) incorporates a dependence on the curvature of ∂Ω. In fact equation (3.5a) is precisely the one dimensional boundary profile established in (2.5) and in particular, its far field asymptotic behavior is given by (2.5b).
We now turn to the problem of establishing a uniformly valid global solution at a general x ∈ Ω. The previous analysis reveals that to leading order, the solution is principally determined by perpendicular distances from the boundary. Therefore to construct the solution at x ∈ Ω, it is necessary to determine all boundary points y ∈ ∂Ω such that the straight line between x and y is contained in Ω and meets the boundary orthogonally at y ∈ ∂Ω. With this in mind, suppose that for some x ∈ Ω there are n boundary points {y 1 , . . . , y n } ∈ ∂Ω such that the straight line l(x, y j ) segment between x and y j is contained in Ω and meets ∂Ω orthogonally at y j , then a two term asymptotic solution u(x, t) is given by In (3.7), u 0 (t) is the solution of the ODE (3.3) and κ(y j ) is the curvature of ∂Ω at y j ∈ ∂Ω. Equation (3.7) represents a uniformly valid asymptotic expansion, derived in the limit φ → 0. Consequently, it is important to remark that the validity of this solution is restricted to regions Ω for which κ = O(1) as φ → 0. In other words, the theory developed here is not valid for regions with rough or spiky boundaries ∂Ω.
The objective now is to deduce the global maxima of the asymptotic formulation (3.7) with the understanding that these critical points provide surrogates for those selected by the dynamics of the full PDE as the location of a finite time singularity.
From the fact that the profile v 0 (η) is strictly increasing in η, the leading order theory predicts that the maximum will be away from ∂Ω and in the region where |x − y j | φ. Therefore, in the case φ 1, we are primarily interested in evaluating (3.7) at points x ∈ Ω for which v 0 (|x − y j |/φ) and v 1 (|x − y j |/φ) can be replaced by their large argument behavior (2.5b). This reduces (3.7) to at points x ∈ Ω such that |x − y j | φ. The form of (3.8) can also be recovered by applying the anzatz to (1.5a) and obtaining equations for Ψ(x) and A(x) in the limit φ → 0. The simplicity of (3.8) allows several deductions to be made regarding the points x ∈ Ω at which u(x, t) attains its maximum value and therefore the points which can be expected to reach singularity before others under the dynamics of the full PDE (2.1). As indicated by the form of (3.8), the distances |x − y j | play a large role in these determinations. As the correction terms to the uniform state are exponentially small, larger values of |x − y j | will have the smallest contribution to (3.8) leaving the smaller values of |x − y j | to dominate the solution. Therefore, the maxima of (3.8) will be located at points x ∈ Ω which are furthest from the boundary, i.e., This result can be heuristically reconciled with the intuitive notion that, since the boundary condition u = 0 on ∂Ω inhibits blow-up, and the solution is monotone away from the boundary, that points furthest from ∂Ω are those more likely to develop a singularity as t → T − ε .
In singularly perturbed elliptic problems, the distance function also plays a key role in the determining the location(s) of solution concentration. For example, it is well known (cf. [36,37,38] and the references therein) that in the problem After expanding (3.11) with v = v 0 + φv 1 + · · · and collecting terms at the relevant order, v 0 and v 1 satisfy (3.12c) As with the previous analysis in the Laplacian case, a uniformly valid asymptotic solution is constructed from the flat outer solution u 0 (t) and the boundary contributions established by (3.11-3.12). Therefore the short time asymptotic solution at x ∈ Ω is given by where {y j } n j=1 ∈ ∂Ω are the points for which the straight line l(x, y j ) between x and y j is contained in Ω and meets ∂Ω orthogonally at y j . An illustration of how the points {y j } n j=1 ∈ ∂Ω depend on any particular x ∈ Ω is given in Fig. 7. The goal now is to predict the location of blow up points by determining the global maxima at t = T ε of (3.13) for a range of bounded two-dimensional regions Ω.
The explicit form of (3.13) indicates that the quantities s j = |x − y j | play a key role in determining whether or not a given point x ∈ Ω will be selected by the dynamics of the full PDE for singularity. As the profile v 0 (η) satisfying (3.12a) has a global max at η = η 0 , it follows from (3.13) that points x ∈ Ω such that |x − y j | = η 0 φ(t; ε) are, to leading order, local maxima. Therefore, at time t and for fixed ε, we define the set ω(t) ⊂ Ω where (3.14) The set ω(T ε ) describes, to leading order in (3.13), points which are local maxima and are therefore more likely to be selected for singularity by the dynamics of the full PDE (3.2b). In general, singularities will not form simultaneously on all points of ω(T ε ). Instabilities along ω(T ε ) in combination with the higher order curvature effects present in (3.13), will increase the value of u(x, t) at certain discrete points along ω(T ε ).
Amongst the points described by ω(T ε ), there may be special points x ∈ Ω which receive multiple boundary contributions, i.e. there are y 1 , y 2 ∈ ∂Ω such that d(x, ∂Ω) = d(x, y 1 ) = d(x, y 2 ). At such points, the leading order terms in the sum (3.13) combine to increase the value of u(x, t) therefore increasing the possibility that that point is selected for singularity. This motivates the definition of the skeleton of the domain, S Ω . The skeleton S Ω is the set of points x ∈ Ω for which there are at least two points y 1 , y 2 ∈ ∂Ω such that d(x, y 1 ) = d(x, y 2 ) and such that the straight lines l(x, y 1 ) and l(x, y 2 ) are both contained in Ω and meet ∂Ω orthogonally. More compactly, l(x, y 1 ), l(x, y 2 ) ∈ Ω, l(x, y 1 ) ⊥ ∂ τ (y 1 ), l(x, y 2 ) ⊥ ∂ τ (y 2 ).
where ∂ τ (y) denotes the unit tangent vector to ∂Ω at y ∈ ∂Ω. It is important to note that for a particular x ∈ S Ω , there may be several sets of boundary points which Figure 8. The two structures ω(t) and S Ω , defined in (3.14) and (3.15) respectively. Panel (a) The set ω(t) consists of the points in Ω that are at the distance η 0 φ(t; ε) from ∂Ω. Panel (b) The skeleton S Ω is the set of points in Ω that are equidistant to two or more points on ∂Ω such that the line segments between x and those points are in Ω and meet ∂Ω orthogonally.
mutually satisfy the conditions specified in (3.15). To make this notion concrete, consider the simple example of the elliptical region For this particular region, (0, 0) ∈ S Ω as both the sets of points {(−a, 0), (a, 0)} and {(0, −b), (0, b)} satisfy the conditions of (3.15). However, as will subsequently become apparent, it is necessary to distinguish between such sets of contributing points by their distances to x ∈ S Ω . For this region, consider the function s Ω : S Ω → R such that So for a skeleton point x ∈ S Ω , which may have multiple sets of contributing boundary points satisfying (3.15), the function s Ω (x) is the shortest distance from x ∈ S Ω to any of its contributing boundary points.
As t increases, the set ω(t) propagates inwards from ∂Ω towards the center of the domain. At time t = T ε , the value of u(x, t) predicted by (3.13) will be quite different depending on whether ω(t) has intersected S Ω or not. To differentiate between these two distinct scenarios, the skeleton arrival time T S is introduced where The value of T S indicates the shortest time at which ω(t) intersects with the skeleton S Ω . The dichotomy of possible singularity locations predicted by (3.13) is now the following; If T ε < T S , then singularities are predicted to occur on ω(T ε ), typically at discrete points selected by the higher order curvature effects of (3.13). On the other hand, for T ε ≥ T S , equation (3.13) predicts that the singularity will occur on S Ω at the point(s) for which s Ω (x) = η 0 φ(T ε ; ε).

4.
Application of the asymptotic theory. In the following subsections, we apply the aforementioned asymptotic theory to predict the singularity set of problems  To apply the theory developed in §3.2 for the fourth order case (1.5b), we first note that the skeleton of the domain is S Ω = {(0, 0)} and therefore T S > 0. The predicted behavior is the following. If T ε < T S , singularities will form on ω(T ε ) which for this example is a ring or radius 1 − η 0 φ(T ε , ε) where φ(t; ε) = εu 1/4 0 and u 0t = f (u 0 ). For T ε > T S , the singularity is predicted to occur at the origin. As the following example will clarify, this situation only materializes in the radially symmetric case -in general the solution develops an instability along ω(t), which results in a single point selected for singularity.
The consideration of radially symmetric solutions to (1.5b), along with the power nonlinearity f (u) = (1 + u) 2 , yields the reduced problem where r = x 2 1 + x 2 2 . In Fig. 9(a), a solution profile of (4.1) close to blow-up ( u ∞ = 10 10 ) is displayed for ε = 0.1. As predicted by the asymptotic theory, the solution is observed to blow-up simultaneously along an inner ring of points. The dependence of the radius of the blow-up ring and ε is illustrated in Fig. 9(b).  However, this blow up ring solution is not stable in a full two dimensional setting.
A similar instability was described in [12] for ring rupture in thin-film equations. In Fig. 10, we display a full 2D (non-radially symmetric) solution of (1.5b) initialized with small amplitude noise. As Fig. 10(a) indicates, the initially noisy data is smoothed out by the dynamics of the PDE and the structure of ω(t) emerges. An instability develops along ω(t) in the angular direction, consequently the dynamics of the PDE selects a single point on ω(t) for blow-up, as shown in Fig. 10(b).
A repetition of this numerical experiment for different realizations of initial data, results in blow-up points which are uniformly distributed along ω(t).

4.2.
Example: Square. We now compare the predictions of the theory for the second order case and the fourth order case for the square geometry Ω = [−1, 1] 2 .
In the second order case (1.5a), the short time asymptotic theory predicts the critical point to be x c = max x∈Ω d(x, ∂Ω), which gives x c = {(0, 0)} in this scenario.
To apply the theory for the fourth order problem, the first step is to construct the skeleton S Ω . For this example, S Ω consists of the x 1 x 2 axes and the diagonals of the square, i.e, An inspection of the skeleton (cf. Fig. 11(a)) reveals that T S = 0 and so the singularities are predicted to form on S Ω .     As the values of s Ω (x) on the diagonal portion of S Ω are lower than those lying on the axes portion, the theory predicts that singularities will form simultaneously at four points on this diagonal for ε < ε c and at the origin for ε > ε c . Moreover, the prediction of the x 1 locations of the singularities are given by where the critical value ε c is implicitly determined by 1 = η 0 φ c (ε c ). The x 2c values are then ±x 1c to give a total of four distinct points in the case ε < ε c .   In Fig. 11(b), we see qualitative agreement between the numerical and asymptotic predictions and good quantitative agreement when ε is small. However, the asymptotic theory does not predict ε c , the threshold between single and multiple blow-up, with high accuracy. This is not surprising since the assumptions which underpin the asymptotic theory: a uniform central region coupled to a propagating boundary effect, do not hold when ε ≈ ε c .
In the presence of a small random perturbation to the initial data, the four point blow up configuration is not generically stable. Indeed, for simulations of (1.5b) initialized with small random noise, the short time solution will develop four peaks, but small discrepancies in their amplitude will result in one being selected for blowup by the PDE dynamics, before the remaining peaks are able able to fully develop.
Over many realizations, one recovers that each of the four possible blow-up locations is selected with uniform probability 1/4. symmetries of the domain. Second, that the multiplicity of singularities can change more than once for ε ∈ (0, ε c ), in contrast to behavior seen in the previous 1D and square example. For the second order problem (1.5a), the asymptotic theory predicts that singularities should form at the origin in the absence of any noise.
In the fourth order case, we begin by considering the partial skeleton of the domain which takes the appearance of an "envelope" (cf. Fig. 13(a)). From the asymptotic theory, we therefore predict that in the absence of noise, problem (1. The asymptotic theory consequently predicts the existence of two critical values ε 1 < ε 2 < ε c such that four singularities occur for ε ∈ (0, ε 1 ), two singularities occur when ε ∈ (ε 1 , ε 2 ) and one singularity occurs when ε ∈ (ε 2 , ε c ). The numerical experiments displayed in Fig. 13 illustrate these three outcomes for problem (1.5b) with the nonlinearity f (u) = e u .
For the rectangular domain, we therefore have that the blow up set of (1.5b) consists of either one, two or four points depending on the value of ε. Moreover, the skeleton theory is able to capture possible singularity configurations of (1.5b) that symmetry    Table. 1 the L ∞ error between the asymptotic and numerical blow up point predictions is shown to be very small, indicating the asymptotic and numerical predictions are in good agreement. The numerical simulations find the maximum value to occur at the same numerical node point for each value of ε -hence the errors are identical for each value of ε. This is exactly as predicted by the leading order asymptotic prediction (3.9), which is independent of ε.
As seen in the previous examples, the application of the leading order asymptotic theory to the fourth order problem (1.5b) is considerably more delicate. The first step is to numerically calculate the skeleton S Ω for the region, which is displayed in Fig. 14(a). An inspection of S Ω indicates that T S > 0, i.e. the skeleton arrival time is positive. Therefore, for T ε < T S the leading order asymptotic theory predicts  Table 1. Accuracy of blow-up point predictions for problem (1.5a) for region (4.4) with f (u) = (1+u) 2 and various ε. The asymptotic prediction x * c = (0.3070, −0.0345) is determined from (3.9) and the estimate x c is obtained from the maximum of the numerical solution profile at u ∞ = 1 × 10 4 . blow-up on discrete points of ω(T ε ) selected by curvature effects. For T ε ≥ T S , the theory predicts blow-up at point(s) x ∈ S Ω for which s Ω (x) = η 0 φ(T ε ; ε).  In Fig. 14(b), the numerically obtained blow-up points for different values of ε are shown overlaid on the skeleton S Ω . The direction of the arrows indicate increasing values of ε. Note that for this example, the relevant portion of S Ω is composed of three distinct branches, labelled I, II and III. The blow-up points corresponding to the smallest values of ε are to the left of segment I. Indeed, T ε < T S for these values of ε and the asymptotic theory predicts that blow-up occurs on discrete points of ω(T ε ) selected by the boundary points of largest curvature (cf. Fig. 15(a)).
As the value of ε increases, T ε exceeds T S , and blow-up occurs initially on segment I of S Ω , in agreement with the asymptotic theory (cf. Fig. 15(b)). In addition, the asymptotic theory predicts that blow-up should occur at x ∈ S Ω such that s Ω (x) = η 0 φ(T ε ; ε). As the shading in Fig. 14(a) represents values of s Ω , we see that the condition s Ω (x) = η 0 φ(T ε ; ε) can be satisfied by multiple x ∈ S Ω so that the leading order theory predicts multiple simultaneous blow up over a range of ε. However, the leading order terms in (3.13) only predict these points to be local maxima while higher order terms in the expansion (3.13), relating to the boundary curvature, will further increase the magnitude of the solution at certain points. It therefore follows that those discrete points with the largest value will be selected by the dynamics of the PDE for blow-up. Consequently, the interaction of the leading and first order terms can result in the blow-up location switching between the three segments of S Ω as ε is increased.
However, as the blow-up point jumps from segment I to II, there is a finite value of ε for which multiple blow-up occurs and similarly with the transition from II to III (cf. Fig. 15(c)-15(e)). From Fig. 15(b), 15(d), we see that when the blow-up occurs on segment I and II of S Ω respectively, the peak on the subsequent segment of S Ω is developing and will eventually form the blow-up point for larger ε.
5. Three Dimensions. The analysis developed so far explains the multiple singularity phenomenon of (1.5b) by means of a propagating non-monotone boundary layer. Propagating boundary effects from distal segments of ∂Ω combine to raise the solution value at certain points in Ω which are in turn selected by the dynamics of the PDE for singularity. This understanding naturally extends to three and higher dimensions, which we now investigate briefly for (1.5b) on the cubic region Ω = [−1, 1] 3 . Extending the analogy of one and two dimensions, we expect there exists an ε c such that as ε increases through ε c , the multiplicity of singularities increases. As the cube has eight corners, the multiplicity can be expected to go from eight to one at this threshold. Again, we assume ε is small enough so that global solutions are not present.
To reduce (1.5b) to a discrete problem in the cubic region Ω = [−1, 1] 3 , a finite difference method is applied with uniform grid spacing h = 0.05. This relatively coarse grid cannot be expected to give accurate quantitative agreement and is mainly useful for observing the emergence of the multiple singularities. In Fig. 16, the numerical solution is visualized on three parallel planes with normal vectors (0, 1, 1). In the left panel, the solution for ε = 0.14 is shown for u ∞ = 5 × 10 2 and can be seen to concentrate on eight distinct points.
It seems natural to conjecture the following from the observed blow-up behavior in one, two and three dimensions: There exists an L c (n) such that for Ω = [−L, L] n and ψ = 0, problem (1.1b) exhibits blow-up at 2 n distinct points whenever L > L c (n)   Figure 15. Solution profiles of (1.5b) on the region (4.4) for f (u) = (1 + u) 2 , u ∞ = 1 × 10 3 and a range of ε. In panel (a), the solution profile is shown for the case where T ε < T S and blow-up occurs on ω(T ε ). The structure of ω(t) is visible as well as the undulations along it due to boundary curvature effects. In panels (b,e,f), solution profiles corresponding to ε values on segments I, II and III of Fig. 14(b) are shown. In panels (c,e), the two peak blow-up solutions very close to the critical ε values corresponding to the transition from I to II and II to III, are shown. It accounts for the absence of the multiple singularity phenomena in (6.2) through the solution profile in a boundary layer near ∂Ω. In the second order problem, the solution profile is monotone increasing towards its limiting value, whilst in the fourth order case there is an oscillatory approach and therefore a global max. There cannot be an overshoot in the second order problem as there is in the fourth order case, because the former has a maximum principle.
In addition, the leading order asymptotic analysis predicts x c = max x∈Ω d(x, ∂Ω) to be the point(s) favored by the dynamics of the PDE for singularity. Underpinning this prediction is the intuition that if u = 0 on ∂Ω, and u is monotone increasing away from ∂Ω, then points furthest from ∂Ω will have larger value and are consequently more likely to be selected for singularity by the dynamics of the PDE. To the author's knowledge, such results have previously been established for radially symmetric regions only [9].
The phenomenon described in the present work fits into a family of very interesting and unexpected solution behaviors associated with higher order PDEs. This behavior is particularly apparent when contrasted against the well understood second order case.
There are many interesting avenues of future work which can potentially emanate from this work. First, a rigorous justification of the phenomena described herein is highly desirable -it is possible that the order-preserving majorizing equation developed in [18] can provide a starting point for the analysis.
Second, in some degenerate examples, for example a stadium composed of a rectangle with semi-circular end pieces (cf. [8]), the predictive power of the asymptotic theory can be reduced by the jump in the curvature along ∂Ω. It would therefore be very interesting to study the blow-up set of (1.1b) as some regularity conditions on the curvature κ of ∂Ω are relaxed. In addition, how robust is the predictive accuracy of the asymptotic theory in such degenerate cases?